C++ TH1D::GetBinLowEdge方法代码示例

std::pair<float,float> getFitRange(TH1D& hist) {
  int maxBin = hist.GetMaximumBin();
  float targetVal = hist.GetBinContent( maxBin ) /2.;


  int nBins = hist.GetNbinsX();
  
  std::pair<float,float> results = {hist.GetBinLowEdge(1),hist.GetBinLowEdge(nBins+1)};

  int diff=0;
  while( maxBin-(++diff) >0) {
    if( hist.GetBinContent(maxBin-diff) < targetVal ){
      results.first = hist.GetBinCenter(maxBin-diff);
      break;
    }
  }

  diff=0;
  while( maxBin+(++diff) <=nBins) {
    if( hist.GetBinContent(maxBin+diff) < targetVal ){
      results.second = hist.GetBinCenter(maxBin+diff);
      break;
    }
  }
  
  
  return results;
}

void loglikdistrib(Int_t ntrials = 10000, Bool_t print = kFALSE)
{
  // compute distribution of log likelihood value
  TH1D * hmc   = gStack[gPadNr][gOrder[gPadNr][0]];
  TH1D * hdata = gStack[gPadNr][gMaxProcess-1];
  Int_t nbins = hmc->GetNbinsX();
  Double_t loglik = loglikelihood(hmc, hdata, 1, nbins);
  TH1D * htest = new TH1D(*hdata);
  TH1D * lldistrib = new TH1D("lldistrib", "log(Likelihood) distribution", 
			      1000, loglik-200, loglik+200);
  setopt(lldistrib);
  for (Int_t n = 0; n < ntrials; n++) {
    // generate poisson around theorie
    for (Int_t i = 1; i <= nbins; i++) {
      htest->SetBinContent(i, gRandom->Poisson(hmc->GetBinContent(i)));
    }
    lldistrib->Fill(loglikelihood(hmc, htest, 1, nbins));
  }
  TCanvas * llcanvas = new TCanvas("llcanvas", "Log(Likelihood) distribution", 
				   40, 40, 800, 600);
  setopt(llcanvas);
  lldistrib->SetFillColor(kYellow);
  lldistrib->Draw();
  lldistrib->GetYaxis()->SetTitle("Anzahl Ereignisse");
  lldistrib->GetXaxis()->SetTitle("-ln L");
  // autozoom
  Int_t lowbin = 1;
  while (lldistrib->GetBinContent(lowbin) == 0)
    lowbin++;
  Int_t highbin = lldistrib->GetNbinsX();
  while (lldistrib->GetBinContent(highbin) == 0)
    highbin--;
  lldistrib->SetAxisRange(lldistrib->GetBinLowEdge(lowbin), 
			  lldistrib->GetBinLowEdge(highbin));
  TH1D * hworse = (TH1D *) lldistrib->Clone();
  for (Int_t nbin = 1; nbin < 501; nbin++) {
    hworse->SetBinContent(nbin, 0);
  }
  hworse->SetFillColor(95);
  hworse->Draw("same");
  Double_t pvalue = lldistrib->Integral(501,1000) / lldistrib->Integral();
  TLatex * tex = new TLatex(0.18, 0.96, Form("-ln L_{obs} = %5.2f", loglik));
  tex->SetNDC();
  tex->SetTextAlign(13);
  tex->Draw();
  tex = new TLatex(0.18, 0.86, Form("CL_{obs} = %.3f", pvalue));
  tex->SetNDC();
  tex->SetTextAlign(13);
  tex->Draw();
  TLine * l = new TLine(loglik, 0, loglik, lldistrib->GetMaximum());
  l->SetLineWidth(3);
  l->SetLineColor(kBlue);
  l->Draw();
  llcanvas->Modified();
  llcanvas->Update();
  if (print)
    llcanvas->Print("lldistrib.pdf");
  cd(gPadNr+1);
}

void Normalize(TH1D &h, double normalization, bool norm_per_avg_width) {
    int nbins = h.GetNbinsX();
    double low = h.GetBinLowEdge(1);
    double high = h.GetBinLowEdge(nbins+1);
    double width = (high-low)/nbins;
    if(norm_per_avg_width) normalization *= width;
    double integral = h.Integral("width");
    h.Scale(normalization/integral);
}

void AdjustDensityForBinWidth(TH1D &h) {
    double entries = h.GetEntries();
    int nbins = h.GetNbinsX();
    double low = h.GetBinLowEdge(1);
    double high = h.GetBinLowEdge(nbins+1);
    double width = (high-low)/nbins;
    for(int bin = 1; bin <= nbins; ++bin) {
        double content = h.GetBinContent(bin);
        double error = h.GetBinError(bin);
        double this_width = h.GetBinWidth(bin);
        double scale = width/this_width;
        h.SetBinContent(bin, content*scale);
        h.SetBinError(bin, error*scale);
    }
    h.SetEntries(entries);
}

// Rebin first histogram to match the second
TH1D *tools::Rebin(const TH1D *h, const TH1D* href) {

  //assert(href->GetNbinsX()<=h->GetNbinsX());
  if (!(href->GetNbinsX()<=h->GetNbinsX())) {
    cout << "Histo has less bins than ref: "
	 << h->GetNbinsX() << " vs " << href->GetNbinsX()
	 << " for " << h->GetName() << endl;
  }
  
  // First, we need to rebin inclusive jets to match b-tagged jets
  TH1D *hre = (TH1D*)href->Clone(Form("%s_rebin",h->GetName()));
  hre->Reset();

  for (int i = 1; i != h->GetNbinsX()+1; ++i) {

    double x = h->GetBinLowEdge(i);
    int j = hre->FindBin(x);
    // Check that h is fully contained within href bin
    if (h->GetBinContent(i)!=0) {

      if (!(h->GetBinLowEdge(i)>=hre->GetBinLowEdge(j) - 1e-5 &&
	    h->GetBinLowEdge(i+1)<=hre->GetBinLowEdge(j+1) + 1e-5)) {
	cerr << Form("Warning, bin edges overlapping: h=[%1.0f,%1.0f],"
		     " hre=[%1.0f,%1.0f] (%s)",
		     h->GetBinLowEdge(i), h->GetBinLowEdge(i+1),
		     hre->GetBinLowEdge(j), hre->GetBinLowEdge(j+1),
		     h->GetName()) << endl;
      }

      double y = ( hre->GetBinContent(j)*hre->GetBinWidth(j)
		   + h->GetBinContent(i)*h->GetBinWidth(i) )
	/ hre->GetBinWidth(j);
      //double ey = ( hre->GetBinError(j)*hre->GetBinWidth(j)
      //	    + h->GetBinError(i)*h->GetBinWidth(i) )
      // / hre->GetBinWidth(j);
      double ey = sqrt( pow(hre->GetBinError(j)*hre->GetBinWidth(j),2)
			+ pow(h->GetBinError(i)*h->GetBinWidth(i),2) )
	/ hre->GetBinWidth(j);
      hre->SetBinContent(j, y);
      hre->SetBinError(j, ey);
    }
  } // for i

  return hre;
} // Rebin

void FindMCToDataScaleFactor() {

  std::string data_filename = "~/data/out/v92/total.root";
  TFile* data_file = new TFile(data_filename.c_str(), "READ");
  std::string data_plotname = "TME_Al50_EvdE/all_particles/SiL_EvdE";
  TH2F* hEvdE_data = (TH2F*) data_file->Get(data_plotname.c_str());
  hEvdE_data->SetDirectory(0);
  data_file->Close();

  std::string MC_filename = "plots.root";
  TFile* MC_file = new TFile(MC_filename.c_str(), "READ");
  std::string MC_plotname = "hAll_EvdE_SiL";
  TH2F* hEvdE_MC = (TH2F*) MC_file->Get(MC_plotname.c_str());
  hEvdE_MC->SetDirectory(0);
  MC_file->Close();

  double energy_slice = 3000;
  int energy_slice_bin = hEvdE_data->GetXaxis()->FindBin(energy_slice);
  TH1D* hDataProjection = hEvdE_data->ProjectionY("_py", energy_slice_bin, energy_slice_bin);  
  energy_slice_bin = hEvdE_MC->GetXaxis()->FindBin(energy_slice);
  TH1D* hMCProjection = hEvdE_MC->ProjectionY("_py", energy_slice_bin, energy_slice_bin);

  TFitResultPtr data_fit_result = hDataProjection->Fit("gaus", "S");
  TFitResultPtr MC_fit_result = hMCProjection->Fit("gaus", "S");

  hDataProjection->SetLineColor(kBlack);
  hDataProjection->GetXaxis()->SetRangeUser(0, energy_slice);
  hDataProjection->Draw();
  hMCProjection->SetLineColor(kRed);
  hMCProjection->GetXaxis()->SetRangeUser(0, energy_slice);
  hMCProjection->Draw("SAME");

  //EvdE_data->Draw("COLZ");
  double data_mean = data_fit_result->Parameter(1);
  double MC_mean = MC_fit_result->Parameter(1);
  double scale_factor = data_mean / MC_mean;
  std::cout << "data / MC = " <<  data_mean << " / " << MC_mean << " = " << scale_factor << std::endl;

  // Go through the MC projection and scale each bin by the scale factor
  int n_bins = hMCProjection->GetNbinsX();
  double min_x = hMCProjection->GetXaxis()->GetXmin();
  double max_x = hMCProjection->GetXaxis()->GetXmax();
  TH1F* hMCProjection_scaled = new TH1F("hMCProjection_scaled", "", n_bins,min_x,max_x);

  for (int i_bin = 1; i_bin <= n_bins; ++i_bin) {
    double old_energy = hMCProjection->GetBinLowEdge(i_bin);
    double new_energy = old_energy * scale_factor;

    double old_bin_content = hMCProjection->GetBinContent(i_bin);
    hMCProjection_scaled->Fill(new_energy, old_bin_content);
  }

  hMCProjection_scaled->Draw();
  hDataProjection->Draw("SAME");
}

void reweighthisto(TH1D* hBPt, TH2D* h2D, TH1D *hRAA=0, bool weightByJpsiRAA=0, bool weightByBRAA=0)
{
   TH1D *htmp = h2D->ProjectionX("htmp");
   TH1D *htmp2 = h2D->ProjectionY("htmp2");
   
   for (int x=1;x<=h2D->GetNbinsX()+1;x++){
      if (htmp->GetBinContent(x)==0) continue;
      double ratio = hBPt->GetBinContent(x)/htmp->GetBinContent(x);
      if (weightByBRAA) ratio *= hRAA->GetBinContent(hRAA->FindBin(htmp->GetBinLowEdge(x)));   
       
      for (int y=1;y<=h2D->GetNbinsY()+1;y++){
         double ratio2 = ratio;
         if (weightByJpsiRAA) ratio2 *= hRAA->GetBinContent(hRAA->FindBin(htmp2->GetBinLowEdge(y)));   
         double val = h2D->GetBinContent(x,y)*ratio2;
         double valError = h2D->GetBinError(x,y)*ratio2;
	 h2D->SetBinContent(x,y,val);
	 h2D->SetBinError(x,y,valError);
      }   
   }
   delete htmp;
   
}

//.........这里部分代码省略.........
	if(fDoPhotonMT2Shape){
	prediction->ChPhotonsMT2 = new Channel("PhotonsMT2", "photon[0].lv.Pt()", cutStreamPhotonMT2.str().c_str(), fTriggerStream.str().c_str(),fSamplesPhotonPt);
	prediction->ChPhotonsMT2->fVerbose =prediction->fVerbose;
	prediction->ChPhotonsMT2->GetShapes("PhotonsMT2", "MET", 40, 0, 800);
	}

	// Znunu Hadronic Search MT2
	if(fDoZnunuMT2Shape){
	prediction->ChZnunuMT2 = new Channel("ZnunuMT2", "misc.MET", cutStreamZnunuMT2.str().c_str(), fTriggerStream.str().c_str(),fSamplesZnunu);
	prediction->ChZnunuMT2->fVerbose =prediction->fVerbose;
	prediction->ChZnunuMT2->GetShapes("ZnunuMT2", "MET", 40, 0, 800);
	}

	if(fDoZnunuMT2Shape && fDoPhotonMT2Shape){
	TH1D* ZnunuToPhotonMT2Ratio = prediction->GetRatio(prediction->ChZnunuMT2->hZJetsToNuNu, prediction->ChPhotonsMT2->hPhotons, 1);
	DrawHisto(ZnunuToPhotonMT2Ratio,ZnunuToPhotonMT2Ratio->GetName(),"EX0");
	}


	
	// Do Pt spectra comparison plot
	if(fDoPtSpectraComparison && fDoPhotonPtShape && fDoZllPtShape){
		*fLogStream<< "************************* produce pr spectra plot ****************************** " << endl;
		TH1D* hZllMC_cp      = prediction->GetScaledHisto(prediction->ChZllPt->hZJetsToLL, fDoPtSpectraComparisonScaling?(prediction->ChZllPt->hData->Integral())/(prediction->ChZllPt->hZJetsToLL->Integral())    :1, 0, 1);
		TH1D* hZllData_cp    = prediction->GetScaledHisto(prediction->ChZllPt->hData,    1, 0, 1) ;
		TH1D* hPhotonMC_cp   = prediction->GetScaledHisto(prediction->ChPhotonPt->hPhotons,fDoPtSpectraComparisonScaling?(prediction->ChPhotonPt->hData->Integral())/(prediction->ChPhotonPt->hPhotons->Integral()):1, 0, 1);
		TH1D* hPhotonData_cp = prediction->GetScaledHisto(prediction->ChPhotonPt->hData,    1, 0, 1);

		if(fDoPtSpectraComparisonAccCorr){
			TFile *f = new TFile("../RootMacros/ZllAcceptance.root", "OPEN");
			TH1D*  hZllAcc = (TH1D*) f->Get("ZJetsToLL_GenZllPtAcc_ZJetsToLL_GenZllPt_Ratio");
			if (hZllAcc==0) {cout << "WARNING: could not get histo ZJetsToLL_GenZllPtAcc_ZJetsToLL_GenZllPt_Ratio" << endl; exit(1);}
			for(int i=1; i<=hZllMC_cp->GetNbinsX(); ++i){
				if(hZllAcc->GetBinLowEdge(i) != hZllMC_cp->GetBinLowEdge(i)) {cout << "Zll Acc Correction: binnin does not match!!" << endl; exit(1);}
				if(hZllMC_cp  ->GetBinContent(i)<=0) continue;
				double acc_eff   = hZllAcc->GetBinContent(i);
				double orig_mc   = hZllMC_cp     ->GetBinContent(i);
				double orig_data = hZllData_cp   ->GetBinContent(i);
				cout << "bin i " << i << " acc eff " << acc_eff << " orig_mc " << orig_mc << " become " << orig_mc/acc_eff 
				     << " orig_data " << orig_data << " becomes " << orig_data/acc_eff << endl;
				hZllMC_cp   ->SetBinContent(i, orig_mc  /acc_eff);
				hZllData_cp ->SetBinContent(i, orig_data/acc_eff);
			}
			delete f;
		}
		
		hZllMC_cp->SetMarkerStyle(22);
		hZllMC_cp->SetLineColor(kOrange);
		hZllMC_cp->SetMarkerColor(kOrange);
		hZllMC_cp->SetMarkerSize(1.2);
		hZllData_cp->SetMarkerStyle(26);
		hZllData_cp->SetMarkerColor(kBlack);
		hZllData_cp->SetLineColor(kBlack);
		hPhotonMC_cp->SetLineColor(kMagenta+2);
		hPhotonMC_cp->SetMarkerColor(kMagenta+2);
		hPhotonMC_cp->SetMarkerStyle(20);
		hPhotonMC_cp->SetMarkerSize(1.2);
		hPhotonData_cp->SetMarkerStyle(4);
		hPhotonData_cp->SetMarkerColor(kBlack);
		hPhotonData_cp->SetLineColor(kBlack);

		TCanvas *col = new TCanvas("ZllToGammaPtSpectra", "", 0, 0, 500, 500);
		col->SetFillStyle(0);
		col->SetFrameFillStyle(0);
//		col->cd();
		gPad->SetFillStyle(0);

//.........这里部分代码省略.........
  leg->SetBorderSize(0);
  leg->SetTextSize(0.04);
  leg->SetTextFont(42);
  leg->SetFillStyle(0);
  leg->AddEntry(h,"Data","pl");
  leg->AddEntry(f,"Fit","l");
  leg->AddEntry(mass,"D^{0}+#bar{D^{#lower[0.2]{0}}} Signal","f");
  leg->AddEntry(massSwap,"K-#pi swapped","f");
  leg->AddEntry(background,"Combinatorial","l");
  leg->Draw("same");

  TLatex* texCms = new TLatex(0.18,0.93, "#scale[1.25]{CMS} Preliminary");
  texCms->SetNDC();
  texCms->SetTextAlign(12);
  texCms->SetTextSize(0.04);
  texCms->SetTextFont(42);
  texCms->Draw();

  TLatex* texCol = new TLatex(0.96,0.93, Form("%s #sqrt{s_{NN}} = 5.02 TeV",collisionsystem.Data()));
  texCol->SetNDC();
  texCol->SetTextAlign(32);
  texCol->SetTextSize(0.04);
  texCol->SetTextFont(42);
  texCol->Draw();

  TLatex* texPt = new TLatex(0.22,0.78,Form("%.1f < p_{T} < %.1f GeV/c",ptmin,ptmax));
  texPt->SetNDC();
  texPt->SetTextFont(42);
  texPt->SetTextSize(0.04);
  texPt->SetLineWidth(2);
  texPt->Draw();

  TLatex* texY = new TLatex(0.22,0.83,"|y| < 1.0");
  texY->SetNDC();
  texY->SetTextFont(42);
  texY->SetTextSize(0.04);
  texY->SetLineWidth(2);
  texY->Draw();

  TLatex* texYield = new TLatex(0.22,0.73,Form("N_{D} = %.0f #pm %.0f",yield,yieldErr));
  texYield->SetNDC();
  texYield->SetTextFont(42);
  texYield->SetTextSize(0.04);
  texYield->SetLineWidth(2);
  texYield->Draw();

  c->SaveAs(Form("plotFits/DMass_expo_%s_%.0f_%.0f.pdf",collisionsystem.Data(),ptmin,ptmax));
  
  TCanvas* cPull = new TCanvas(Form("cPull_%.0f_%.0f",ptmin,ptmax),"",600,700);
  TH1D* hPull = (TH1D*)h->Clone("hPull");
  for(int i=0;i<h->GetNbinsX();i++)
    {
      Double_t nfit = f->Integral(h->GetBinLowEdge(i+1),h->GetBinLowEdge(i+1)+h->GetBinWidth(i+1))/h->GetBinWidth(i+1);
      hPull->SetBinContent(i+1,(h->GetBinContent(i+1)-nfit)/h->GetBinError(i+1));
      hPull->SetBinError(i+1,0);
    }
  hPull->SetMinimum(-4.);
  hPull->SetMaximum(4.);
  hPull->SetYTitle("Pull");
  hPull->GetXaxis()->SetTitleOffset(1.);
  hPull->GetYaxis()->SetTitleOffset(0.65);
  hPull->GetXaxis()->SetLabelOffset(0.007);
  hPull->GetYaxis()->SetLabelOffset(0.007);
  hPull->GetXaxis()->SetTitleSize(0.12);
  hPull->GetYaxis()->SetTitleSize(0.12);
  hPull->GetXaxis()->SetLabelSize(0.1);
  hPull->GetYaxis()->SetLabelSize(0.1);
  hPull->GetYaxis()->SetNdivisions(504);
  TLine* lPull = new TLine(1.7, 0, 2., 0);
  lPull->SetLineWidth(1);
  lPull->SetLineStyle(7);
  lPull->SetLineColor(1);
  TPad* pFit = new TPad("pFit","",0,0.3,1,1);
  pFit->SetBottomMargin(0);
  pFit->Draw();
  pFit->cd();
  h->Draw("e");
  background->Draw("same");
  mass->Draw("same");
  massSwap->Draw("same");
  f->Draw("same");
  leg->Draw("same");
  texCms->Draw();
  texCol->Draw();
  texPt->Draw();
  texY->Draw();
  texYield->Draw();
  cPull->cd();
  TPad* pPull = new TPad("pPull","",0,0,1,0.3);
  pPull->SetTopMargin(0);
  pPull->SetBottomMargin(0.3);
  pPull->Draw();
  pPull->cd();
  hPull->Draw("p");
  lPull->Draw();
  cPull->cd();
  cPull->SaveAs(Form("plotFits/DMass_expo_%s_%.0f_%.0f_Pull.pdf",collisionsystem.Data(),ptmin,ptmax));

  return mass;
}

void plotAnaMult4(char *infname="ana.root")
{    
   TFile *inf = new TFile(infname);
   
   TTree *tData = (TTree*) inf->Get("Data_tree");
   TTree *tMC   = (TTree*) inf->Get("MC_tree");
   TTree *tPP   = (TTree*) inf->Get("PP_tree");

   TCanvas *c1 = new TCanvas("c","",1200,700);
//   c->Divide(4,1);
   TCut recoCut = "leadingJetPt>120&&subleadingJetPt>50&&dphi>5*3.14159265358979/6.&&abs(leadingJetEta)<1.6&&abs(subleadingJetEta)<1.6";
   TCut genCut = "genleadingJetPt>120&&gensubleadingJetPt>50&&genDphi>5*3.14159265358979/6.&&abs(genleadingJetEta)<1.6&&abs(gensubleadingJetEta)<1.6";

   
   const int nPtBin=5;
   Float_t PtBins[nPtBin+1] = {0.001,0.5,1,1.5,2,3.2};
   
   Int_t centBin[5] = {200,100,60,20,0};
   makeMultiPanelCanvas(c1,4,2,0.0,0.0,0.2,0.2,0.02);
   
   for (int i=0;i<4;i++) {
      c1->cd(i+1);
      TH1D * empty=new TH1D("empty","",100,0,3.19);
//      empty->Fill(0.5,1000); 
      empty->SetMaximum(39.99); 
      empty->SetMinimum(0.001); 
      empty->SetNdivisions(105,"X");
      empty->GetXaxis()->SetTitleSize(28);
      empty->GetXaxis()->SetTitleFont(43); 
      empty->GetXaxis()->SetTitleOffset(1.8);
      empty->GetXaxis()->SetLabelSize(22);
      empty->GetXaxis()->SetLabelFont(43);
      empty->GetYaxis()->SetTitleSize(28);
      empty->GetYaxis()->SetTitleFont(43); 
      empty->GetYaxis()->SetTitleOffset(1.8);
      empty->GetYaxis()->SetLabelSize(22);
      empty->GetYaxis()->SetLabelFont(43);
      empty->GetXaxis()->CenterTitle();
      empty->GetYaxis()->CenterTitle();

      empty->SetXTitle("#Delta#eta_{1,2}");
      empty->SetYTitle("Multiplicity Difference");
   
      TProfile *pData = new TProfile(Form("pData%d",i),"",nPtBin,PtBins);
      TProfile *pMC = new TProfile(Form("pMC%d",i),"",nPtBin,PtBins);
      TProfile *pPP = new TProfile(Form("pPP%d",i),"",nPtBin,PtBins);
      TProfile *pGen = new TProfile(Form("pGen%d",i),"",nPtBin,PtBins);
      TCut centCut = Form("hiBin>=%d&&hiBin<%d",centBin[i+1],centBin[i]);
      tData->Draw(Form("-multDiff:abs(leadingJetEta-subleadingJetEta)>>pData%d",i),recoCut&&centCut);
      tPP->Draw(Form("-multDiff:abs(leadingJetEta-subleadingJetEta)>>pPP%d",i),recoCut);
      tMC->Draw(Form("-multDiff:abs(leadingJetEta-subleadingJetEta)>>pMC%d",i),recoCut&&centCut);  
      tMC->Draw(Form("-genMultDiff:abs(genleadingJetEta-gensubleadingJetEta)>>pGen%d",i),genCut&&centCut);  

      pMC->SetLineColor(2);
      pMC->SetMarkerColor(2);
      pPP->SetLineColor(4);
      pPP->SetMarkerColor(4);
      
//      pData->SetAxisRange(0,50,"Y");
//      pData->SetAxisRange(0,0.49,"X");
      empty->Draw();

      double diff=0;
      if (i==0) diff=0.1;
      drawText(Form("%d-%d %%",(int)(0.5*centBin[i+1]),(int)(0.5*centBin[i])),0.22+diff,0.65);
      if (i==0) drawText("PbPb #sqrt{s_{NN}}=2.76 TeV 150/#mub",0.22+diff,0.85);
      if (i==0) drawText("pp      #sqrt{s_{NN}}=2.76 TeV 5.3/pb",0.22+diff,0.75);
      if (i==3) drawText("CMS Preliminary",0.3+diff,0.85);
            
      Float_t sys[4]={1,1,2.5,3};
      

      for (int j=1;j<=pData->GetNbinsX();j++)
      {
         TBox *b = new TBox(pData->GetBinLowEdge(j),pData->GetBinContent(j)-sys[i],pData->GetBinLowEdge(j+1),pData->GetBinContent(j)+sys[i]);
	 //b->SetFillColor(kGray);
	 b->SetFillStyle(0);
	 b->SetLineColor(1);
	 b->Draw();
         TBox *b2 = new TBox(pPP->GetBinLowEdge(j),pPP->GetBinContent(j)-1,pPP->GetBinLowEdge(j+1),pPP->GetBinContent(j)+1);
	 //b2->SetFillColor(65);
	 b2->SetFillStyle(0);
	 b2->SetLineColor(4);
	 b2->Draw();
      }
      pData->Draw("same");
      pMC->Draw("same");
      pPP->Draw("same");
      pGen->SetMarkerColor(4);
      pGen->SetMarkerStyle(24);
      pData->Draw("same");
      
//      pGen->Draw(" same");

      c1->cd(5+i);
      TH1D *empty2 = (TH1D*)empty->Clone("empty2");
      empty2->SetYTitle("PbPb - pp");
      empty2->SetMinimum(-5);
      empty2->SetMaximum(+29.99);
      empty2->Draw();
//.........这里部分代码省略.........

void compare_emu(std::string mcfile1, std::string mcfile2, 
		std::string var1, std::string var2="",
		float xmin=-9999.0, float xmax=-9999.0,
		std::string output="test",
		std::string header="Before unfolding",
		std::string mcName1="e",
		std::string mcName2="#mu",
		bool logScale=false
		)
{
  
  setTDRStyle();
  gStyle->SetOptStat(0);

  TH1F* h1;
  TH1F* h2;

  char tempName[300];
  if(var2 ==  "" )var2=var1;
  
  // first get the histogram files
  TFile *fmc1 = TFile::Open(mcfile1.data());
  TFile *fmc2   = TFile::Open(mcfile2.data());

  h1  = (TH1F*)(fmc1->Get(var1.data()));
  h2    = (TH1F*)(fmc2->Get(var2.data()));

  TH1D* hscale =(TH1D*) h1->Clone("hscale");
  hscale->SetYTitle(Form("%s/%s",mcName1.data(),mcName2.data()));

  h1->GetXaxis()->SetNdivisions(5);
  h1->GetYaxis()->SetDecimals();

  h2->GetXaxis()->SetNdivisions(5);
  h2->GetYaxis()->SetDecimals();

  hscale->GetXaxis()->SetNdivisions(5);
  hscale->GetYaxis()->SetDecimals();


  h1->SetLineColor(2);
  h1->SetMarkerColor(2);
  h1->SetMarkerSize(1);
  h1->SetMarkerStyle(24);


  h2->SetLineColor(4);
  h2->SetMarkerColor(4);
  h2->SetMarkerSize(1);
  h2->SetMarkerStyle(21);

  // if normalizing to the same area, set the scale 

  int binLo = -1;
  int binHi = -1;
  int nbins = h1->GetNbinsX();
  if(xmin>-9999.0 && xmax>-9999.0)
    {

      binLo = h1->FindBin(xmin);
      binHi = h1->FindBin(xmax)-1;

    }

  else
    {
      binLo = 1;
      binHi = nbins;
      xmin = h1->GetBinLowEdge(1);
      xmax = h1->GetBinLowEdge(nbins+1);
    }



  h1->Sumw2();
  h1->Scale(1.0/h1->Integral(binLo, binHi));

  h2->Sumw2();
  h2->Scale(1.0/h2->Integral(binLo, binHi));

  cout << "h2 integral = " << h2->Integral() << endl;
  cout << "h1 integral = "   << h1->Integral() << endl;;

  // get the ratio
  hscale->Divide(h1, h2, 1,1);
  
  cout << "===========================================================" << endl;
  cout << "Central value " << endl;
  for(int i=1;i<=hscale->GetNbinsX();i++)
    cout << "Bin " << i << " ( " << hscale->GetBinLowEdge(i) << "~" 
	 << hscale->GetBinLowEdge(i+1) << " ): " 
	 << h1->GetBinContent(i) << "/" << h2->GetBinContent(i) << " = " 
	 << hscale->GetBinContent(i) << " +- " << hscale->GetBinError(i) 
	 << endl;

  cout << "Uncertainty " << endl;

  for(int i=1;i<=hscale->GetNbinsX();i++)
    cout << "Bin " << i << " ( " << hscale->GetBinLowEdge(i) << "~" 
	 << hscale->GetBinLowEdge(i+1) << " ): " 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	jettrk_.trkmul.push_back(trkcorr[2]);
	jettrk_.trksec.push_back(trkcorr[3]);
	// ------------------------
	// calculate particle jet correlations
	// ------------------------
	Double_t pdphi[2]={ 9999,9999 };
	Double_t pdr[2]={ 9999,9999 };
	Double_t pdrbg[2]={ 9999,9999 };
	for (Int_t j=0; j<2; ++j) {
	  // signal
	  pdphi[j] = reco::deltaPhi(p_[i].phi(),anaJets_[j].phi());
	  pdr[j] = reco::deltaR(p_[i].eta(),p_[i].phi(),anaJets_[j].eta(),anaJets_[j].phi());
	  // bcksub
	  // * If don't do event selection, make eta reflection the default bkg axis
	  if (cut.BkgSubType.Contains("EtaRefl"||!doEvtSel_)) {
	    pdrbg[j] = reco::deltaR(p_[i].eta(),p_[i].phi(),-1*anaJets_[j].eta(),anaJets_[j].phi());
	  }
	  // monitor histograms
	  hPJDPhi[j]->Fill(pdphi[j],trackWeight);
	  jettrk_.pdr[j].push_back(pdr[j]);
	  jettrk_.pdrbg[j].push_back(pdrbg[j]);
	}

	Float_t trkEnergy=p_[i].pt();
	// ------------------------
	// met calculation
	// ------------------------
	Float_t pptx=cos(pdphi[0])*trkEnergy*trackWeight;
	metx+=pptx;
	Float_t ppty=sin(pdphi[0])*trkEnergy*trackWeight;
	mety+=ppty;
	if (fabs(pdr[0])<0.8||fabs(pdr[1])<0.8) metConex+=pptx; else metOutOfConex+=pptx;
	//for (int i=0;i<hPt->GetNbinsX()+2;++i) cout << "Bin " << i << " ledge: " << hPt->GetBinLowEdge(i) << endl;
	if (trkEnergy>=hPt->GetBinLowEdge(1)&&trkEnergy<hPt->GetBinLowEdge(2)) {
	  metx0+=pptx;
	  mety0+=ppty;
	  if (fabs(pdr[0])<0.8||fabs(pdr[1])<0.8) metConex0+=pptx; else metOutOfConex0+=pptx;
	} else if (trkEnergy>=hPt->GetBinLowEdge(2)&&trkEnergy<hPt->GetBinLowEdge(3)) {
	  metx1+=pptx;
	  mety1+=ppty;
	  if (fabs(pdr[0])<0.8||fabs(pdr[1])<0.8) metConex1+=pptx; else metOutOfConex1+=pptx;
	} else if (trkEnergy>=hPt->GetBinLowEdge(3)&&trkEnergy<hPt->GetBinLowEdge(4)) {
	  metx2+=pptx;
	  mety2+=ppty;
	  if (fabs(pdr[0])<0.8||fabs(pdr[1])<0.8) metConex2+=pptx; else metOutOfConex2+=pptx;
	} else if (trkEnergy>=hPt->GetBinLowEdge(4)&&trkEnergy<hPt->GetBinLowEdge(5)) {
	  metx3+=pptx;
	  mety3+=ppty;
	  if (fabs(pdr[0])<0.8||fabs(pdr[1])<0.8) metConex3+=pptx; else metOutOfConex3+=pptx;
	} else if (trkEnergy>=hPt->GetBinLowEdge(5)&&trkEnergy<hPt->GetBinLowEdge(6)) { 
	  metx4+=pptx;
	  mety4+=ppty;
	  if (fabs(pdr[0])<0.8||fabs(pdr[1])<0.8) metConex4+=pptx; else metOutOfConex4+=pptx;
	} else if (trkEnergy>=hPt->GetBinLowEdge(6)&&trkEnergy<hPt->GetBinLowEdge(7)) { 
	  metx5+=pptx;
	  mety5+=ppty;
	  if (fabs(pdr[0])<0.8||fabs(pdr[1])<0.8) metConex5+=pptx; else metOutOfConex5+=pptx;
	}

	// =====================================================
	// Calculate Cone Sums
	// =====================================================
	// cone sum
	for (Int_t j=0; j<2; ++j) {
	  for (Int_t b=0; b<numPtBins; ++b) {
	    if (trkEnergy>=hPt->GetBinLowEdge(b+1)&&trkEnergy<hPt->GetBinLowEdge(b+2)) {

//.........这里部分代码省略.........
      }
      else {
	cout<< "Using inclusive template "<<inclusiveTemplateName<<" for region: HT "<<ht_LOW<<"-"<<ht_HI<<" and NJ "<<njets_LOW<<"-"<<njets_HI<<endl; 
	// Get the template
	lastbin_hybrid = makeHybridTemplate(srname, h_MT2Template, inclusiveTemplateName , fData, fZinv, fDY, lastmt2val_hybrid);
	cout<<"lastbin_hybrid "<<lastbin_hybrid<<" and lastmt2val_hybrid "<<lastmt2val_hybrid<<endl;
	if (h_MT2Template!=0) h_MT2Template->Print();
	else cout<<"h_MT2Template is 0, using hybrid within this region (no external templates)"<<endl;
      }

   


    }

  


    // If there is no template for this region, go back to the standard hybrid
    if (doHybridSimple || (doHybridInclusiveTemplate && h_MT2Template == 0)) {
      // hybrid method: use nominal MC CR yield histogram to determine how many MT2 bins to use
      //  by default: use all MT2 bins integrated (no bin-by-bin).
      //  choose the last bin to try to have at least hybrid_nevent_threshold integrated events

      // Calculate last bin on local histogram
      for ( int ibin=1; ibin <= hDY->GetNbinsX(); ++ibin ) {
	float top = 0, integratedYield = 0;
	//if (hDataEM) top = hDataEM->Integral(ibin,-1)*rSFOF;
	integratedYield = hDY->Integral(ibin,-1) - top;
	if (integratedYield < hybrid_nevent_threshold) {
	  if (ibin == 1) lastbin_hybrid = 1;
	  else {
	    lastbin_hybrid = ibin-1;
	    lastmt2val_hybrid = hDY->GetBinLowEdge(ibin);
	  }
	  break;
	}
      }
      cout<<"lastbin_hybrid for doHybridSimple: "<<lastbin_hybrid<<endl;
    }
    
    

    TH1D* ratio = (TH1D*) hZinv->Clone("ratio");
    ratio->Divide(hDY);
    double errNum, errDen;
    float ratioValue = hZinv->IntegralAndError(1,-1,errNum) / hDY->IntegralAndError(1,-1,errDen);
    float ratioErr = ratioValue*sqrt(pow(errNum/hZinv->Integral(), 2) + pow(errDen/hDY->Integral(),2));

    TH1D* CRyield = (TH1D*) hData->Clone("h_mt2binsCRyield");

    TH1D* purityMC = (TH1D*) hDY->Clone("h_mt2binsPurityMC");
    if (hTop) purityMC->Add(hTop, -1);
    purityMC->Divide(hDY);

    TH1D* purityData = (TH1D*) hData->Clone("h_mt2binsPurityData");
    if (hDataEM) purityData->Add(hDataEM, -1*rSFOF); 
    purityData->Divide(purityData, hData, 1, 1, "B");
    
    TH1D* Stat = (TH1D*) CRyield->Clone("h_mt2binsStat");
    Stat->Multiply(purityData);
    Stat->Multiply(ratio);

    TH1D* Syst = (TH1D*) Stat->Clone("h_mt2binsSyst");
    TH1D* pred = (TH1D*) Stat->Clone("h_mt2bins");
    for ( int ibin = 0; ibin <= Stat->GetNbinsX(); ++ibin) { 

//_______________________________________
void FitMCshape(AliDielectronSignalBase *sig)
{
  TFile mcFile(mcLineShapeFile);
  if (!mcFile.IsOpen()) {
    printf("mcMinv_LHC10e2 not found!!!\n");
    return;
  }
  TH1D *hMmc = (TH1D*)mcFile.Get("mcMinv");
  if (!hMmc) {
    printf("mcMinv not found!!\n");
    return;
  }
  hMmc->SetDirectory(0);
  hMmc->SetName("mcMinv");
  mcFile.Close();
  
  TH1* hMsub=sig->GetSignalHistogram();
  Double_t mb1=sig->GetIntegralMin();
  Double_t mb2=sig->GetIntegralMax();
  
  Double_t effInt = 0.;
  for(Int_t iBin=hMmc->FindBin(mb1); iBin<hMmc->FindBin(mb2); iBin++) {
    effInt += hMmc->GetBinContent(iBin);
  }
  effInt/=hMmc->Integral();
  printf("MC signal fraction in range %4.2f-%4.2f GeV: %5.3f \n",hMmc->GetBinLowEdge(hMmc->FindBin(mb1)),hMmc->GetBinLowEdge(hMmc->FindBin(mb2)+1),effInt);
 
  Float_t mcScale1=(hMsub->GetXaxis()->GetBinWidth(1)/hMmc->GetXaxis()->GetBinWidth(1))*
    hMsub->Integral(hMsub->FindBin(mb1),hMsub->FindBin(mb2))/
    hMmc->Integral(hMmc->FindBin(mb1),hMmc->FindBin(mb2));
  
  printf("1st guess of MC scale factor: %6.3f\n",mcScale1);
  
  Float_t mcScale=0.;
  Float_t chi2_min=100000.;
  Int_t iMin=0;
  Int_t ndf=0;
  
  for(Int_t i=0; i<20; i++){
    Float_t chi2=0.;
    Float_t scale=(0.4+0.05*(Float_t)i)*mcScale1;
    ndf=0;
    for(Int_t ib=1; ib<=hMsub->GetXaxis()->GetNbins(); ib++){
      Float_t data=(Float_t)hMsub->GetBinContent(ib);
      Float_t err=(Float_t)hMsub->GetBinError(ib);
      Float_t mc=scale*((Float_t)hMmc->GetBinContent(hMmc->FindBin(hMsub->GetBinCenter(ib))));
      if (err>0) {
        chi2 += ((data-mc)*(data-mc))/(err*err);
        ndf++;
      } else {
        //printf("bin %d Err: %6.3f, chi2: %6.1f\n",ib,err,chi2);
      }
    }
      //printf("%d scale factor: %6.3f, chi2: %6.1f\n",i,scale,chi2);
    if(chi2 < chi2_min){
      chi2_min = chi2;
      mcScale = scale;
      iMin=i;
    }
  }
  //Float_t chi2dof=chi2_min/(Float_t)(hMinv->GetXaxis()->GetNbins()-1);
  Float_t chi2dof=chi2_min/((Float_t)(ndf-1));
  printf("MC fit (i=%d): chi2/dof: %6.3f/%d, Scale: %7.4f \n",iMin,chi2_min,(ndf-1),mcScale);
  hMmc->SetTitle(Form("%f",chi2dof));
  
  //mcScale=IntData/IntMC;printf("Int Data, MC: %10.1f %10.1f, MC scale: %6.3f\n",IntData,IntMC,mcScale);
  
  hMmc->Scale(mcScale);
  hMmc->SetOption("sameHISTC");
  hMsub->GetListOfFunctions()->Add(hMmc);
}

void Smooth (TString sel) 
{

  const int nbins = 50;

  const int nvars=17;
  TString var[nvars] = { "C8", "M8", "C6", "M6", "MEt", "MEtSig", "CorrSumEt", "GoodHt", 
		      "M45bestall", "Chi2mass", "Chi2extall", "Mbbnoh", "DPbbnoh", 
		      "SumHED4", "SumHED6", "DP12", "MEtDP2" };
  TString pippo[nvars];
  TString pippotot[nvars];
  TString pippotth[nvars];
  TString pippototS[nvars];
  TString pippotthS[nvars];
  for ( int i=0; i<nvars; i++ ) { 
    pippo[i] = var[i]+sel; 
    pippotot[i]=var[i]+sel+"_bgr";
    pippotth[i]=var[i]+sel+"_sig";
    pippototS[i]=var[i]+sel+"_bgrS";
    pippotthS[i]=var[i]+sel+"_sigS";
  }
  
  const int nqcdsamples=8;
  TFile * QCD[nqcdsamples];     
  QCD[0] = new TFile("./root/et30_eta2.5/TDAna_QCD30-50_tk3.root");
  QCD[1] = new TFile("./root/et30_eta2.5/TDAna_QCD50-80_tk3.root");
  QCD[2] = new TFile("./root/et30_eta2.5/TDAna_QCD80-120_tk3.root");
  QCD[3] = new TFile("./root/et30_eta2.5/TDAna_QCD120-170_tk3.root");
  QCD[4] = new TFile("./root/et30_eta2.5/TDAna_QCD170-230_tk3.root");
  QCD[5] = new TFile("./root/et30_eta2.5/TDAna_QCD230-300_tk3.root");
  QCD[6] = new TFile("./root/et30_eta2.5/TDAna_QCD300-380_tk3.root");
  QCD[7] = new TFile("./root/et30_eta2.5/TDAna_QCD380incl_tk3.root");
  double QCDxs[nqcdsamples] = { 155929000., 20938850., 2949713., 499656., 100995.,  23855., 6391., 2821.};
  double NQCD[nqcdsamples] = { 86000., 78000., 104000., 96000., 100000., 102000., 112000., 102000.};


  const int nwsamples=11;
  TFile * W[nwsamples];
  W[0] = new TFile ("./root/et30_eta2.5/TDAna_W0w_tk3.root");
  W[1] = new TFile ("./root/et30_eta2.5/TDAna_W10w_tk3.root");
  W[2] = new TFile ("./root/et30_eta2.5/TDAna_W11w_tk3.root");
  W[3] = new TFile ("./root/et30_eta2.5/TDAna_W20w_tk3.root");
  W[4] = new TFile ("./root/et30_eta2.5/TDAna_W21w_tk3.root");
  W[5] = new TFile ("./root/et30_eta2.5/TDAna_W30w_tk3.root");
  W[6] = new TFile ("./root/et30_eta2.5/TDAna_W31w_tk3.root");
  W[7] = new TFile ("./root/et30_eta2.5/TDAna_W40w_tk3.root");
  W[8] = new TFile ("./root/et30_eta2.5/TDAna_W41w_tk3.root");
  W[9] = new TFile ("./root/et30_eta2.5/TDAna_W50w_tk3.root");
  W[10] = new TFile ("./root/et30_eta2.5/TDAna_W51w_tk3.root");
  double Wxs[nwsamples] = { 45000., 9200., 250., 2500., 225., 590., 100., 125., 40., 85., 40. };
  double NW[nwsamples] = { 88000., 40000., 100530., 99523., 105255., 79000., 
			   88258., 83038., 30796., 59022., 41865. };

  TFile * TTH = new TFile("./root/et30_eta2.5/TDAna_ttH_120_tk3.root");
  double TTHxs = 0.667 ;
  double NTTH = 62000.; // 1652000.; // 62000.;

  const int nttsamples=5;
  TFile * TT[nttsamples];
  TT[0] = new TFile("./root/et30_eta2.5/TDAna_TT0_tk3.root");
  TT[1] = new TFile("./root/et30_eta2.5/TDAna_TT1_tk3.root");
  TT[2] = new TFile("./root/et30_eta2.5/TDAna_TT2_tk3.root");
  TT[3] = new TFile("./root/et30_eta2.5/TDAna_TT3_tk3.root");
  TT[4] = new TFile("./root/et30_eta2.5/TDAna_TT4_tk3.root");
  // double TTxs[5] = { 619., 176., 34.,  6., 1.5 };  // from web
  double TTxs[nttsamples] = { 434., 162., 43., 10., 1.9 };     // from note
  double NTT[nttsamples] = { 57900., 66000., 98159., 14768., 5304. };

  double Lumfactor = 100000; // 100/fb of luminosity assumed in histograms

  TH1D * Histo_TOT[nvars];
  TH1D * Histo_TTH[nvars];
  TH1D * Histo_TOTS[nvars];
  TH1D * Histo_TTHS[nvars];
  for ( int i=0; i<nvars; i++ ) {
    cout << i << endl;
    TH1D * H = dynamic_cast<TH1D*>(TTH->Get(pippo[i]));
    double minx=H->GetBinLowEdge(1);
    double maxx=nbins*H->GetBinWidth(1);
    Histo_TOT[i] = new TH1D ( pippotot[i]," ", nbins, minx, maxx );
    Histo_TTH[i] = new TH1D ( pippotth[i]," ", nbins, minx, maxx );
    Histo_TOTS[i] = new TH1D ( pippototS[i]," ", nbins, minx, maxx );
    Histo_TTHS[i] = new TH1D ( pippotthS[i]," ", nbins, minx, maxx );
  }

  cout << "Starting loop on variables needing smoothing" << endl;

  // Loop on variables 
  // -----------------
  for ( int ivar=0; ivar<nvars; ivar++ ) {

    // Extract sum histograms with the right normalization and errors
    // --------------------------------------------------------------
    double totWW[nwsamples][nbins]={0.};
    double totW[nbins]={0.};
    double s2_totW[nbins]={0.};
    double totNW[nwsamples][nbins]={0.};
    for ( int i=0; i<nwsamples; i++ ) {
      cout << "Processing W file #" << i << " ..." << endl;
      TH1D * Histo = dynamic_cast<TH1D*>(W[i]->Get(pippo[ivar]));
//.........这里部分代码省略.........