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

void plotweight(int isub){
    TFile *f = TFile::Open("/phenix/u/xuq/xuq/taxi/Run15pAu200FVTXClusAna503/8845/data/435527_0.root");
    int ivz = 0;
    TH1F* h = (TH1F*)f->Get(Form("phiweight_0_0_0_%d_0_%d",ivz,isub));
    TH1F* hc = (TH1F*)h->Clone("hc");
    h->Rebin(4);
    hc->Rebin(4);
    for(int ibin=0;ibin<=h->GetNbinsX();ibin++){
        if(h->GetBinContent(ibin)!=0)
        hc->SetBinContent(ibin,h->Integral()/h->GetNbinsX()/h->GetBinContent(ibin));
        else 
        hc->SetBinContent(ibin,1.);
    }
    TCanvas *c1 = new TCanvas("c1");
    hc->Draw();
    TCanvas *c2 = new TCanvas("c2");
    h->Draw();
}

mipana()
{
  ScanConfigFile();

  //  int nBins = (int)(((int)(MAX_BIN - MIN_BIN +1))/((int)BIN_SIZE));
  //  cout << "nBins = " << nBins  << endl;
  char inFileName[256];
  char outFileName[256];
  char tmpHistoName[256];

  //  TH2F *entries;


  for(int node = 0; node < nNodes; node  ++)
    {
  
      TH1F *tmpHist;
      sprintf(inFileName,"%s_ClusterEnergies.root", nodes[node]);
      sprintf(outFileName, "%s_FilteredClusterEnergies.root", nodes[node]);
      TFile *infile  = new TFile(inFileName, "read");
      TFile *outFile = new TFile(outFileName, "recreate"); 

      for(int z =Z_MIN; z < Z_MAX; z ++  )
	{
	  for(int x = X_MIN; x < X_MAX; x ++)
	    {
	      sprintf(tmpHistoName, "ClusterEnergies3x3_2_%d_%d_1", z, x);
	  
	      tmpHist = (TH1F*)infile->Get(tmpHistoName);

	      if(tmpHist != 0)
		{
		  cout << "rebinning histogram " << tmpHistoName <<endl;
 
		  outFile->cd();
		  tmpHist->Rebin(4);
		  tmpHist->Write();
		  infile->cd();
		}
	      else
		{
		  //	      cout << "ERROR: histogram " << tmpHistoName << "  doesnt exists" <<endl; 
		}

	      //do something
	    }
	}
  
      infile->Close();
      outFile->Close();
    }
}

TH1* his( TH1* input, Int_t nbins, Double_t xlow, Double_t xhigh ) {
  if ( !input ) { return 0; }
  TString name( TString(input->GetName()) + "Clone" );
  TH1F* output = new TH1F( name, "", input->GetNbinsX(), xlow, xhigh );
  for ( Int_t bin = 1; bin < input->GetNbinsX()+1; ++bin ) {
    Double_t centre = input->GetBinLowEdge(bin) + input->GetBinWidth(bin)/2.;
    if ( centre > xlow && 
	 centre < xhigh && 
	 input->GetBinContent(bin) > 0. ) {
      output->Fill( centre, input->GetBinContent(bin) );
    }
  }
  output->Sumw2();
  output->Rebin(Int_t(input->GetNbinsX()/nbins));
  return output;
}

void BackgroundFit() {
    TFile* inputFile = new TFile("0_MuScleFit.root", "READ");
    TH1F* histo = (TH1F*)inputFile->FindObjectAny("hRecBestRes_Mass");
    histo->Rebin(30);
    histo->Scale(1/histo->GetEntries());

    // TF1 * functionToFit = lorentzianFit();
    // TF1 * functionToFit = crystalBallFit();
    // TF1 * functionToFit = powerLawFit();
    // TF1 * functionToFit = lorentzianAndPowerLaw();
    // TF1 * functionToFit = exponentialFit();
    TF1 * functionToFit = lorenzianAndExponentialFit();

    histo->Fit(functionToFit, "M", "", 42, 160);

}

TH1F* readHist(TString nameHist,TString nameFile, TString nameDir, int rebin)
{
 TFile* file = new TFile(nameFile);
 file->ls();
 TDirectory *dir=(TDirectory*)file->Get(nameDir);
 
 TH1F* hist = (TH1F*)dir->Get(nameHist);
 // if (hist==0) return;
 hist->GetSumw2();
 // hist->SetLineWidth(2);
 if(rebin>0) hist->Rebin(rebin);
 hist->GetXaxis()->SetTitleSize(.055);
 hist->GetYaxis()->SetTitleSize(.055);
 hist->GetXaxis()->SetLabelSize(.05);
 hist->GetYaxis()->SetLabelSize(.05);
 hist->SetStats(kFALSE);
 return hist;
}

void ProbsFitter()
{
  TFile * inputFile = new TFile("Sherpa_nocuts.root", "READ");
  TH1F * histo = (TH1F*)inputFile->FindObjectAny("HistAllZMass");
  histo->Scale(1/histo->GetEntries());
  histo->Rebin(6);
  TCanvas * canvas = new TCanvas("canvas", "canvas", 1000, 800);
  canvas->Divide(2,3);

  // gStyle->SetOptFit(1);

  TF1 * residuals1 = iterateFitter(histo, 0, 0, canvas);
//   TF1 * residuals2 = iterateFitter(histo, 1, residuals1, canvas);
//   iterateFitter(histo, 2, residuals2, canvas);

//   canvas->cd(3);
//   TH1F * hclone = (TH1F*)histo->Clone();
//   TF1 * combinedFunctionToFit = combinedFit(residualFitFunction->GetParameters());
//   hclone->Fit(combinedFunctionToFit, "MN", "", 60, 120);
//   hclone->Draw();
//   combinedFunctionToFit->Draw("same");
//   combinedFunctionToFit->SetLineColor(kRed);

}

//.........这里部分代码省略.........
  hm4->Add(pm4t);
  hm4->Add(pm4s);

  THStack *hm5 = new THStack("hm5","JetsPt");
  hm5->Add(pm5q);
  hm5->Add(pm5t);
  hm5->Add(pm5s);


  //************************


  THStack *nj3t = new THStack("nj3t","3 j Evs tight Jets");
  nj3t->Add(njt3q);
  nj3t->Add(njt3t);
  nj3t->Add(njt3s);

  THStack *nj3l = new THStack("nj3l","3 j Evs loose Jets");
  nj3l->Add(njl3q);
  nj3l->Add(njl3t);
  nj3l->Add(njl3s);

  THStack *nj4t = new THStack("nj4t","4 j Evs tight Jets");
  nj4t->Add(njt4q);
  nj4t->Add(njt4t);
  nj4t->Add(njt4s);

  THStack *nj4l = new THStack("nj4l","4 j Evs loose Jets");
  nj4l->Add(njl4q);
  nj4l->Add(njl4t);
  nj4l->Add(njl4s);

  //Rebinning the mass plots
  Lhotmassq->Rebin(2);
  Lhotmasst->Rebin(2);
  Lhotmasss->Rebin(2);
  Thotmassq->Rebin(2);
  Thotmasst->Rebin(2);
  Thotmasss->Rebin(2);



  THStack *nj4l2m = new THStack("nj4l2m","4 j mass Evs loose Jets");
  nj4l2m->Add(Lhotmassq);
  nj4l2m->Add(Lhotmasst);
  nj4l2m->Add(Lhotmasss);

  THStack *nj4T2m = new THStack("nj4T2m","4 j mass Evs Tight Jets");
  nj4T2m->Add(Thotmassq);
  nj4T2m->Add(Thotmasst);
  nj4T2m->Add(Thotmasss);


  
 


  TCanvas *c1 = new TCanvas("c1","step1",1200,600);
  c1->Divide(2,1); 
  c1->cd(1);
  hpt->Draw();
  hpt->GetXaxis()->SetTitle("Jet Pt GeVc^{-1}");
  
  c1->cd(2);
  hm1->Draw();
  hm1->GetXaxis()->SetTitle("Jet Prunned Mass GeVc^{-2}"); 

void plot_histograms()
{
   gROOT->Reset();
   gROOT->SetStyle("Plain");
   gROOT->LoadMacro("./get_th1f.C");

   TH1F* hecore;
   TH1F* hecores;
   TH1F* hecore_first;
   TH1F* hecore_second;
   TH1F* hecore_seconds;
   TH1F* hecore_secondt;
   TH1F* hecore_secondts;

   TFile* finput = new TFile("../../data/histograms/h.sii.ln2o.root");
   //TFile* finput = new TFile("histograms/h.susie.k1.root");
   //TFile* finput = new TFile("histograms/h.si.k1.root");
   //TFile* finput = new TFile("histograms/h.siii.g2.root");
   //   TFile* finput = new TFile("histograms/h.siii.feb2010.root");
   finput->cd();
   hecore = (TH1F*) get_th1f("hecore")->Clone();
   hecores = (TH1F*) get_th1f("hecores")->Clone();
   hecore_first = (TH1F*) get_th1f("hecore_first")->Clone();
   hecore_second = (TH1F*) get_th1f("hecore_second")->Clone();
   hecore_seconds = (TH1F*) get_th1f("hecore_seconds")->Clone();
   hecore_secondt = (TH1F*) get_th1f("hecore_secondt")->Clone();
   hecore_secondts = (TH1F*) get_th1f("hecore_secondts")->Clone();
//   hecore->Rebin(10);

   hecore->Rebin(10);
   hecores->Rebin(10);

   hecore_first->Rebin(100);
   hecore_second->Rebin(100);
   hecore_seconds->Rebin(100);
   hecore_secondt->Rebin(100);
   hecore_secondts->Rebin(100);

   TCanvas* c1 = new TCanvas("c1","c1",10,10,900,600);
   c1->cd(); gPad->SetLogy();
   hecore->GetXaxis()->SetTitle("Ecore [keV]");
   hecore->GetYaxis()->SetTitle("Entries/10keV");
   hecore->Draw();
   hecores->SetLineColor(kRed);
   //hecores->SetFillColor(kRed);
   hecores->Draw("same");
   TLegend* lone = new TLegend(0.5,0.5,0.85,0.75);
   lone->SetFillStyle(0);
   lone->AddEntry(hecore,"all events", "l");
   lone->AddEntry(hecores,"single-segment events", "l");
   lone->Draw();
   c1->Update();
   c1->Print("../../data/plots/siio/ecore_ss.eps");


   TCanvas* c2 = new TCanvas("c2","c2",10,10,900,600);
   c2->cd();
   hecore_first->Draw();
   c2->Update();

   TCanvas* c3 = new TCanvas("c3","c3",10,10,900,600);
   c3->cd(); gPad->SetLogy();
   hecore_second->SetTitle("second event");
   hecore_second->GetXaxis()->SetTitle("Ecore [keV]");
   hecore_second->GetYaxis()->SetTitle("Entries/100keV");
   hecore_second->SetLineWidth(2);
   hecore_second->SetMinimum(0.1);
   hecore_second->Draw();
   //hecore_seconds->SetLineColor(kRed);
   //hecore_seconds->SetFillColor(kRed);
   //hecore_seconds->Draw("same");
   hecore_secondt->SetLineColor(kBlue);
   hecore_secondt->SetLineWidth(2);
   hecore_secondt->Draw("same");
   hecore_secondts->SetLineColor(kRed);
   hecore_secondts->SetLineWidth(2);
   hecore_secondts->Draw("same");
   TLegend* ltwo = new TLegend(0.4,0.5,0.85,0.75);
   ltwo->SetFillStyle(0);
   ltwo->AddEntry(hecore_second,"Bi214 consecutive events","l");
   ltwo->AddEntry(hecore_secondt,"5x164 #mus time window","l");
   ltwo->AddEntry(hecore_secondts,"5x164 #mus and single-segment","l");
   ltwo->Draw();
   c3->Update();
   c3->Print("../../data/plots/siio/ecore_bi214.eps");

}

//================================================
void DeltaZVsPos(const Int_t save = 0)
{
  THnSparseF *hn = (THnSparseF*)f->Get(Form("mhTrkDzDy_%s",trigName[kTrigType]));
  TList *list = new TList;

  // dz vs BL
  TH2F *hTrkDzVsBL = (TH2F*)hn->Projection(1,3);
  c = draw2D(hTrkDzVsBL,Form("%s: #Deltaz of matched track-hit pairs",trigName[kTrigType]));
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_BL_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_BL_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  list->Clear();
  TString legName[30];
  TH1F *hTrkDzInBL[30];
  Int_t counter = 0;
  for(Int_t i=0; i<30; i++)
    {
      hTrkDzInBL[i] = (TH1F*)hTrkDzVsBL->ProjectionY(Form("hDeltaZ_BL%d",i+1),i+1,i+1);
      if(hTrkDzInBL[i]->GetEntries()>0)
	{
	  legName[counter] = Form("Module %d",i+1);
	  hTrkDzInBL[i]->SetLineColor(color[counter]);
	  list->Add(hTrkDzInBL[i]);
	  counter ++;
	}
    }
  c = drawHistos(list,"TrkDzInBL",Form("%s: #Deltaz of matched track-hit pairs in backleg;#Deltaz (cm)",trigName[kTrigType]),kTRUE,-100,100,kTRUE,0,1.2*hTrkDzInBL[1]->GetMaximum(),kFALSE,kTRUE,legName,kTRUE,"",0.15,0.25,0.2,0.88,kFALSE,0.04,0.04,kFALSE,1,kTRUE,kFALSE);
  TLine *line = GetLine(0,0,0,1.1*hTrkDzInBL[1]->GetMaximum(),1);
  line->Draw();
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_BL_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_BL_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  // dz vs Mod
  TH2F *hTrkDzVsMod = (TH2F*)hn->Projection(1,4);
  c = draw2D(hTrkDzVsMod,Form("%s: #Deltaz of matched track-hit pairs",trigName[kTrigType]));
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_Mod_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_Mod_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  TH1F *hMthMod = (TH1F*)hTrkDzVsMod->ProjectionX("hMthMod");
  hMthMod->Sumw2();
  hMthMod->Scale(1./hMthMod->Integral());
  TH2F *hMtdHitMap = (TH2F*)f->Get(Form("mhMtdHitMap_%s",trigName[kTrigType]));
  TH1F *htmp = (TH1F*)hMtdHitMap->ProjectionY("hHitMod_finebin");
  htmp->Rebin(12);
  TH1F *hMtdHitMod = new TH1F(Form("hMtdHitMod_%s",trigName[kTrigType]),"# of MTD hits per module;module",5,1,6);
  for(int i=0; i<hMtdHitMod->GetNbinsX(); i++)
    {
      hMtdHitMod->SetBinContent(i+1,htmp->GetBinContent(i+1));
      hMtdHitMod->SetBinError(i+1,htmp->GetBinError(i+1));
    }
  hMtdHitMod->Scale(1./hMtdHitMod->Integral());
  list->Clear();
  list->Add(hMthMod);
  list->Add(hMtdHitMod);
  TString legName3[2] = {"Matched good hits","All good hits"};
  c = drawHistos(list,"MtdHitMod",Form("%s: MTD hits per module;module;probability",trigName[kTrigType]),kFALSE,0,5,kTRUE,0,0.5,kFALSE,kTRUE,legName3,kTRUE,"",0.15,0.25,0.6,0.88,kTRUE);
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sCompMtdHitMod_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sCompMtdHitMod_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  list->Clear();
  TString legName2[5];
  TH1F *hTrkDzInMod[5];
  for(Int_t i=0; i<5; i++)
    {
      hTrkDzInMod[i] = (TH1F*)hTrkDzVsMod->ProjectionY(Form("hDeltaZ_Mod%d",i+1),i+1,i+1);
      legName2[i] = Form("Module %d",i+1);
      list->Add(hTrkDzInMod[i]);
    }
  c = drawHistos(list,"TrkDzInMod",Form("%s: #Deltaz of matched track-hit pairs in module;#Deltaz (cm)",trigName[kTrigType]),kTRUE,-100,100,kTRUE,0,1.2*hTrkDzInMod[3]->GetMaximum(),kFALSE,kTRUE,legName2,kTRUE,"",0.15,0.25,0.6,0.88,kTRUE);
  TLine *line = GetLine(0,0,0,hTrkDzInMod[3]->GetMaximum()*1.05,1);
  line->Draw();
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_Mod_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_Mod_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }
}

void plottingmacro()
{
 setTDRStyle();
 gROOT->ForceStyle();
 initOptions();

 std::vector<Sample> s = samples();
 Sample data(1,"fake data","S1.root",0,true,1000);

 for(size_t i=0;i< s.size();i++) if(s[i].data) {data=s[i];break;}
 data.file()->ls(); 
 for(size_t i=0;i< s.size();i++) s[i].dump(data.lumi());

 std::vector<std::string> names;

 TList * subs = data.file()->GetListOfKeys();
 for(size_t i=0;i< subs->GetSize();i++)
  {
    TList * objs = ((TDirectoryFile *) data.file()->Get(subs->At(i)->GetName()))->GetListOfKeys();
     for(size_t j=0;j< objs->GetSize();j++)
     {
         names.push_back(subs->At(i)->GetName()+std::string("/")  + objs->At(j)->GetName());
 //      std::cout << subs->At(i)->GetName() << "/"  << objs->At(j)->GetName() << std::endl;
         //TODO: select plots via regexp
     }
    
  }

 for(size_t i = 0 ; i < names.size() ; i++) 
  {
   std::map<std::string,TH1F *> grouped;
   TString n=names[i];
   if(!n.Contains(TRegexp("V.*RegionH.*mu.*HiggsMass"))) continue;
   TCanvas *c = new TCanvas();
   c->SetLogy(true);
   c->SetTitle(names[i].c_str());
   TH1F *hd = ((TH1F*)data.file()->Get(names[i].c_str()));
   Options o=options[names[i]];
   hd->Rebin(o.rebin);
   hd->SetMarkerStyle(21);
   hd->Draw("E1");
   hd->SetYTitle(o.yaxis.c_str());
   THStack * sta = new THStack("sta",hd->GetTitle());
   TLegend * l = new TLegend(o.legendx1,o.legendy1,o.legendx2,o.legendy2); //0.7,0.1,0.9,0.6);
  
   l->AddEntry(hd, "Data","LP");

   for(size_t j=0;j< s.size() ;j++) 
   { 
       if(!s[j].data) 
      {
       TH1F * h = ((TH1F*)s[j].file()->Get(names[i].c_str()));
       h->Scale(s[j].scale(data.lumi()));
       h->SetLineColor(s[j].color);
       h->SetFillColor(s[j].color);
       h->Rebin(options[names[i]].rebin);
       if(grouped.find(s[j].name)==grouped.end()) {
          grouped[s[j].name]=(TH1F *)h->Clone(("_"+names[i]).c_str());
          l->AddEntry(h,s[j].name.c_str(),"F");
       }
       else
       {
        grouped[s[j].name]->Add(h);
       }
       sta->Add(h);
//     h->Draw("same");
      }
   }
   sta->Draw("same");
   hd->Draw("E1same");
   hd->GetYaxis()->SetRangeUser(options[names[i]].min,options[names[i]].max);
   l->Draw();


   std::cout << names[i] << " d: " <<  hd->Integral() << " ";
   THStack * sta2 = new THStack("sta2",hd->GetTitle());
   float tot=0;
   float toterr2=0;

   for(std::map<std::string,TH1F *>::iterator it = grouped.begin(); it != grouped.end();it++)
   {
             std::cout << it->first << " " << it->second->Integral() << " | " ;
             if(it->second->GetEntries() > 0) {
             float er=1.*sqrt(it->second->GetEntries())/it->second->GetEntries()*it->second->Integral();
             toterr2+=er*er;
             }
	     tot+=it->second->Integral();
             sta2->Add(it->second);
   }   
    std::cout << " Tot: " << tot << "+-" << sqrt(toterr2) <<  " SF: " << hd->Integral()/tot << std::endl;
    c = new TCanvas();
    sta2->Draw();
    hd->Draw("E1,same");
    sta2->GetYaxis()->SetRangeUser(options[names[i]].min,options[names[i]].max);
    hd->GetYaxis()->SetRangeUser(options[names[i]].min,options[names[i]].max);
    l->Draw();

   

  }
//.........这里部分代码省略.........

void singlejeteff(){


Double_t xbins[13] = {10,20,30,35,40,45,50,60,80,90,100,150,200};
int nbins = 12;

 TFile* g = new TFile("oct9/effoct11.root");
        g->TFile::Open();
          TH1F *h   = (TH1F*)g->Get("Denom");
	  h->Sumw2();
          h->Rebin(2,"hnew");




          TH1F *h20   = (TH1F*)g->Get("a4tt/hista4tt15");
          //h20->Sumw2();  h20->Rebin(nbins,"hnew20",xbins); 
          h20->Rebin(2,"hnew20");


          TH1F *h20h   = (TH1F*)g->Get("a4ttjes/hista4ttjes35");
          //h20h->Sumw2();  h20h->Rebin(nbins,"hnew20h",xbins); 
          h20h->Rebin(2,"hnew20h");

    
     
      
       TGraphAsymmErrors *h5=new TGraphAsymmErrors(hnew20,hnew);
         h5->Draw("AP");
         h5->SetMarkerStyle(20); //h5->SetMarkerSize(1.6);
         h5->SetMaximum(1.2); //to set the y axis maximum
         h5->SetMinimum(0.1); // to set the y axis minimum
         h5->GetXaxis()->SetRangeUser(0.,120.);
         h5->GetXaxis()->SetTitle("Leading jet p_{T,offline} (EM+JES) [GeV]");
         h5->GetYaxis()->SetTitle("Efficiency");


     
        TGraphAsymmErrors *h6=new TGraphAsymmErrors(hnew20h,hnew);
        h6->Draw("P,same");
        h6->SetMarkerStyle(22);h6->SetMarkerColor(2);h6->SetLineColor(2);h6->SetLineStyle(2);   



       

TPaveText *pt = new TPaveText(0.506182,0.4016529,0.75163,0.65,"blNDC");
   pt->SetBorderSize(0);
   pt->SetFillColor(0);
   pt->SetTextFont(42);
   pt->SetTextSize(0.05);            
   pt->AddText("| #eta_{offline} | <2.8");
   pt->Draw("same");




  TLegend *leg = new TLegend(0.506182,0.2016529,0.75163,0.4,"","brNDC");leg->SetFillColor(kWhite);  leg->SetBorderSize(0);
  leg->SetTextSize(0.05);
  leg->SetTextFont(42);
  leg->AddEntry(h5,"L2FS EM J15");    
  leg->AddEntry(h6,"L2FS EM+JES J35");  
  leg->SetBorderSize(0);
  leg->SetTextSize(0.05);   
  leg->SetTextFont(42);   
  leg->Draw("same");

         //c1->SetGridx();         c1->SetGridy();
                        //  c1->SetGridx();         c1->SetGridy();         //c1->SetLogx(1);
 //h6->Draw("HIST,C,same");h6->SetLineColor(2); 

       //c1->SetGridx();         c1->SetGridy();
  

          //TH1F *h5 = new TH1F("h5","L2FS_6j25  Efficiency",nbins,xbins);
         // TH1F *h5 = new TH1F("h5","L2FS_6j25  Efficiency",125,0,500);
         // h5->Divide(hnew20,hnew,1,1,"B");

}

void plotFR_QCD(){

  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetStatFormat("6.4f");
  gStyle->SetFitFormat("6.4f");
  int BoxValue = 11111111; //4680;  
  gStyle->SetOptFit(11);
  gStyle->SetOptDate(0);
  gStyle->SetOptTitle(0);
  //gStyle->SetOptStat(BoxValue);
  gStyle->SetOptStat(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetCanvasColor(0); //(10);
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadBottomMargin(0.15);
  gStyle->SetPalette(0);
  TPaveLabel pl;
  TLatex lt;
  lt.SetTextFont(70);
  lt.SetTextAlign(12);
  lt.SetTextSize(0.07);
  lt.SetTextColor(1);

  TPaveText* tText1 = new TPaveText(0.70, 0.90, 0.90, 0.95, "brNDC");
  tText1->SetBorderSize(0);
  tText1->SetFillColor(0);
  tText1->SetFillStyle(0);
  TText *t1 = tText1->AddText("(13 TeV)");
  tText1->SetTextSize(0.035);
  //tText1->Draw(); 
  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  TPaveText* tText2 = new TPaveText(0.2, 0.90, 0.4, 0.95, "brNDC");
  tText2->SetBorderSize(0);
  tText2->SetFillColor(0);
  tText2->SetFillStyle(0);
  TText *t2 = tText2->AddText("CMS Spring15 Simulation");
  tText2->SetTextSize(0.035);
  //tText2->Draw(); 


  TLegend *leg = new TLegend(0.17, 0.77, 0.4, 0.85);
  leg->SetBorderSize(0.0);
  leg->SetMargin(0.3);
  leg->SetFillColor(0);
  leg->SetFillStyle(10);
  leg->SetLineColor(0);
  Float_t tsize2 = 0.03;
  leg->SetTextSize(tsize2); 

  
  TFile *_file0 = TFile::Open("ZprimetoMuMu-MC-CMSSW745_FR_QCD_25ns.root");
  

  TH1F *num; _file0->GetObject("h1_Num_Pt_w",num);
  TH1F *den; _file0->GetObject("h1_Den_Pt_w",den);

  TH1F *numMB; _file0->GetObject("h1_Num_Pt_Barrel_w",numMB);
  TH1F *denMB; _file0->GetObject("h1_Den_Pt_Barrel_w",denMB);

  TH1F *numME; _file0->GetObject("h1_Num_Pt_EndCap_w",numME);
  TH1F *denME; _file0->GetObject("h1_Den_Pt_EndCap_w",denME);

  
  double Nbins = num->GetNbinsX();
  cout << "Nbins=" << Nbins << endl;

  int nRebin=10;

  num->Rebin(nRebin);
  den->Rebin(nRebin);

  numMB->Rebin(nRebin);
  denMB->Rebin(nRebin);

  numME->Rebin(nRebin);
  denME->Rebin(nRebin);


  // 
  TH1F *fake=new TH1F("fake","fake",Nbins/nRebin,0.,2000.);
  TH1F *fakeMB=new TH1F("fakeMB","fakeMB",Nbins/nRebin,0.,2000.);
  TH1F *fakeME=new TH1F("fakeME","fakeME",Nbins/nRebin,0.,2000.);
   
  cout << "Bins=" << Nbins/nRebin << endl;

  int* arraysize = new int[1];
  arraysize[0]=Nbins/nRebin;
  Float_t xMB[arraysize[0]],yMB[arraysize[0]],exlMB[arraysize[0]],exhMB[arraysize[0]],eylMB[arraysize[0]],eyhMB[arraysize[0]];


  for (unsigned int i=1; i<=Nbins/nRebin;i++){
    // All Muon
    if (denMB->GetBinContent(i)>0.) {
      fakeMB->SetBinContent(i,double(numMB->GetBinContent(i)/denMB->GetBinContent(i)));

      xMB[i-1]=fakeMB->GetBinCenter(i);
      yMB[i-1]=fakeMB->GetBinContent(i);
      exlMB[i-1]=0.;
      exhMB[i-1]=0.;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      if (strcmp(pu, "0PU") == 0)  pileup = "NoPileUp";
      if (strcmp(pu, "50PU") == 0) pileup = "50PileUp";
      if (strcmp(pu, "140PU") == 0) pileup = "140PileUp";
      std::cout << " pileup " << pileup <<" dec " << det << std::endl;

      HTSample *HT = new HTSample("LPC_MHTCUT/", "Wlv*_14TEV_HT", pileup, det);
      HT->InitSample(3000*1000);
      c1->cd();
      c1->Update();

      TH1F* reco = NULL;
      TH1F* gen = NULL;
      //TH2D* reco = NULL;
      //TH2D* gen = NULL;


      std::cout << " " << vLep.at(j) << std::endl;
      if (vLep.at(j) != "MET" && vLep.at(j) != "MHT")
      {
        TString reconame = "AppMatched"+vLep.at(j);
        TString genname = "AppGen"+vLep.at(j);
        //reco = (TH2D*)HT->GetTH2D(reconame.Data(), cut);
        //gen = (TH2D*)HT->GetTH2D(genname.Data(), cut);
        reco = (TH1F*)HT->GetTH1(reconame.Data(), cut);
        gen = (TH1F*)HT->GetTH1(genname.Data(), cut);
      } else {
        //reco = (TH1F*)HT->GetTH1(vLep.at(j), cut);
        //gen = (TH1F*)HT->GetTH1("GenMet");
      }


      if (vLep.at(j).find("Pt") != std::string::npos)
      {
        reco->Rebin(5);
        gen->Rebin(5);
        
      }
      //TH1F* reco = (TH1F*)HT->GetTH1("JetEta", cut);
      //TH1F* gen = (TH1F*)HT->GetTH1("GenJetEta", cut);
      //TH1F* reco = (TH1F*)HT->GetTH1("MuonPt", cut);
      //reco->Rebin(5);
      //TH1F* gen = (TH1F*)HT->GetTH1("GenMuonPt", cut);
      //gen->Rebin(5);
      //TH1F* reco = (TH1F*)HT->GetTH1("ElePt", cut);
      //
      //TH1F* gen = (TH1F*)HT->GetTH1("GenElePt", cut);

      reco->Divide(gen);
      reco->SetTitle("");
      //gen->Draw();
      //gen->SetLineColor(1);
      reco->SetLineWidth(3);
      reco->SetLineColor(VColor.at(i));

      ymin = ymin < reco->GetMinimum() ? ymin : reco->GetMinimum();
      ymax = ymax > reco->GetMaximum() ? ymax : reco->GetMaximum();
      if (i == 0)
      {
        yaxis = reco->GetYaxis();
        xaxis = reco->GetXaxis();
        reco->Draw();
      }
      else
      {
        reco->Draw("same");
      }

//.........这里部分代码省略.........
		addedhistos.AddAt(hthis,i-1);
		nbkg++;
		//cout << "Substituing bkg " << i << " + " << i-1 << " in addedhistos["<< i-1 <<"]" << endl;
	    }	    
	} // end of: if adding histograms
    }
    cout << " nbkg = " << nbkg << endl; 
    
    // Rebin histos if necessary, but first calculate KS:
    TH1F *hbkg = (TH1F*) addedhistos[nbkg];
    double KS = h0->KolmogorovTest(hbkg);
    double chi2ndf = h0->Chi2Test(hbkg, "UWUFOFCHI2/NDF");
    //cout << title.Data() << " KS = " << KS << " chi2/NDF = " << chi2ndf << endl;
    // Rebin? Set nrebin = 0 to NOT do rebinning. 
    // Will rebin only histos whose maximum x axis value exceeds 20. 
    // Anything with less will most probably be already made of integers, so no
    // need to rebin that! 
    Int_t nbinsx = h0->GetXaxis()->GetNbins();
    Int_t nbinsy = 100;
    Int_t nrebin = 5;
    if ( nbinsx > 750 && nbinsx <= 1000) nrebin = 30;
    if ( nbinsx > 400 && nbinsx <= 750 ) nrebin = 25;//20
    if ( nbinsx > 300 && nbinsx <= 400 ) nrebin = 25;//15
    if ( nbinsx > 200 && nbinsx <= 300 ) nrebin = 25;//15
    if ( nbinsx > 150 && nbinsx <= 200 ) nrebin = 10;//10
    if ( nbinsx > 100 && nbinsx <= 150 ) nrebin = 10;//10
    if ( nbinsx > 50 && nbinsx <= 100 )  nrebin = 10;//10
    if ( nbinsx > 20 && nbinsx <= 50 )   nrebin = 2;
    if ( nbinsx <= 20 ) nrebin = 1;  
    
    printf(" Saw nbins =%4i, rebinning by nrebin =%2i to final %3i bins \n",nbinsx,nrebin,int(nbinsx/nrebin));	

    if ( nrebin != 0 ) {
	h0->Rebin(nrebin); // data
	for (Int_t i = 0; i<=nbkg; i++){
	    TH1F * h = (TH1F*) addedhistos[i];
	    h->Rebin(nrebin);
	}
	for (Int_t i = 0; i<nsig; i++){
	    TH1F * h = (TH1F*) signalhistos[i];
	    h->Rebin(nrebin);
	}
    }

    // default text size: 0.045
    // make it bigger for the paper
    float textSize = 0.045;
    if(params->GetValue("Histo.Preliminary","yes")==TString("paper")) textSize=0.07;
    if(params->Defined("Histo.TextSize")) textSize=params->GetValue("Histo.TextSize",0.07);
    
    // Now, check largest dimensions so that we can plot all histograms at once.
    Float_t xmin=9999., xmax=-9999., ymin=9999., ymax=-9999.;
    for(Int_t i = 0; i<=nbkg; i++){
	TH1F * h = (TH1F*) addedhistos[i];
	ostringstream baseSrcName;
	baseSrcName << "Files." << i+1 << ".";
	TString bSrcName(baseSrcName.str().c_str());
	
	TAxis *axis = h->GetXaxis();
	if( axis->GetXmin() < xmin ) xmin = axis->GetXmin();
	if( axis->GetXmax() > xmax ) xmax = axis->GetXmax();
	if( h->GetMinimum() < ymin ) ymin = h->GetMinimum();
	if( h->GetMaximum() > ymax ) ymax = h->GetMaximum();
    }
    ymax = TMath::Nint(ymax*1.25+1); // Make enough room for the big legend
    TString title = h0->GetTitle();

//.........这里部分代码省略.........
  gsfGsfHisto.push_back("histoGsfGsfCorr");
  ////////////////////////////////////////////////////////////////////////////

  vector<TString> canvasName;
  canvasName.push_back("ratioFR");
  canvasName.push_back("ratioFRNoHeep");
  canvasName.push_back("ratioFRCorrDY");
  canvasName.push_back("ratioFRCorrFull");

  vector<TString> canvasTitle;
  canvasTitle.push_back("Fake rate ratio GSF-GSF / HEEP-GSF - uncorrected");
  canvasTitle.push_back("Fake rate ratio GSF-GSF / HEEP-GSF - non HEEP");
  canvasTitle.push_back("Fake rate ratio GSF-GSF / HEEP-GSF - DY corrected");
  canvasTitle.push_back("Fake rate ratio GSF-GSF / HEEP-GSF - full corrected");

  vector<TString> acroSuffix;
  acroSuffix.push_back("");
  acroSuffix.push_back("BB");
  acroSuffix.push_back("BE");
  acroSuffix.push_back("EE");

  vector<TString> suffix;
  suffix.push_back("");
  suffix.push_back(" EB-EB");
  suffix.push_back(" EB-EE");
  suffix.push_back(" EE-EE");

  sStream.str("");
  sStream << inputFilePrefix << lumi << "pb-1.root";
  TFile input(sStream.str().c_str(), "read");
  input.cd();

  cout << endl << "Input file: " << sStream.str() << endl;

  // to keep the histogram when the file is closed
  TH1::AddDirectory(kFALSE);
  TH1::SetDefaultSumw2(kTRUE);
  
  for (unsigned int j = 0; j < 4; ++j) {
    if (!plotClosureTest[j]) continue;
    for (unsigned int p = 0; p < 4; ++p) {
      if (!plotHisto[p]) continue;
  
      TCanvas *c0 = new TCanvas(canvasName[j] + acroSuffix[p], canvasTitle[j] + suffix[p], 100, 100, 800, 600);
      c0->cd();
      c0->SetBorderMode(0);
      c0->SetFrameBorderMode(0);
      c0->SetFillColor(0);
      c0->SetFrameFillColor(0);
      gStyle->SetOptStat(0);
      gStyle->SetPadTickY(1);
  
      // get the histograms 
      input.cd(folderGsfGsfHisto[j]);
      TH1F *numHisto = (TH1F *)gDirectory->Get(gsfGsfHisto[j] + acroSuffix[p]);
      input.cd(folderHeepGsfHisto[j]);
      TH1F *denomHisto = (TH1F *)gDirectory->Get(heepGsfHisto[j] + acroSuffix[p]);
  
      TH1F *numHistoRebinned = (TH1F *)numHisto->Rebin(nBins, "numHistoRebinned" + acroSuffix[p], binArray);
      TH1F *denomHistoRebinned = (TH1F *)denomHisto->Rebin(nBins, "denomHistoRebinned" + acroSuffix[p], binArray);
  
      TH1F *ratioHisto = new TH1F("histoRatioCorr" + acroSuffix[p], canvasTitle[j] + suffix[p], nBins, binArray);
      ratioHisto->Divide(numHistoRebinned, denomHistoRebinned);
  
      ratioHisto->SetLineColor(4);
      ratioHisto->SetMarkerColor(4);
      ratioHisto->SetMarkerStyle(20);
  
      ratioHisto->Draw();
      ratioHisto->Fit("pol0", "+", "lep", 120, 500);
  
      sStream.str("");
      sStream << "#sqrt{s} = 7TeV,  #int L dt = " << lumi << "pb^{-1}";
      TPaveLabel *label0 = new TPaveLabel(0.6, 0.7, 0.9, 0.8, sStream.str().c_str(), "brNDC");
      label0->SetFillColor(0);
      label0->SetFillStyle(0);
      label0->SetBorderSize(0);
      label0->SetTextSize(0.30);
      label0->SetTextFont(font);
      label0->Draw("sames");
      TPaveLabel *label1 = new TPaveLabel(0.6, 0.8, 0.9, 0.9, "CMS preliminary", "brNDC");
      label1->SetFillColor(0);
      label1->SetFillStyle(0);
      label1->SetBorderSize(0);
      label1->SetTextSize(0.40);
      label1->SetTextFont(font);
      label1->Draw("sames");
  
      //TLegend *legend = new TLegend(0.38, 0.6, 0.53, 0.9);
      //legend->SetTextSize(0.03);
      //legend->SetTextFont(font);
      //legend->SetBorderSize(0);
      //legend->SetFillStyle(0);
      //legend->AddEntry(ratioHisto, "Ratio GSF-GSF / HEEP-GSF (non HEEP)", "lep");
      //legend->Draw("sames");

    } // end loop over eta ranges
  } // end loop over corrections
  input.Close();
}