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

//.........这里部分代码省略.........
	}
	double spread = seperation/2;
	double spreadsq2 = 2*spread*spread;
	//Swap kaon and pion numbers

	for (int i = 1; i < hpion->GetNbinsX()/2; i++)
	{
		double t_swap = hpion->GetBinContent(i);
		hpion->SetBinContent(i,hpion->GetBinContent(hpion->GetNbinsX() - i + 1));
		hpion->SetBinContent(hpion->GetNbinsX() - i + 1, t_swap);
	}
	for (int i = 1; i < hkaon->GetNbinsX()/2; i++)
	{
		double t_swap = hkaon->GetBinContent(i);
		hkaon->SetBinContent(i,hkaon->GetBinContent(hkaon->GetNbinsX() - i + 1));
		hkaon->SetBinContent(hkaon->GetNbinsX() - i + 1, t_swap);
	}

	double titlesize=1.2*.04;

	int rebin = 20;
	rebin = 160;
	hpion->Rebin(rebin);
	hkaon->Rebin(rebin);

	hpion->SetAxisRange(hmin,hmax);
	hkaon->SetAxisRange(hmin,hmax);

	hkaon->GetXaxis()->SetTitle("Loglikelihood difference");
	hkaon->GetYaxis()->SetTitle("A.U.");
	hkaon->GetXaxis()->SetTitleSize(titlesize);
	hkaon->GetYaxis()->SetTitleSize(titlesize);

	hpion->SetStats(false);
	hkaon->SetStats(false);

	hpion->SetLineColor(kCyan);
	//hpion->SetFillColorAlpha(kRed,.5);

	hkaon->SetLineColor(kBlue);
	//hkaon->SetFillColorAlpha(kBlue,.5);

	TLegend *leg_ll = new TLegend(.6,.6,.8,.8);
	leg_ll->AddEntry(hpion,"Pion");
	leg_ll->AddEntry(hkaon,"Kaon");
	leg_ll->SetBorderSize(0);


	hkaon->SetTitle("log(P(Pi)/P(K)) for actual Pi (red) and K (blue) at 5 GeV");



	TH1F *pion_veto_eff = new TH1F(*hpion);
	TH1F *kaon_missid = new TH1F(*hkaon);

	pion_veto_eff->SetName("pion_veto_eff");
	pion_veto_eff->SetTitle("");

	kaon_missid->SetName("kaon_missid");
	kaon_missid->SetTitle("");



	for (int i = 0; i < pion_veto_eff->GetNbinsX(); i++)
	{
		pion_veto_eff->SetBinContent(i,hpion->Integral(0,i));

void MakeSystPlot(const TString& channel, TFile * input, RooWorkspace * ws, const RooArgList * obs, 
                  const Int_t p, const Int_t up, const Int_t down)
{
    TString pdfname = g_pdfname;
    pdfname.ReplaceAll("$CHANNEL", channel);

    const TString& process = g_processes[p];
    TString systUp = g_systematics[up];
    systUp.ReplaceAll("$CHANNEL", channel);
    systUp.ReplaceAll("$PROCESS", process);
    
    TString systDown = g_systematics[down];
    systDown.ReplaceAll("$CHANNEL", channel);
    systDown.ReplaceAll("$PROCESS", process);

    if (process != "Wj0b" && process != "Wj1b" && process != "Wj2b") {
        if (systUp.Contains("WJModel") || systDown.Contains("WJModel") || systUp.Contains("WJSlope") || systDown.Contains("WJSlope"))
            return;
    }

    if (process != "Zj0b" && process != "Zj1b" && process != "Zj2b") {
        if (systUp.Contains("ZJModel") || systDown.Contains("ZJModel") || systUp.Contains("ZJSlope") || systDown.Contains("ZJSlope"))
            return;
    }
    
    if (process != "TT") {
        if (systUp.Contains("TTModel") || systDown.Contains("TTModel"))
            return;
    }
    
    TH1F * h = (TH1F *) input->Get(channel + "/" + process);
    TH1F * hUp = (TH1F *) input->Get(channel + "/" + process + "_" + systUp);
    TH1F * hDown = (TH1F *) input->Get(channel + "/" + process + "_" + systDown);
    if ((h->Integral() > 0. && hUp->Integral() <= 0.) || h->Integral() <= 0.) {
        TString clonename = hUp->GetName();
        delete hUp;
        hUp = (TH1F*) h->Clone(clonename);
    }
    if ((h->Integral() > 0. && hDown->Integral() <= 0.) || h->Integral() <= 0.) {
        TString clonename = hDown->GetName();
        delete hDown;
        hDown = (TH1F*) h->Clone(clonename);
    }
    if (process == "WH" || process == "ZH" || process == "VH") {
        if (systUp.Contains("eff_b"))
            systUp.ReplaceAll("eff_b", "eff_b_SIG");
        if (systDown.Contains("eff_b"))
            systDown.ReplaceAll("eff_b", "eff_b_SIG");
    }
    RooDataHist * dhUp = new RooDataHist(process + "_" + systUp, "", *obs, hUp);
    ws->import(*dhUp);
    RooDataHist * dhDown = new RooDataHist(process + "_" + systDown, "", *obs, hDown);
    ws->import(*dhDown);
    
    h->SetStats(0);
    h->SetTitle("; BDT");
    h->SetLineColor(1);
    h->SetLineWidth(2);
    h->SetFillColor(0);
    h->SetMarkerStyle(20);
    h->SetMinimum(0.01);
    h->GetXaxis()->CenterTitle();
    
    hUp->SetLineColor(g_upcol);
    hUp->SetLineWidth(2);
    hUp->SetFillColor(0);
    
    hDown->SetLineColor(g_downcol);
    hDown->SetLineWidth(2);
    hDown->SetFillColor(0);
    
    h->Draw("e1");
    hUp->Draw("hist same");
    hDown->Draw("hist same");
    h->Draw("e1 same");
    
    TLegend * leg = new TLegend(0.35, 0.20, 0.92, 0.35);
    leg->SetFillColor(0);
    leg->SetFillStyle(0);
    leg->SetLineColor(0);
    leg->SetShadowColor(0);
    leg->SetTextFont(62);
    //leg->SetTextSize(0.015);
    leg->AddEntry(h, process, "pl");
    leg->AddEntry(hUp, systUp, "l");
    leg->AddEntry(hDown, systDown, "l");
    leg->Draw();

    gPad->RedrawAxis();
    gPad->Modified();
    gPad->Update();
    gPad->Print(pdfname);
    
    delete dhUp;
    delete dhDown;
    delete leg;
    
    return;
}

void Demo_TryExtrapolationInXT_Exp0_LogLogFits(Bool_t xt=kTRUE, Float_t expo=0.)
{
        SetStyle();
        gStyle->SetOptFile(0);
        gStyle->SetOptStat(0);
        gStyle->SetOptFit(0);
 
        xt=kTRUE;

	//define dummy histogram and some style parameters
	TH1F *dum;
 	dum = new TH1F("dum","",160,5e-4,0.3); 
	dum->SetMinimum(1e-14);
	dum->SetMaximum(1);
	dum->SetTitle(Form(";x_{T};#sqrt{s}^{%0.1f}  E d^{3}#sigma/dp^{3}",expo));
	dum->SetLineWidth(0);
	dum->SetStats(0);
	dum->GetXaxis()->CenterTitle();
	dum->GetYaxis()->CenterTitle();
	dum->GetXaxis()->SetTitleSize(0.05);
	dum->GetYaxis()->SetTitleSize(0.05);
	dum->GetXaxis()->SetTitleOffset(1.17);
	dum->GetYaxis()->SetTitleOffset(1.3);

	gROOT->LoadMacro("/net/hidsk0001/d00/scratch/krajczar/ppRefForpPb_PilotRun/interpolation_HIN10005_kk/data_table_to_graph.C");

	//get 7 TeV points
	TGraphErrors *cms_7000_g = data_table_to_graph("cms",7000,xt,expo);
	cms_7000_g->SetMarkerColor(kBlack);
	
	TF1 *cms_7000_fit = new TF1("cms_7000_fit","[0]*pow(1.0+(x/[1]),[2])",10.*2./7000.,0.1);//Fit from 10 GeV/c
	cms_7000_fit->SetLineWidth(1);
	cms_7000_fit->SetParameters(3e22,2.5e-4,-7);
	cms_7000_g->Fit(cms_7000_fit,"REMW0");
	
	//get 2.36 TeV points
//	TGraphErrors *cms_2360_g = data_table_to_graph("cms",2360,xt,expo);
//	cms_2360_g->SetMarkerColor(kMagenta+3);
	
//	TF1 *cms_2360_fit = new TF1("cms_2360_fit","[0]*pow(1.0+(x/[1]),[2])",2e-3,0.1);
//	cms_2360_fit->SetLineWidth(1);
//	cms_2360_fit->SetParameters(3e22,2.5e-4,-7);
//	cms_2360_g->Fit(cms_2360_fit,"REMW0");

	//get 2.76 TeV points (KK using existing txt files)
	TGraphErrors *cms_2760_g = data_table_to_graph("cms",2760,xt,expo);
	cms_2760_g->SetMarkerColor(kMagenta+3);
	TF1 *cms_2760_fit = new TF1("cms_2760_fit","[0]*pow(1.0+(x/[1]),[2])",10.*2./2760.,0.1);//Fit from 10 GeV/c
	cms_2760_fit->SetLineColor(kMagenta+3);
	cms_2760_fit->SetLineWidth(1);
	cms_2760_fit->SetParameters(3e22,2.5e-4,-7);
	cms_2760_g->Fit(cms_2760_fit,"REMW0");
	
	//get 1.96 TeV points
	TGraphErrors *cdf_1960_g = data_table_to_graph("cdf",1960,xt,expo);
	cdf_1960_g->SetMarkerColor(kOrange-3);
	cdf_1960_g->SetMarkerStyle(30);
	
	//TGraphErrors *cdfold_1960_g = data_table_to_graph("cdfold",1960,xt);
	//cdfold_1960_g->SetMarkerColor(kBlue);
	//cdfold_1960_g->SetMarkerStyle(30);
	
	TF1 *cdf_1960_fit = new TF1("cdf_1960_fit","[0]*pow(1.0+(x/[1]),[2])",2.*10./1960.,0.1);//Fit from 10 GeV/c
	cdf_1960_fit->SetLineColor(kOrange-3);
	cdf_1960_fit->SetLineWidth(1);
	cdf_1960_fit->SetParameters(3e22,2.5e-4,-7);
	cdf_1960_g->Fit(cdf_1960_fit,"REMW0");
	
	//get 1.8 TeV points
	TGraphErrors *cdf_1800_g = data_table_to_graph("cdf",1800,xt,expo);
	cdf_1800_g->SetMarkerColor(kGreen+3);
	cdf_1800_g->SetMarkerStyle(28);

	TF1 *cdf_1800_fit = new TF1("cdf_1800_fit","[0]*pow(1.0+(x/[1]),[2])",2.*10./1800.,0.1);//Fit from 10 GeV/c
	cdf_1800_fit->SetLineColor(kGreen+3);
	cdf_1800_fit->SetLineWidth(1);
	cdf_1800_fit->SetParameters(3e22,2.5e-4,-7.2);
	cdf_1800_fit->FixParameter(2,-7.2);
	cdf_1800_g->Fit(cdf_1800_fit,"REMW0");
	
	//get 0.9 TeV points
	TGraphErrors *cms_900_g = data_table_to_graph("cms",900,xt,expo);
	cms_900_g->SetMarkerColor(kRed);
	
	TF1 *cms_900_fit = new TF1("cms_900_fit","[0]*pow(1.0+(x/[1]),[2])",2.*10./900.,0.01);//Fit from 10 GeV/c
	cms_900_fit->SetLineColor(kRed);
	cms_900_fit->SetLineWidth(1);
	cms_900_fit->SetParameters(3e22,2.5e-4,-7);
	cms_900_g->Fit(cms_900_fit,"REMW0");
	
//	TGraphErrors *ua1_900_g = data_table_to_graph("ua1",900,xt,expo);
//	ua1_900_g->SetMarkerColor(kCyan+1);
//	ua1_900_g->SetMarkerStyle(26);

	// get 0.63 TeV points
	TGraphErrors *cdf_630_g = data_table_to_graph("cdf",630,xt,expo);
	cdf_630_g->SetMarkerColor(kOrange+3);
	cdf_630_g->SetMarkerStyle(27);

	//draw graphs to canvas
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  
  RooPlot* lFrame6 = lToy->plotNLL(0,1000,100);
  lFrame6->SetTitle("-log(L)");
  lFrame6->SetXTitle("-log(L)");
  TCanvas* lC05 = new TCanvas("logl","logl",600,600) ;
  lFrame6->Draw() ;
  lC05->SaveAs(("logL_toyfits_"+iBkg+"_" + iDir + "_" + iEnergy+".png").c_str());

  RooPlot* lFrame7 = lToy->plotParam(lNSig);
  lFrame7->SetTitle("distribution of values of N_{sig} after toy fit");
  lFrame7->SetXTitle("N_{sig}");
  TCanvas* lC06 = new TCanvas("Nsig","Nsig",600,600) ;
  lFrame7->Draw() ;
  lC06->SaveAs(("NSig_toyfits_"+iBkg+"_" + iDir + "_" + iEnergy+".png").c_str());
  
  RooPlot* lFrame8 = lToy->plotParam(lNB);
  lFrame8->SetTitle("distribution of values of N_{bkg} after toy fit");
  lFrame8->SetXTitle("N_{bkg}");
  TCanvas* lC07 = new TCanvas("Nbkg","Nbkg",600,600) ;
  lFrame8->Draw() ;
  lC07->SaveAs(("Nbkg_toyfits_"+iBkg+"_" + iDir + "_" + iEnergy+".png").c_str());
 

  if(iRebin) { 
    const int lNBins = lData->GetNbinsX();
    double *lAxis    = getAxis(lData);
    lH0     = rebin(lH0    ,lNBins,lAxis);
    lH      = rebin(lH     ,lNBins,lAxis);
    lHUp    = rebin(lHUp   ,lNBins,lAxis);
    lHDown  = rebin(lHDown ,lNBins,lAxis);
    lHUp1   = rebin(lHUp1  ,lNBins,lAxis);
    lHDown1 = rebin(lHDown1,lNBins,lAxis);
  }

  // we dont need this bin errors since we do not use them (fit tails replaces bin-by-bin error!), therefore i set all errors to 0, this also saves us from modifying the add_bbb_error.py script in which I otherwise would have to include a option for adding bbb only in specific ranges
  int lMergeBin = lH->GetXaxis()->FindBin(iFirst);
  for(int i0 = lMergeBin; i0 < lH->GetNbinsX()+1; i0++){
    lH->SetBinError  (i0,0);
    lHUp->SetBinError  (i0,0);
    lHDown->SetBinError  (i0,0);
    lHUp1->SetBinError  (i0,0);
    lHDown1->SetBinError  (i0,0);
  }


  TFile *lOutFile =new TFile("Output.root","RECREATE");
  cloneFile(lOutFile,lFile,iDir+"/"+iBkg);
  lOutFile->cd(iDir.c_str());
  lH     ->Write();
  lHUp   ->Write(); 
  lHDown ->Write(); 
  lHUp1  ->Write(); 
  lHDown1->Write(); 

  // Debug Plots
  lH0->SetStats(0);
  lH->SetStats(0);
  lHUp->SetStats(0);
  lHDown->SetStats(0);
  lHUp1->SetStats(0);
  lHDown1->SetStats(0);
  lH0    ->SetLineWidth(1); lH0->SetMarkerStyle(kFullCircle);
  lH     ->SetLineColor(kGreen);
  lHUp   ->SetLineColor(kRed);
  lHDown ->SetLineColor(kRed);
  lHUp1  ->SetLineColor(kBlue);
  lHDown1->SetLineColor(kBlue);
  TCanvas *lC0 = new TCanvas("Can","Can",800,600);
  lC0->Divide(1,2); lC0->cd();  lC0->cd(1)->SetPad(0,0.2,1.0,1.0); gPad->SetLeftMargin(0.2) ; 
  lH0->Draw();
  lH     ->Draw("hist sames");
  lHUp   ->Draw("hist sames");
  lHDown ->Draw("hist sames");
  lHUp1  ->Draw("hist sames");
  lHDown1->Draw("hist sames");
  gPad->SetLogy();
  
  TLegend* leg1;
  /// setup the CMS Preliminary
  leg1 = new TLegend(0.7, 0.80, 1, 1); 
  leg1->SetBorderSize( 0 );
  leg1->SetFillStyle ( 1001 );
  leg1->SetFillColor (kWhite);
  leg1->AddEntry( lH0 , "orignal",  "PL" );
  leg1->AddEntry( lH , "cental fit",  "L" );
  leg1->AddEntry( lHUp , "shift1 up",  "L" );
  leg1->AddEntry( lHDown , "shift1 down",  "L" );
  leg1->AddEntry( lHUp1 , "shift2 up",  "L" );
  leg1->AddEntry( lHDown1 , "shift2 down",  "L" );
  leg1->Draw("same");


  lC0->cd(2)->SetPad(0,0,1.0,0.2); gPad->SetLeftMargin(0.2) ;
  drawDifference(lH0,lH,lHUp,lHDown,lHUp1,lHDown1);
  lH0->SetStats(0);
  lC0->Update();
  lC0->SaveAs((iBkg+"_"+"CMS_"+iName+"1_" + iDir + "_" + iEnergy+".png").c_str());
  //lFile->Close();
  return;
}

void addNuisance(std::string iFileName,std::string iChannel,std::string iBkg,std::string iEnergy,std::string iName,std::string iDir,bool iRebin=true,bool iVarBin=false,int iFitModel=1,double iFirst=150,double iLast=1500) { 
  std::cout << "======> " << iDir << "/" << iBkg << " -- " << iFileName << std::endl;  
  if(iVarBin) addVarBinNuisance(iFileName,iChannel,iBkg,iEnergy,iName,iDir,iRebin,iFitModel,iFirst,iLast);
  if(iVarBin) return;

  TFile *lFile = new TFile(iFileName.c_str());
  TH1F  *lH0   = (TH1F*) lFile->Get((iDir+"/"+iBkg).c_str());
  TH1F  *lData = (TH1F*) lFile->Get((iDir+"/data_obs").c_str());

  //Define the fit function
  RooRealVar lM("m","m" ,0,5000);   //lM.setBinning(lBinning);
  RooRealVar lA("a","a" ,50,  0.1,100);
  RooRealVar lB("b","b" ,0.0 , -10.5,10.5); //lB.setConstant(kTRUE);
  RooDataHist *pH0  =  new RooDataHist("Data","Data" ,RooArgList(lM),lH0);
  RooGenericPdf *lFit  = 0; lFit = new RooGenericPdf("genPdf","exp(-m/(a+b*m))",RooArgList(lM,lA,lB));
  if(iFitModel == 1) lFit = new RooGenericPdf("genPdf","exp(-a*pow(m,b))",RooArgList(lM,lA,lB));
  if(iFitModel == 1) {lA.setVal(0.3); lB.setVal(0.5);}
  if(iFitModel == 2) lFit = new RooGenericPdf("genPdf","a*exp(b*m)",RooArgList(lM,lA,lB));
  if(iFitModel == 3) lFit = new RooGenericPdf("genPdf","a/pow(m,b)",RooArgList(lM,lA,lB));
  RooFitResult  *lRFit = 0;
  double lFirst = iFirst;
  double lLast  = iLast;
  //lRFit = lFit->chi2FitTo(*pH0,RooFit::Save(kTRUE),RooFit::Range(lFirst,lLast));
  lRFit = lFit->fitTo(*pH0,RooFit::Save(kTRUE),RooFit::Range(lFirst,lLast),RooFit::Strategy(0)); 
  TMatrixDSym lCovMatrix   = lRFit->covarianceMatrix(); 
  TMatrixD  lEigVecs(2,2);    lEigVecs = TMatrixDSymEigen(lCovMatrix).GetEigenVectors();
  TVectorD  lEigVals(2);      lEigVals = TMatrixDSymEigen(lCovMatrix).GetEigenValues();
  cout << " Ve---> " << lEigVecs(0,0) << " -- " << lEigVecs(1,0) << " -- " << lEigVecs(0,1) << " -- " << lEigVecs(1,1) << endl;
  cout << " Co---> " << lCovMatrix(0,0) << " -- " << lCovMatrix(1,0) << " -- " << lCovMatrix(0,1) << " -- " << lCovMatrix(1,1) << endl;
  double lACentral = lA.getVal();
  double lBCentral = lB.getVal();
  lEigVals(0) = sqrt(lEigVals(0));
  lEigVals(1) = sqrt(lEigVals(1));
  cout << "===> " << lEigVals(0) << " -- " << lEigVals(1) << endl;
  
  TH1F* lH     = (TH1F*) lFit->createHistogram("fit" ,lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));
  lA.setVal(lACentral + lEigVals(0)*lEigVecs(0,0));
  lB.setVal(lBCentral + lEigVals(0)*lEigVecs(1,0));
  TH1F* lHUp   = (TH1F*) lFit->createHistogram("Up"  ,lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));
  lA.setVal(lACentral - lEigVals(0)*lEigVecs(0,0));
  lB.setVal(lBCentral - lEigVals(0)*lEigVecs(1,0));
  TH1F* lHDown = (TH1F*) lFit->createHistogram("Down",lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));

  lA.setVal(lACentral + lEigVals(1)*lEigVecs(0,1));
  lB.setVal(lBCentral + lEigVals(1)*lEigVecs(1,1));
  TH1F* lHUp1   = (TH1F*) lFit->createHistogram("Up1",lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));
  lA.setVal(lACentral - lEigVals(1)*lEigVecs(0,1));
  lB.setVal(lBCentral - lEigVals(1)*lEigVecs(1,1));
  TH1F* lHDown1 = (TH1F*) lFit->createHistogram("Down1",lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));

  std::string lNuisance1 =  iBkg+"_"+"CMS_"+iName+"1_" + iChannel + "_" + iEnergy;
  std::string lNuisance2 =  iBkg+"_"+"CMS_"+iName+"2_" + iChannel + "_" + iEnergy;
  lHUp    = merge(lNuisance1 + "Up"   ,lFirst,lH0,lHUp);
  lHDown  = merge(lNuisance1 + "Down" ,lFirst,lH0,lHDown);
  lHUp1   = merge(lNuisance2 + "Up"   ,lFirst,lH0,lHUp1);
  lHDown1 = merge(lNuisance2 + "Down" ,lFirst,lH0,lHDown1);
  lH      = merge(lH0->GetName()      ,lFirst,lH0,lH);

  if(iRebin) { 
    const int lNBins = lData->GetNbinsX();
    double *lAxis    = getAxis(lData);
    lH0     = rebin(lH0    ,lNBins,lAxis);
    lH      = rebin(lH     ,lNBins,lAxis);
    lHUp    = rebin(lHUp   ,lNBins,lAxis);
    lHDown  = rebin(lHDown ,lNBins,lAxis);
    lHUp1   = rebin(lHUp1  ,lNBins,lAxis);
    lHDown1 = rebin(lHDown1,lNBins,lAxis);
  }

  // we dont need this bin errors since we do not use them (fit tails replaces bin-by-bin error!), therefore i set all errors to 0, this also saves us from modifying the add_bbb_error.py script in which I otherwise would have to include a option for adding bbb only in specific ranges
  int lMergeBin = lH->GetXaxis()->FindBin(iFirst);
  for(int i0 = lMergeBin; i0 < lH->GetNbinsX()+1; i0++){
    lH->SetBinError  (i0,0);
    lHUp->SetBinError  (i0,0);
    lHDown->SetBinError  (i0,0);
    lHUp1->SetBinError  (i0,0);
    lHDown1->SetBinError  (i0,0);
  }


  TFile *lOutFile =new TFile("Output.root","RECREATE");
  cloneFile(lOutFile,lFile,iDir+"/"+iBkg);
  lOutFile->cd(iDir.c_str());
  lH     ->Write();
  lHUp   ->Write(); 
  lHDown ->Write(); 
  lHUp1  ->Write(); 
  lHDown1->Write(); 

  // Debug Plots
  lH0->SetStats(0);
  lH->SetStats(0);
  lHUp->SetStats(0);
  lHDown->SetStats(0);
  lHUp1->SetStats(0);
  lHDown1->SetStats(0);
  lH0    ->SetLineWidth(1); lH0->SetMarkerStyle(kFullCircle);
  lH     ->SetLineColor(kGreen);
  lHUp   ->SetLineColor(kRed);
  lHDown ->SetLineColor(kRed);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	HedpNoCut->Add(HedpNoCutR,HedpNoCutL);
	HedpNoCut->SetTitle(titledp);
	HedpNoCut->Draw();
	HedpCut->Add(HedpCutR,HedpCutL);
	HedpCut->SetLineColor(kBlack);
	HedpCut->SetFillColor(kViolet);
	HedpCut->Draw("same");

	pad2->cd();
	cout << "Drawing Target Cut..." << endl;
	TString titleTarget = "Target Cut";
	TH1F *HeReactZNoCut = new TH1F("HeReactZNoCut",titleTarget,400,-0.3,0.3);
	TH1F *HeReactZCut = new TH1F("HeReactZCut",titleTarget,400,-0.3,0.3);
	chainHeR->Draw("ReactPt_R.z>>HeReactZNoCutR(400,-0.3,0.3)", "", "");
	chainHeR->Draw("ReactPt_R.z>>HeReactZCutR(400,-0.3,0.3)", targetR, "");
	chainHeL->Draw("ReactPt_L.z>>HeReactZNoCutL(400,-0.3,0.3)", "", "");
	chainHeL->Draw("ReactPt_L.z>>HeReactZCutL(400,-0.3,0.3)", targetL, "");
	HeReactZNoCut->Add(HeReactZNoCutR,HeReactZNoCutL);
	HeReactZNoCut->SetTitle(titleTarget);
	HeReactZNoCut->Draw();
	HeReactZCut->Add(HeReactZCutR,HeReactZCutL);
	HeReactZCut->SetLineColor(kBlack);
	HeReactZCut->SetFillColor(kViolet);
	HeReactZCut->Draw("same");

	pad3->cd();
	cout << "Drawing Theta and Phi..." << endl;
	TString titleThetaPhiNoCut = "Theta and Phi, No Cut";
	TH2F *HeThetaPhiNoCut = new TH2F("HeThetaPhiNoCut",titleThetaPhiNoCut,100,-0.05,0.05,100,-0.1,0.1);
	chainHeR->Draw("ExTgtCor_R.th:ExTgtCor_R.ph>>HeThetaPhiNoCutR(100,-0.05,0.05,100,-0.1,0.1)", targetR, "");
	chainHeL->Draw("ExTgtCor_L.th:ExTgtCor_L.ph>>HeThetaPhiNoCutL(100,-0.05,0.05,100,-0.1,0.1)", targetL, "");
	HeThetaPhiNoCut->Add(HeThetaPhiNoCutR,HeThetaPhiNoCutL);
	HeThetaPhiNoCut->SetTitle(titleThetaPhiNoCut);
	HeThetaPhiNoCut->SetStats(kFALSE);
	HeThetaPhiNoCut->Draw("COLZ");

	pad4->cd();
	cout << "Drawing Theta and Phi Cut..." << endl;
	TString titleThetaPhi = "Theta and Phi Cut";
	TH2F *HeThetaPhiCut = new TH2F("HeThetaPhiCut",titleThetaPhi,100,-0.05,0.05,100,-0.1,0.1);
	chainHeR->Draw("ExTgtCor_R.th:ExTgtCor_R.ph>>HeThetaPhiCutR(100,-0.05,0.05,100,-0.1,0.1)", targetR && thetaphiR, "");
	chainHeL->Draw("ExTgtCor_L.th:ExTgtCor_L.ph>>HeThetaPhiCutL(100,-0.05,0.05,100,-0.1,0.1)", targetL && thetaphiL, "");
	HeThetaPhiCut->Add(HeThetaPhiCutR,HeThetaPhiCutL);
	HeThetaPhiCut->SetTitle(titleThetaPhi);
	HeThetaPhiCut->SetStats(kFALSE);
	HeThetaPhiCut->Draw("COLZ");

	pad5->cd();
	TString titleQ2 = "Q2";
	TH1F *histQ2 = new TH1F("histQ2",titleQ2,400,0,1.6);
	cout << "Drawing " << titleQ2 << "..." << endl;
	chainHeR->Draw("PriKineR.Q2>>histQ2R(400,0,1.6)", cutsR, "");	
	chainHeL->Draw("PriKineL.Q2>>histQ2L(400,0,1.6)", cutsL, "");
	histQ2->Add(histQ2R,histQ2L);
	histQ2->SetTitle(titleQ2);
	histQ2->SetStats(kFALSE);
	histQ2->Draw();

	pad6->cd();
	TString titleNu = "Nu";
	cout << "Drawing " << titleNu << "..." << endl;
	TH1F *histNu = new TH1F("histNu",titleNu,100,0.3,0.8);
	chainHeR->Draw("PriKineR.nu>>histNuR(100,0.3,0.8)", cutsR, "");
	chainHeL->Draw("PriKineL.nu>>histNuL(100,0.3,0.8)", cutsL, "");
	histNu->Add(histNuR,histNuL);
	histNu->SetTitle(titleNu);

void ootpu_comparison(TString files, TString var, TString title, int nbins, float low, float high, TString comments="") {

  TChain* chain = new TChain("reduced_tree");
  chain->Add(files);

  TH1::SetDefaultSumw2();

  TH1F* hA = new TH1F("hA",title, nbins, low, high);
  TH1F* hB = new TH1F("hB",title, nbins, low, high);
  TH1F* hC = new TH1F("hC",title, nbins, low, high);
  TH1F* hD = new TH1F("hD",title, nbins, low, high);
  hA->SetStats(0);
  hA->GetYaxis()->SetLabelSize(0.04);
  //hA->GetYaxis()->SetTitleOffset(1.3);
  hA->GetXaxis()->SetTitleOffset(1.1);
  hA->GetXaxis()->SetTitleFont(132);
  hA->GetXaxis()->SetTitleSize(0.042);
  hA->GetXaxis()->SetLabelSize(0.04);
  hA->SetLineWidth(2);

  //  hA->StatOverflows(true);
  //  hB->StatOverflows(true);
  //  hC->StatOverflows(true);
  //  hD->StatOverflows(true);

  int n1(0), n2(0), n3(0), n4(0), n5(0);
  if (files.Contains("20bx25")) {
    n1=5;
    n2=17;
    n3=21;
    n4=25;
    n5=40;
  }
  else if (files.Contains("S14")) {
    n1=0;
    n2=25;
    n3=40;
    n4=55;
    n5=120;
  }
  else { // default: 8 TeV scenario
    n1=0;
    n2=15;
    n3=22;
    n4=32;
    n5=70;
  }


  TString mu("num_gen_muons==1&&muon_reco_match>=0");
  //if (files.Contains("PU_S10")) {
  TString cuts = Form("(%s)&&(eoot_pu>=%d&&eoot_pu<%d)",mu.Data(),n1,n2);
  cout << "Cuts: " << cuts.Data() << endl;
  chain->Project("hA", var, cuts);
  cuts = Form("(%s)&&(eoot_pu>=%d&&eoot_pu<%d)",mu.Data(),n2,n3);
  cout << "Cuts: " << cuts.Data() << endl;
  chain->Project("hB", var, cuts);
  cuts = Form("(%s)&&(eoot_pu>=%d&&eoot_pu<%d)",mu.Data(),n3,n4);
  cout << "Cuts: " << cuts.Data() << endl;
  chain->Project("hC", var, cuts);
  cuts = Form("(%s)&&(eoot_pu>=%d)",mu.Data(),n4);
  cout << "Cuts: " << cuts.Data() << endl;
  chain->Project("hD", var, cuts);
  // }
  // else {
  //  }
  
  float avg1(hA->GetMean());
  float avg2(hB->GetMean());
  float avg3(hC->GetMean());
  float avg4(hD->GetMean());

  hA->Scale(1/hA->Integral());
  hB->Scale(1/hB->Integral());
  hC->Scale(1/hC->Integral());
  hD->Scale(1/hD->Integral());

  hA->SetLineColor(1);
  hB->SetLineColor(2);
  hC->SetLineColor(3);
  hD->SetLineColor(4);

  hA->SetLineWidth(2);
  hB->SetLineWidth(2);
  hC->SetLineWidth(2);
  hD->SetLineWidth(2);

  TCanvas* c1 = new TCanvas();
  float max = TMath::Max(hA->GetMaximum(), hB->GetMaximum());
  if (hC->GetMaximum()>max) max = hC->GetMaximum();
  if (hD->GetMaximum()>max) max = hD->GetMaximum();

  hA->SetMaximum(max*1.1);
  hA->Draw("e1");
  hB->Draw("e1,same");
  hC->Draw("e1,same");
  hD->Draw("e1,same");

  TLegend* leg = new TLegend(0.42,0.6,0.9,0.9);
  leg->SetFillStyle(0);
//.........这里部分代码省略.........