C++ TH1F::SetMaximum方法代码示例

void multifit() {
   const Int_t np = 49;
   Float_t x[np] = {1.913521, 1.953769, 2.347435, 2.883654, 3.493567,
                    4.047560, 4.337210, 4.364347, 4.563004, 5.054247,
                    5.194183, 5.380521, 5.303213, 5.384578, 5.563983,
                    5.728500, 5.685752, 5.080029, 4.251809, 3.372246,
                    2.207432, 1.227541, 0.8597788,0.8220503,0.8046592,
                    0.7684097,0.7469761,0.8019787,0.8362375,0.8744895,
                    0.9143721,0.9462768,0.9285364,0.8954604,0.8410891,
                    0.7853871,0.7100883,0.6938808,0.7363682,0.7032954,
                    0.6029015,0.5600163,0.7477068,1.188785, 1.938228,
                    2.602717, 3.472962, 4.465014, 5.177035};

   TH1F *h = new TH1F("h","Example of several fits in subranges",np,85,134);
   h->SetMaximum(7);

   for (int i=0;i<np;i++) {
      h->SetBinContent(i+1,x[i]);
   }

   Double_t par[9];
   TF1 *g1    = new TF1("g1","gaus",85,95);
   TF1 *g2    = new TF1("g2","gaus",98,108);
   TF1 *g3    = new TF1("g3","gaus",110,121);
   TF1 *total = new TF1("total","gaus(0)+gaus(3)+gaus(6)",85,125);
   total->SetLineColor(2);
   h->Fit(g1,"R");
   h->Fit(g2,"R+");
   h->Fit(g3,"R+");
   g1->GetParameters(&par[0]);
   g2->GetParameters(&par[3]);
   g3->GetParameters(&par[6]);
   total->SetParameters(par);
   h->Fit(total,"R+");
}

void Deconvolution_wide_boost() {
   Int_t i;
   const Int_t nbins = 256;
   Double_t xmin     = 0;
   Double_t xmax     = nbins;
   Double_t source[nbins];
   Double_t response[nbins];
   gROOT->ForceStyle();

   TH1F *h = new TH1F("h","Deconvolution",nbins,xmin,xmax);
   TH1F *d = new TH1F("d","",nbins,xmin,xmax);

   TString dir  = gROOT->GetTutorialDir();
   TString file = dir+"/spectrum/TSpectrum.root";
   TFile *f     = new TFile(file.Data());
   h = (TH1F*) f->Get("decon3");
   h->SetTitle("Deconvolution of closely positioned overlapping peaks using boosted Gold deconvolution method");
   d = (TH1F*) f->Get("decon_response_wide");

   for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
   for (i = 0; i < nbins; i++) response[i]=d->GetBinContent(i + 1);

   h->SetMaximum(200000);
   h->Draw("L");
   TSpectrum *s = new TSpectrum();
   s->Deconvolution(source,response,256,200,50,1.2);

   for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
   d->SetLineColor(kRed);
   d->Draw("SAME L");
}

void seism() {

   TStopwatch sw; sw.Start();
   //set time offset
   TDatime dtime;
   gStyle->SetTimeOffset(dtime.Convert());

   TCanvas *c1 = new TCanvas("c1","Time on axis",10,10,1000,500);
   c1->SetFillColor(42);
   c1->SetFrameFillColor(33);
   c1->SetGrid();

   Float_t bintime = 1; //one bin = 1 second. change it to set the time scale
   TH1F *ht = new TH1F("ht","The ROOT seism",10,0,10*bintime);
   Float_t signal = 1000;
   ht->SetMaximum( signal);
   ht->SetMinimum(-signal);
   ht->SetStats(0);
   ht->SetLineColor(2);
   ht->GetXaxis()->SetTimeDisplay(1);
   ht->GetYaxis()->SetNdivisions(520);
   ht->Draw();

   for (Int_t i=1;i<2300;i++) {
      //======= Build a signal : noisy damped sine ======
      Float_t noise  = gRandom->Gaus(0,120);
      if (i > 700) noise += signal*sin((i-700.)*6.28/30)*exp((700.-i)/300.);
      ht->SetBinContent(i,noise);
      c1->Modified();
      c1->Update();
      gSystem->ProcessEvents(); //canvas can be edited during the loop
   }
   printf("Real Time = %8.3fs, Cpu Time = %8.3fs\n",sw.RealTime(),sw.CpuTime());
}

void DrawInsituCorr_vs_Eta(Str jetAlgo) {
  Nbins=sizeof(Ebins)/sizeof(double);
  myJES = new JetCalibrationTool(jetAlgo,_jesFile);

  TH1F *h = new TH1F("","",100,-4.8,4.8);
  h->SetXTitle("Jet #eta_{det}");
  h->SetYTitle("Residual correction for data only");
  h->SetMaximum(1.1); h->SetMinimum(0.85);

  h->Draw();
  TLine *l = new TLine(); 
  l->SetLineWidth(2); l->SetLineColor(kGray+2); l->SetLineStyle(4);
  l->DrawLine(-4.65,1,4.65,1);
  
  for (int i=0;i<Nbins;++i) {
    Graph *g = GetInsituGraph(Ebins[i]);
    FormatGraph(g,i); g->Draw("PL");
    double x=0.52, y=0.3-0.05*i;
    if (i>2) { x=0.74; y=0.3-0.05*(i-3); }
    DrawLabel(Form("p_{T} = %.0f GeV",Ebins[i]),x,y,i);
  }

  l->SetLineWidth(2);
  l->SetLineColor(kGray+2);
  l->SetLineStyle(1);

  tex->SetNDC(); tex->SetTextAlign(12);
  //tex->DrawLatex(0.18,0.975,myJES->AbsoluteJES_Description());
  tex->DrawLatex(0.18,0.9,GetJetDesc(jetAlgo));
}

void DrawJMS_vs_Eta(Str jetAlgo) {
  double massEbins[] = {50, 100, 250, 500, 750, 1500};
  double massNbins = sizeof(massEbins)/sizeof(double);
  myJES = new JetCalibrationTool(jetAlgo,_jesFile);
  TH1F *h = new TH1F("","",100,-4.8,4.8);
  h->SetXTitle("Jet #eta_{det}"); h->SetYTitle("Jet mass response");
  h->SetMaximum(1.4); h->SetMinimum(0.6);
  h->Draw();
  
  TLine *l = new TLine(); 
  l->SetLineWidth(2); l->SetLineColor(kGray+2); l->SetLineStyle(4);
  l->DrawLine(-4.65,1,4.65,1);

  for (int i=0;i<massNbins;++i) {
    Graph *g = GetJMSGraph(massEbins[i]);
    FormatGraph(g,i);
    g->Draw("PL");
    double x=0.52, y=0.3-0.05*i;
    if (i>2) { x=0.74; y=0.3-0.05*(i-3); }
    DrawLabel(Form("E = %.0f GeV",massEbins[i]),x,y,i);
  }
  
  tex->SetNDC(); tex->SetTextAlign(12);
  tex->DrawLatex(0.18,0.975,myJES->AbsoluteJES_Description());
  tex->DrawLatex(0.18,0.9,jetAlgo);
}

void DrawInsituCorr_vs_Pt(Str jetAlgo) {
  myJES = new JetCalibrationTool(jetAlgo,_jesFile);
  TH1F *h = new TH1F("","",100,15,2500);
  h->SetXTitle("Jet p_{T}");
  h->SetYTitle("Residual correction for data only");
  h->SetMaximum(1.1); h->SetMinimum(0.85);
  h->Draw();
  TLine *l = new TLine(); 
  l->SetLineWidth(2); l->SetLineColor(kGray+2); l->SetLineStyle(4);
  l->DrawLine(16,1,2400,1);
  for (int i=0;i<6;++i) {
    double eta=4.0*i/5;
    Graph *g = GetInsituGraphVsPt(eta);
    FormatGraph(g,i); g->Draw("PL");
    double x=0.52, y=0.3-0.05*i;
    if (i>2) { x=0.74; y=0.3-0.05*(i-3); }
    DrawLabel(Form("#eta = %.1f",eta),x,y,i);
  }

  l->SetLineWidth(2);
  l->SetLineColor(kGray+2);
  l->SetLineStyle(1);

  tex->SetNDC(); tex->SetTextAlign(12);
  //tex->DrawLatex(0.18,0.975,myJES->AbsoluteJES_Description());
  tex->DrawLatex(0.18,0.9,GetJetDesc(jetAlgo));
}

TH1F* CreateEffHisto(TH1F* hGen, TH1F* hRec)
{
// create an efficiency histogram

  Int_t nBins = hGen->GetNbinsX();
  TH1F* hEff = (TH1F*) hGen->Clone("hEff");
  hEff->SetTitle("");
  hEff->SetStats(kFALSE);
  hEff->SetMinimum(0.);
  hEff->SetMaximum(110.);
  hEff->GetYaxis()->SetTitle("#epsilon [%]");

  for (Int_t iBin = 0; iBin <= nBins; iBin++) {
    Double_t nGen = hGen->GetBinContent(iBin);
    Double_t nRec = hRec->GetBinContent(iBin);
    if (nGen > 0) {
      Double_t eff = nRec/nGen;
      hEff->SetBinContent(iBin, 100. * eff);
      Double_t error = sqrt(eff*(1.-eff) / nGen);
      if (error == 0) error = 0.0001;
      hEff->SetBinError(iBin, 100. * error);
    } else {
      hEff->SetBinContent(iBin, -100.);
      hEff->SetBinError(iBin, 0);
    }
  }

  return hEff;
}

void PlotMinEtFromSim(Bool_t isPhos = kFALSE){
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  Float_t min = 0;
  float max = 1;
  TString filename, detname;
  if(isPhos){
    min = 0.655;
    max = 0.785;
    detname = "PHOS";
    filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.PHOS.LHC11a10a_bis.Run139465.root";
  }
  else{
    min = 0.58;
    max = 0.725;
    filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root";
    detname = "EMCal";
  }
  
  TFile *f = TFile::Open(filename, "READ");
  TList *l = dynamic_cast<TList*>(f->Get("out1"));
  TH1F *fHistSimulatedGammaEnergyAboveThreshold = l->FindObject("fHistSimulatedGammaEnergyAboveThreshold");
  TH1F *fHistSimulatedGammaEnergy = l->FindObject("fHistSimulatedGammaEnergy");
  SetStyles(fHistSimulatedGammaEnergyAboveThreshold,20,TColor::kRed);
  fHistSimulatedGammaEnergyAboveThreshold->Divide(fHistSimulatedGammaEnergy);

    TCanvas *c1 = new TCanvas("c1","Simulation",600,400);
    c1->SetTopMargin(0.02);
    c1->SetRightMargin(0.03);
    c1->SetLeftMargin(0.11745);
    c1->SetBottomMargin(0.11745);
    c1->SetBorderSize(0);
    c1->SetFillColor(0);
    c1->SetFillColor(0);
    c1->SetBorderMode(0);
    c1->SetFrameFillColor(0);
    c1->SetFrameBorderMode(0);
    fHistSimulatedGammaEnergyAboveThreshold->SetMaximum(max +0.1);
    fHistSimulatedGammaEnergyAboveThreshold->SetMinimum(min-0.1);
    fHistSimulatedGammaEnergyAboveThreshold->GetXaxis()->SetTitle("Centrality bin");
    fHistSimulatedGammaEnergyAboveThreshold->GetYaxis()->SetTitle("f_{minEt}");
    fHistSimulatedGammaEnergyAboveThreshold->GetYaxis()->SetLabelSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->GetXaxis()->SetLabelSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->GetYaxis()->SetTitleSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->GetXaxis()->SetTitleSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->Draw();
    TLine *lineMin = new TLine(-0.5,min,19.5,min);
    lineMin->Draw();
    TLine *lineMax = new TLine(-0.5,max,19.5,max);
    lineMax->Draw();
    lineMin->SetLineColor(TColor::kBlue);
    lineMax->SetLineColor(TColor::kBlue);
    lineMin->SetLineStyle(2);
    lineMax->SetLineStyle(2);

    TString outfile = "/tmp/MinEtFromSim"+detname+".png";
    c1->SaveAs(outfile.Data());
}

void draw_ratio(std::vector<TH1F*> h,
	  TString name, TString xTitle,
	  double xmin, double xmax,
	  TString legHeader = "", bool legRIGHT = true, bool legTOP = true,
	  bool logX = false, bool stat = false, int rebin = -1, int orbin = -1,
	  TString option = "", int nclus = 99) {  //double ymin_ratio, double ymax_ratio,
  TCanvas* can = new TCanvas(name+"_ratio",name+"_ratio",900,450);
  can->cd();

  double legxmin = (legRIGHT ? 0.55 : 0.18);
  double legxmax = legxmin+0.25;
  double legymin = (legTOP ? 0.70 : 0.15);
  double legymax = legymin+0.15;
  TLegend* leg = new TLegend(legxmin,legymin,legxmax,legymax);
  if (legHeader!="") leg->SetHeader(legHeader);
  leg->SetTextSize(0.04);
  leg->SetFillColor(0);
  leg->SetLineColor(0);

  TString options = (option=="" ? "pe" : option);

  if (rebin>0) h[h.size()-1]->Rebin(rebin); //to rebin benchmark before divide
  for (size_t i=0; i<h.size(); i++) {
    //if(h[i]->GetNbinsX() != orbin) cout << "WARNING: orbin for " << h[i]->GetName() << " are " << h[i]->GetNbinsX() << endl; //debug - shift of h[][] wrt clu[][]
    if (rebin>0 && i<(h.size()-1)) h[i]->Rebin(rebin);
    TH1F* ratio = (TH1F*)h[i]->Clone("ratio_"+name);
    ratio->Sumw2();
    ratio->Divide(h[h.size()-1]); //benchmark is at the end.

    if (logX) gPad->SetLogx();
    ratio->SetMarkerStyle(20+i);
    ratio->SetMarkerSize(1.0); //1.2
    ratio->GetXaxis()->SetRangeUser(xmin,xmax);
    ratio->SetMinimum(-0.1);
    ratio->SetMaximum(4);
    if (i==0){    //just for the first one
     ratio->GetXaxis()->SetLabelSize(0.05);
     ratio->GetXaxis()->SetTitle(xTitle);
     ratio->GetXaxis()->SetTitleOffset(1);
     ratio->GetXaxis()->SetTitleSize(0.06);
     ratio->GetYaxis()->SetTitle("ratio");
     ratio->GetYaxis()->SetTitleSize(0.06);
     ratio->GetYaxis()->SetTitleOffset(0.7);
     ratio->GetYaxis()->SetLabelSize(0.05);
    }
    string nam = "";
    nam = translate(clu[nclus][0].c_str());
    if(i==(h.size()-1)) leg->AddEntry(h[i],nam.c_str(),"l"); //to print only benchmark (first in the list)
    if (i==1) options = options + (stat ? "sames" : "same"); //once is enought
    ratio->Draw(options);
  }
  leg->Draw("same");
  drawPrivate(0.04);
  can->Update();
  can->SaveAs(Outfolder+name+"_ratio.png");
}

void Background_order() {
   Int_t i;
   const Int_t nbins = 4096;
   Double_t xmin     = 0;
   Double_t xmax     = 4096;
   Double_t source[nbins];
   gROOT->ForceStyle();

   TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
   TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax);
   TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax);
   TH1F *d4 = new TH1F("d4","",nbins,xmin,xmax);

   TString dir  = gROOT->GetTutorialsDir();
   TString file = dir+"/spectrum/TSpectrum.root";
   TFile *f     = new TFile(file.Data());
   TH1F *back = (TH1F*) f->Get("back2");
   back->SetTitle("Influence of clipping filter difference order on the estimated background");
   back->SetAxisRange(1220,1460);
   back->SetMaximum(3000);
   back->Draw("L");

   TSpectrum *s = new TSpectrum();

   for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
   s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
                 TSpectrum::kBackOrder2,kFALSE,
                 TSpectrum::kBackSmoothing3,kFALSE);
   for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]);
   d1->SetLineColor(kRed);
   d1->Draw("SAME L");

   for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
   s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
                 TSpectrum::kBackOrder4,kFALSE,
                 TSpectrum::kBackSmoothing3,kFALSE);
   for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]);
   d2->SetLineColor(kBlue);
   d2->Draw("SAME L");

   for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
   s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
                 TSpectrum::kBackOrder6,kFALSE,
                 TSpectrum::kBackSmoothing3,kFALSE);
   for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,source[i]);
   d3->SetLineColor(kGreen);
   d3->Draw("SAME L");

   for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
   s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
                 TSpectrum::kBackOrder8,kFALSE,
                 TSpectrum::kBackSmoothing3,kFALSE);
   for (i = 0; i < nbins; i++) d4->SetBinContent(i + 1,source[i]);
   d4->SetLineColor(kMagenta);
   d4->Draw("SAME L");
}

//------------------------------------//
// Main
//------------------------------------//
void AnaEField()
{

  // Make Canvas 
  TCanvas* can = makeCanvas("can");
  can->SetLogy();
  can->SetGridx();
  can->SetGridy();

  // Define the function here with parameters:
  // [0] = L
  // [1] = freq
  // [2] = k
  // [3] = n
  // [4] = sqrt(2pi) * mu * mu0
  TString func = "[4]*[0]*[1]*sin(x/57.2957)*exp(-pow([2]*[0],2)*pow(cos(x/57.2957)-1/[3],2)/2)";
  //TString func = "[4]*[1]*sin(x/57.2957)*exp(-pow([2]*[0],2)*pow(cos(x/57.2957)-1/[3],2)/2)";
  //TString func = "[4]*[0]*[1]*exp(-pow([2]*[0],2)*pow(cos(x/57.2957)-1/[3],2)/2)";

  // Now fix some of the variables:
  double freq = 1e9; //600e6; // Hz
  double c    = 3e8;   // m/s
  double n    = 1.3; //1.78;
  double k    = 2*TMath::Pi() * n * freq/c; // 1/m
  double Const= sqrt(2*TMath::Pi()) * 4 * TMath::Pi() * 1e-7 * 1;

  // Make a dummy histogram
  TH1F* h = makeHist("h",100,0,90,"Angle [deg]", "Field Strength",kBlack,0);
  h->Fill(1e10);
  h->SetMinimum(10);
  h->SetMaximum(3000);
  h->Draw();

  // Make two instances of the function
  TF1* f0 = new TF1("f0",func.Data(),0,90);
  setParam(f0,1.2,freq,k,n,Const,kBlack);

  TF1* f1 = new TF1("f1",func.Data(),0,90);
  setParam(f1,0.1,freq,k,n,Const,kBlue);
 
  // Draw
  f0->Draw("same");
  f1->Draw("same");
  
  // Add Legend
  TLegend* leg = makeLegend(0.15,0.3,0.7,0.9);
  leg->SetHeader("f = 600 MHz");
  leg->AddEntry(f0,"L = 1.2m","l");
  leg->AddEntry(f1,"L = 0.1m","l");
  leg->Draw("same");
  
  //can->SaveAs("../plots/quickAnalytic_600MHz.png");

}

TCanvas * plot (TH1F* histoDataIn, TString legendData, TH1F* histoSimulationIn, TString legendSimulation,
                TString & canvasName, Float_t maximum = 0.15, TString xAxisTitle = "#eta",
                TString yAxisTitle = "Number of Clusters", TString error = "", bool useLegend = true,
                TString text = "", Float_t textX = 0.7, Float_t textY = 0.4, Float_t rebin = 0 ) {

  TH1F * histoData = (TH1F*)histoDataIn->Clone();
  TH1F * histoSimulation = (TH1F*)histoSimulationIn->Clone();

  // histoData->Sumw2();
  histoData->Scale(1/(histoData->Integral()));
  histoSimulation->Scale(1/(histoSimulation->Integral()));

  // Create also the legend and add the histograms
  TLegend * legend = new TLegend( 0.55, 0.65, 0.76, 0.82 );
  legend->AddEntry( histoData, xAxisTitle );
  legend->AddEntry( histoSimulation, yAxisTitle, "F" );

  cout << "histoData = " << histoData << endl;
  cout << "histoSimulation = " << histoSimulation << endl;

  TCanvas * c = new TCanvas( canvasName, canvasName, 1000, 800 );
  c->Draw();

  histoSimulation->SetMaximum(maximum);

  histoSimulation->SetFillColor(kRed);
  // histoSimulation->SetLineWidth(0);
  histoSimulation->SetLineColor(kRed);
  histoSimulation->SetXTitle(xAxisTitle);
  histoSimulation->SetYTitle(yAxisTitle);
  histoSimulation->SetTitleOffset(1.6,"Y");
  histoSimulation->SetTitle();

  histoData->SetLineStyle(1);
  histoData->SetLineWidth(2.5);

  histoSimulation->Draw();
  histoData->Draw("same");

  legend->SetFillColor(kWhite);
  if (useLegend) legend->Draw("same");

  if ( text != "" ) {
    TPaveText * pt = new TPaveText(textX, textY, textX+0.2, textY+0.17, "NDC" ); // "NDC" option sets coords relative to pad dimensions
    pt->SetFillColor(0); // text is black on white
    pt->SetTextSize(0.08); 
    pt->SetBorderSize(0);
    pt->SetTextAlign(12);
    pt->AddText(text);
    pt->Draw("same");       //to draw your text object
  }

  return c;
};

CheckEnabledChannels(const Char_t *runlist)
{

  ifstream is(runlist);
  Char_t buf[4096];
  Int_t run[1024];
  Int_t nrun = 0;
  while(!is.eof()) {
    is.getline(buf, 4096);
    if (is.eof()) break;
    run[nrun] = atoi(buf);
    printf("added run number %d\n", run[nrun]);
    nrun++;
  }
  printf("%d runs added\n", nrun);
  is.close();

  TH1F *hActive = new TH1F("hActive", "active channels;run;fraction", nrun, 0, nrun);
  TH1F *hReadout = new TH1F("hReadout", "good readout;run;fraction", nrun, 0, nrun);
  for (Int_t irun = 0; irun < nrun; irun++) {
    hr = CheckEnabledChannels(run[irun], kTRUE);
    ha = CheckEnabledChannels(run[irun], kFALSE);
    hReadout->SetBinContent(irun + 1, hr->Integral());
    hActive->SetBinContent(irun + 1, ha->Integral());
    hReadout->GetXaxis()->SetBinLabel(irun + 1, Form("%d", run[irun]));
    delete hr; delete ha;
  }
  
  hReadout->SetMarkerStyle(20);
  hReadout->SetMarkerColor(4);
  hActive->SetMarkerStyle(25);
  hActive->SetMarkerColor(2);
  hReadout->Sumw2();
  hActive->Sumw2();
  hReadout->Divide(hReadout, hActive, 1., 1., "B");
  hActive->Scale(1. / 152928.);
  hReadout->SetMinimum(0.);
  hReadout->SetMaximum(1.);
  hReadout->Draw("E");
  hActive->Draw("E, same");
  TLegend *l = gPad->BuildLegend();
  l->SetFillStyle(0);

}

void CompareGammaPoisson(float N=100)
{ 

    gStyle->SetOptStat(111111110);
    gStyle->SetStatW(0.3);                
    gStyle->SetStatH(0.25);                

    TRandom1 *frand = new TRandom1();
    TRandom3 rand(1234); // gamma

    int Nbins = (int)(N+TMath::Sqrt(N)*10);

    TH1F *hpoisson    = new TH1F("hpoisson",    "Poisson",   Nbins, -0.5,  Nbins-0.5);
    TH1F *hgamma      = new TH1F("hgamma",      "Gamma",     Nbins, -0.5,  Nbins-0.5);
    hpoisson->Sumw2();
    hgamma->Sumw2();

    for(int i=0; i<Npseudo; i++)  
    { 
        int checkpoint = (int)Npseudo/10;
        //if((i%checkpoint)==0) cout << "Generated " << i << "/" << Npseudo << " experiments " << endl;
        
        float Npoisson      = frand->Poisson(N);
        float Ngamma        = gsl_ran_gamma(2*N+1,0.5,rand);;

        hpoisson->Fill(Npoisson);
        hgamma->Fill(Ngamma);

    }
    hpoisson->Scale(1./hpoisson->Integral());
    hgamma->Scale(1./hgamma->Integral());

    h1cosmetic(hpoisson, Form("Poisson(black) and Gamma(red), N=%i",(int)N), kBlack, 3, 0, "");
    h1cosmetic(hgamma,   "Gamma", kRed, 3, 0, "");

    TCanvas *c = new TCanvas("c","c",400,400);
    c->cd(1);
    hpoisson->SetMaximum(hpoisson->GetMaximum()>hgamma->GetMaximum()?hpoisson->GetMaximum()*1.2:hgamma->GetMaximum()*1.2);
    hpoisson->Draw("HIST");
    hgamma->Draw("HIST SAME");
    c->Print(Form("fig/ComparePoissonGamma_N%i.pdf", (int)N));
}

void drawEmptyHisto ( float Xmin, float Xmax, TString xTitle, float Ymin, float Ymax, TString yTitle, TString name )
{
    printf ( "Setting empty histo: X: %.2f-%.2f\tY: %.2f-%.2f\n",Xmin,Xmax,Ymin,Ymax );
    printf("Total luminosity to be plotted: %.3f\n",Xmax);
    printf("Luminosity scale used: %.2f\n",lumiScale);

    Xmax*=1.02;	  // Increasing max lumi to have free space
    gStyle->SetOptStat ( 0 );
    //    TGaxis::SetMaxDigits ( 4 );

    TH1F *histo = new TH1F ( name, xTitle, 1, Xmin, Xmax );
    TAxis *xAx = histo->GetXaxis();
    TAxis *yAx = histo->GetYaxis();
    xAx->SetTitle ( xTitle );
    //    xAx->SetTitleOffset ( 1.1 );
    //    xAx->SetNdivisions ( 510 );
    yAx->SetTitle ( yTitle );
    //    yAx->SetTitleOffset ( 1.15 );
    //    yAx->SetNdivisions ( 510 );
    histo->SetMinimum ( Ymin );
    histo->SetMaximum ( Ymax );
    histo->Draw ( "AXIS" );
}