C++ TCanvas::Divide方法代码示例

void All(){
	string varname = "topness";
	//string signal = "T2tt_850_100";
	string signal = "T2tt_650_100";
	//string region = "baseline";
	//string region = "baselineSmallDM";
	string region = "baselineSmallDM300";
	//string region = "baselineSmallDM2b";
	string channel = "allChannels";
	/*
	Optimization("topness", signal, region, channel);
	Optimization("topness_m1", signal, region, channel);
	Optimization("topness_m2", signal, region, channel);
	Optimization("topness_m3", signal, region, channel);
	Optimization("topness_m4", signal, region, channel);
	Optimization("topness_m5", signal, region, channel);
	Optimization("topness_m6", signal, region, channel);
	Optimization("topness_m7", signal, region, channel);
	Optimization("topness_m8", signal, region, channel);
	Optimization("topness_m9", signal, region, channel);
	Optimization("topness_m10", signal, region, channel);
	Optimization("topness_m11", signal, region, channel);
	Optimization("topness_m12", signal, region, channel);
	Optimization("topness_m13", signal, region, channel);
*/
	OptRes res_650_ref = Optimization(varname, "T2tt_650_100", region, channel,"plots/config_def.txt/1DStack.root");
	OptRes res_850_ref = Optimization(varname, "T2tt_850_100", region, channel,"plots/config_def.txt/1DStack.root");
	Shape(varname, region, channel,"plots/config_def.txt/1DSuperimposed.root");
	//Optimization("topness", signal, region, channel,"plots/config_def.txt/pow2/1DStack.root");
	//Optimization("topness", signal, region, channel,"plots/config_def.txt/pow2/b1nonb/1DStack.root");
	/*
	OptRes res_650_comb3 = Optimization(varname, "T2tt_650_100", region, channel,"plots/config_def.txt/comb3_all/1DStack.root");
	OptRes res_850_comb3 = Optimization(varname, "T2tt_850_100", region, channel,"plots/config_def.txt/comb3_all/1DStack.root");
	OptRes res_650_comb4 = Optimization(varname, "T2tt_650_100", region, channel,"plots/config_def.txt/comb4_all/1DStack.root");
	OptRes res_850_comb4 = Optimization(varname, "T2tt_850_100", region, channel,"plots/config_def.txt/comb4_all/1DStack.root");
	OptRes res_650_comb5 = Optimization(varname, "T2tt_650_100", region, channel,"plots/config_def.txt/comb5_all/1DStack.root");
	OptRes res_850_comb5 = Optimization(varname, "T2tt_850_100", region, channel,"plots/config_def.txt/comb5_all/1DStack.root");
	OptRes res_650_comb4_b1 = Optimization(varname, "T2tt_650_100", region, channel,"plots/config_def.txt/comb4/1DStack.root");
	OptRes res_850_comb4_b1 = Optimization(varname, "T2tt_850_100", region, channel,"plots/config_def.txt/comb4/1DStack.root");
	OptRes res_650_comb3_4 = Optimization(varname, "T2tt_650_100", region, channel,"plots/config_def.txt/comb3_4/1DStack.root");
	OptRes res_850_comb3_4 = Optimization(varname, "T2tt_850_100", region, channel,"plots/config_def.txt/comb3_4/1DStack.root");
		
		TCanvas* cROC = new TCanvas();
		TLegend* leg = new TLegend(0.6,0.7,0.8,0.9);
		cROC->Divide(2);
		cROC->cd(1);
		res_650_ref.gROC->SetTitle("ROC - T2tt_650-100");
		res_650_ref.gROC->Draw("ACP");
		res_650_comb3.gROC->SetLineColor(kGreen);
		res_650_comb3.gROC->Draw("same");
		res_650_comb4.gROC->SetLineColor(kRed);
		res_650_comb4.gROC->Draw("same");
		res_650_comb5.gROC->SetLineColor(kBlue);
		res_650_comb5.gROC->Draw("same");
		res_650_comb4_b1.gROC->SetLineColor(kOrange);
		res_650_comb4_b1.gROC->Draw("same");
		res_650_comb3_4.gROC->SetLineColor(kViolet);
		res_650_comb3_4.gROC->Draw("same");
		leg->AddEntry(res_650_ref.gROC,"default","l");
		leg->AddEntry(res_650_comb3.gROC,"3 jets","l");
		leg->AddEntry(res_650_comb4.gROC,"4 jets","l");
		leg->AddEntry(res_650_comb5.gROC,"5 jets","l");
		leg->AddEntry(res_650_comb4_b1.gROC,"4 jets (nb<=1)","l");
		leg->AddEntry(res_650_comb3_4.gROC,"4 jets (nb<=1) && 3 jets (nb>=2)","l");
		leg->Draw("same");
		cROC->cd(2);
		res_850_ref.gROC->SetTitle("ROC - T2tt_850-100");
		res_850_ref.gROC->Draw("ACP");
		res_850_comb3.gROC->SetLineColor(kGreen);
		res_850_comb3.gROC->Draw("same");
		res_850_comb4.gROC->SetLineColor(kRed);
		res_850_comb4.gROC->Draw("same");
		res_850_comb5.gROC->SetLineColor(kBlue);
		res_850_comb5.gROC->Draw("same");
		res_850_comb4_b1.gROC->SetLineColor(kOrange);
		res_850_comb4_b1.gROC->Draw("same");
		res_850_comb3_4.gROC->SetLineColor(kViolet);
		res_850_comb3_4.gROC->Draw("same");
		leg->Draw("same");
		cROC->Draw();
		cROC->Print("combjets.png");
 	return 0;	
*/	
	//for(int i=12;i<13;i++){
	//for(int i=45;i<60;i++){
	for(int i=20;i<22;i++){
	//for(int i=1;i<45;i++){
	//for(int i=1;i<9;i++){
		char dir[120]; 
		sprintf(dir,"plots/config2_%d.txt/",i);
		//sprintf(dir,"plots/config2_%d.txt/comb3_4/",i);
		//sprintf(dir,"plots/config_def.txt/comb3_4/");

		
		//sprintf(dir,"plots/tp2_config_%d.txt/pow2/1b1nonb/",i);
		//sprintf(dir,"plots/tp2_config_%d.txt/pow2/loopnonb_all/",i);
		//sprintf(dir,"plots/tp2_config_%d.txt/pow2/1DStack.root",i);
		//sprintf(dir,"plots/tp2_config_%d.txt/pow2/loopnonb/1DStack.root",i);
		//sprintf(dir,"plots/tp2_config_%d.txt/pow2/loopnonb_all/1DStack.root",i);
		//sprintf(dir,"plots/tp2_config_%d.txt/pow2/comb4/1DStack.root",i);
//.........这里部分代码省略.........

void plot_CL_chi2_roofit(char * filename, double min, double max, double initial, double ndof_min, double ndof_max, char * plot, char * var = "chi2")
{

	//gStyle->SetOptStat(0);
	//gStyle->SetOptFit(1);
	//gStyle->SetStatFontSize(0.02);
	TFile * _file0 = TFile::Open(filename);
	TTree * t = (TTree*)_file0->Get("tuple");

	RooRealVar * chi2 = new RooRealVar(var, "#chi^{2}", min, max);
	RooRealVar * ndof = new RooRealVar("ndof", "ndof", initial, ndof_min, ndof_max);
	RooChiSquarePdf * pdf = new RooChiSquarePdf("pdf", "pdf", *chi2, *ndof);

	RooDataSet * data = new RooDataSet("data", "data", RooArgSet(*chi2), RooFit::Import(*t));

	pdf->fitTo(*data);

	char formula[30];
	sprintf(formula, "TMath::Prob(%s,ndof)", var);

	RooFormulaVar * CL_ndof_eff_formula = new RooFormulaVar("CL","CL(#chi^{2})",formula, RooArgList(*chi2, *ndof));
	RooRealVar * CL_ndof_eff = (RooRealVar*) data->addColumn(*CL_ndof_eff_formula);
	CL_ndof_eff->setRange(0, 1);
	RooUniform * uniform = new RooUniform("uniform", "uniform", *CL_ndof_eff);
	uniform->fitTo(*data);

	//RooFormulaVar * CL_ndof_min_formula = new RooFormulaVar("CL","CL(#chi^{2})","TMath::Prob(chi2,39)", RooArgList(*chi2));
	//RooRealVar * CL_ndof_min = (RooRealVar*) data->addColumn(*CL_ndof_min_formula);
	//CL_ndof_min->setRange(0, 1);

	RooPlot * frame0 = chi2->frame(RooFit::Bins(25));
	data->plotOn(frame0);
	pdf->plotOn(frame0);
	pdf->paramOn(frame0, RooFit::Format("NELU", RooFit::AutoPrecision(2)), RooFit::Layout(0.6,0.95,0.75)); 
	data->statOn(frame0, RooFit::Format("NELU", RooFit::AutoPrecision(2)), RooFit::Layout(0.6,0.95,0.95));

	RooPlot * frame1 = CL_ndof_eff->frame(RooFit::Bins(10));
	data->plotOn(frame1);
	uniform->plotOn(frame1);	

        TCanvas * c = new TCanvas("c","c",1200, 600);
        c->Divide(2,1);
        c->cd(1);
	frame0->Draw();
        c->cd(2);
	frame1->Draw();

/*
        char buf[30];
        sprintf(buf, "TMath::Prob(chi2,%f)>>h1", f1->GetParameter(0));
        cout << buf << endl;
        c->Modified();
        c->Update();

        c->cd(2);
        t->Draw("TMath::Prob(chi2,ndof-8)>>h0");
        t->Draw(buf);
        h1->Draw();
        h1->Fit("pol0");
        h0->Draw("same");
        h1->GetXaxis()->SetTitle("CL(#chi^{2})");
        h1->GetYaxis()->SetTitle("Number of toys / 0.1");
        h1->SetMinimum(0);
        h1->SetMaximum(2*t->GetEntries()/nbins);
*/
        c->SaveAs(plot);



}

//.........这里部分代码省略.........
  // Visualize 1-sigma error encoded in fit result 'r' as orange band using linear error propagation
  // This results in an error band that is by construction symmetric
  //
  // The linear error is calculated as
  // error(x) = Z* F_a(x) * Corr(a,a') F_a'(x)
  //
  // where     F_a(x) = [ f(x,a+da) - f(x,a-da) ] / 2, 
  // 
  //         with f(x) = the plotted curve 
  //              'da' = error taken from the fit result
  //        Corr(a,a') = the correlation matrix from the fit result
  //                Z = requested significance 'Z sigma band'
  //
  // The linear method is fast (required 2*N evaluations of the curve, where N is the number of parameters), 
  // but may not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and 
  // Gaussian approximations made
  //
  model.plotOn(frame,VisualizeError(*r,1),FillColor(kOrange)) ;


  // Calculate error using sampling method and visualize as dashed red line. 
  //
  // In this method a number of curves is calculated with variations of the parameter values, as sampled 
  // from a multi-variate Gaussian p.d.f. that is constructed from the fit results covariance matrix. 
  // The error(x) is determined by calculating a central interval that capture N% of the variations 
  // for each valye of x, where N% is controlled by Z (i.e. Z=1 gives N=68%). The number of sampling curves 
  // is chosen to be such that at least 100 curves are expected to be outside the N% interval, and is minimally 
  // 100 (e.g. Z=1->Ncurve=356, Z=2->Ncurve=2156)) Intervals from the sampling method can be asymmetric, 
  // and may perform better in the presence of strong correlations, but may take (much) longer to calculate
  model.plotOn(frame,VisualizeError(*r,1,kFALSE),DrawOption("L"),LineWidth(2),LineColor(kRed)) ;
  
  // Perform the same type of error visualization on the background component only.
  // The VisualizeError() option can generally applied to _any_ kind of plot (components, asymmetries, efficiencies etc..)
  model.plotOn(frame,VisualizeError(*r,1),FillColor(kOrange),Components("bkg")) ;
  model.plotOn(frame,VisualizeError(*r,1,kFALSE),DrawOption("L"),LineWidth(2),LineColor(kRed),Components("bkg"),LineStyle(kDashed)) ;

  // Overlay central value
  model.plotOn(frame) ;
  model.plotOn(frame,Components("bkg"),LineStyle(kDashed)) ;
  d->plotOn(frame) ;
  frame->SetMinimum(0) ;


  // V i s u a l i z e   p a r t i a l   f i t   e r r o r 
  // ------------------------------------------------------

  // Make plot frame
  RooPlot* frame2 = x.frame(Bins(40),Title("Visualization of 2-sigma partial error from (m,m2)")) ;
  
  // Visualize partial error. For partial error visualization the covariance matrix is first reduced as follows
  //        ___                   -1
  // Vred = V22  = V11 - V12 * V22   * V21
  //
  // Where V11,V12,V21,V22 represent a block decomposition of the covariance matrix into observables that
  // are propagated (labeled by index '1') and that are not propagated (labeled by index '2'), and V22bar
  // is the Shur complement of V22, calculated as shown above  
  //
  // (Note that Vred is _not_ a simple sub-matrix of V)

  // Propagate partial error due to shape parameters (m,m2) using linear and sampling method
  model.plotOn(frame2,VisualizeError(*r,RooArgSet(m,m2),2),FillColor(kCyan)) ;
  model.plotOn(frame2,Components("bkg"),VisualizeError(*r,RooArgSet(m,m2),2),FillColor(kCyan)) ;
  
  model.plotOn(frame2) ;
  model.plotOn(frame2,Components("bkg"),LineStyle(kDashed)) ;
  frame2->SetMinimum(0) ;
 

  // Make plot frame
  RooPlot* frame3 = x.frame(Bins(40),Title("Visualization of 2-sigma partial error from (s,s2)")) ;
  
  // Propagate partial error due to yield parameter using linear and sampling method
  model.plotOn(frame3,VisualizeError(*r,RooArgSet(s,s2),2),FillColor(kGreen)) ;
  model.plotOn(frame3,Components("bkg"),VisualizeError(*r,RooArgSet(s,s2),2),FillColor(kGreen)) ;
  
  model.plotOn(frame3) ;
  model.plotOn(frame3,Components("bkg"),LineStyle(kDashed)) ;
  frame3->SetMinimum(0) ;


  // Make plot frame
  RooPlot* frame4 = x.frame(Bins(40),Title("Visualization of 2-sigma partial error from fsig")) ;
  
  // Propagate partial error due to yield parameter using linear and sampling method
  model.plotOn(frame4,VisualizeError(*r,RooArgSet(fsig),2),FillColor(kMagenta)) ;
  model.plotOn(frame4,Components("bkg"),VisualizeError(*r,RooArgSet(fsig),2),FillColor(kMagenta)) ;
  
  model.plotOn(frame4) ;
  model.plotOn(frame4,Components("bkg"),LineStyle(kDashed)) ;
  frame4->SetMinimum(0) ;


  
  TCanvas* c = new TCanvas("rf610_visualerror","rf610_visualerror",800,800) ;
  c->Divide(2,2) ;
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4)  ; frame->Draw() ;
  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.6) ; frame2->Draw() ;
  c->cd(3) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.6) ; frame3->Draw() ;
  c->cd(4) ; gPad->SetLeftMargin(0.15) ; frame4->GetYaxis()->SetTitleOffset(1.6) ; frame4->Draw() ;
}

void CsIProj()
{
  TFile *file = new TFile("../root/NZ_55_New.root");
  TFile *gates = new TFile("../gates/zlines.root");
  TFile *gates2 = new TFile("../gates/zlines_new.root");

  ofstream ofile("CsI_55A_New.dat");
  
  TCanvas *mycan = (TCanvas*)gROOT->FindObjectAny("mycan");
  if(!mycan)
    {
      mycan = new TCanvas("mycan","mycan");
      mycan->Divide(1,2);
    }
  
  
  ostringstream outstring;
  string name;
  int p1= 30, p2=50; //+- fit limits up to 2 peaks. May be different.
  int const num_par = 5; //number of peaks times 2(pol1)+3(gaus).
  
  for(int ic =0;ic<56;ic++)
    {
      
      outstring.str("");
      outstring << "dEE/dEE_" << ic;
      name = outstring.str();
      
      mycan->cd(1);
      TH2I *hist = (TH2I*)file->Get(name.c_str());
      hist->Draw("col");  
      hist->GetXaxis()->SetRangeUser(200.,1800.); 
      hist->GetYaxis()->SetRangeUser(5.,50.); 
   
      
      if(ic <16 || ic > 31)
	TCutG *mycut = (TCutG*)gates->Get(Form("Zline_%i_2_4",ic));
      else 
	TCutG *mycut = (TCutG*)gates2->Get(Form("Zline_%i_2_4",ic));
      mycut->Draw();

      file->cd();
      outstring.str("");
      outstring << "CsI/CsIGate/ECsI_" << ic << "_Gate";
      name = outstring.str();
      gPad->SetLogz();

      mycan->cd(2);
      TH1I * proj = (TH1I*)file->Get(name.c_str());
      proj->Draw();
      proj->Rebin(4);
      proj->GetXaxis()->SetRangeUser(700.,1800.);

      mycan->Modified();
      mycan->Update();

      TMarker * mark;
      mark=(TMarker*)mycan->WaitPrimitive("TMarker"); //Get the Background limits
      int bkg_lo = mark->GetX();
      delete mark;  
      mark=(TMarker*)mycan->WaitPrimitive("TMarker");
      int bkg_hi = mark->GetX();
      delete mark;
      mark=(TMarker*)mycan->WaitPrimitive("TMarker"); // Get the 1st peak initial guess
      int peak1 = mark->GetX();
      delete mark;
      
      
      double par[num_par] = {0.};
      double out[num_par] = {0.}; 
      int peak1_lo = peak1 - p1, peak1_hi = peak1 + p1; // Peak center and limits
      
      
      TF1 *l1 = new TF1("l1", "pol1", bkg_lo, bkg_hi);
      TF1 *g1 = new TF1("g1", "gaus", peak1_lo,peak1_hi);
      
      TF1 *total = new TF1("total", "pol1(0)+gaus(2)", bkg_lo,bkg_hi);
      
      proj->Fit(l1,"R");
      proj->Fit(g1,"R+");
      
      l1->GetParameters(&par[0]);
      g1->GetParameters(&par[2]);
      
      total->SetParameters(par);
      proj->Fit(total,"R");
      total->GetParameters(out);
      
      
      ofile << ic << " " << out[3] << endl;
      
      outstring.str("");
      outstring << "55A_" << ic;
      name = outstring.str();
      total->SetName(name.c_str());
      total->Draw("same");
      mycan->Modified();
      mycan->Update();
      
      bool IsGood = 0;
//.........这里部分代码省略.........

void StAssociator(Int_t nevents=1,
const char *MainFile="/afs/rhic.bnl.gov/star/data/samples/*.geant.root")
{

  // Dynamically link needed shared libs
  gSystem->Load("St_base");
  gSystem->Load("StChain");
  gSystem->Load("StBFChain");

  gSystem->Load("St_Tables");
  gSystem->Load("StUtilities");
  gSystem->Load("StIOMaker");
  gSystem->Load("StarClassLibrary");
  gSystem->Load("StDetectorDbMaker");


  gSystem->Load("StTpcDb");
  gSystem->Load("StEvent");
//   gSystem->Load("StEventMaker"); //not needed if event.root branch present
  gSystem->Load("StEmcUtil"); 
  gSystem->Load("StEEmcUtil");
  
  gSystem->Load("StMcEvent");
  gSystem->Load("StMcEventMaker");
  gSystem->Load("StAssociationMaker");
  gSystem->Load("StMcAnalysisMaker");
    
  chain = new StChain("StChain"); 
  chain->SetDebug();
   
  // Now we add Makers to the chain...

  StIOMaker* ioMaker = new StIOMaker("IO","r",MainFile,"bfcTree");
  ioMaker->SetDebug();
  ioMaker->SetIOMode("r");
  ioMaker->SetBranch("*",0,"0");                 //deactivate all branches
  ioMaker->SetBranch("geantBranch",0,"r"); //activate geant Branch
  ioMaker->SetBranch("eventBranch",0,"r"); //activate geant Branch
//   ioMaker->SetBranch("dstBranch",0,"r"); //activate Event Branch
//   ioMaker->SetBranch("runcoBranch",0,"r"); //activate runco Branch

  // Note, the title "events" is used in the Association Maker, so don't change it.
  // StEventMaker is not needed for event.root files
//   StEventMaker*       eventReader   = new StEventMaker("events","title");
//   eventReader->doPrintMemoryInfo = kFALSE;
  StMcEventMaker*     mcEventReader = new StMcEventMaker; // Make an instance...
  StAssociationMaker* associator    = new StAssociationMaker;

  // If you need to use L3 TRIGGER uncomment the line:
  // associator->useL3Trigger();
  //associator->SetDebug();
  // For tracks created with the Sti package (ITTF), uncomment the next line:
  //associator->useInTracker();
  // Switch to use the distance or id association.
  associator->useDistanceAssoc();
  //associator->useIdAssoc();
  // Note: useDistanceAssoc and useIdAssoc are mutually exclusive
  // and they set and unset the same flag.
  // The flag will be set by the call done at the end.
  //associator->doPrintMemoryInfo = kTRUE;
  StMcAnalysisMaker*  examples      = new StMcAnalysisMaker;

  // Define the cuts for the Associations

  StMcParameterDB* parameterDB = StMcParameterDB::instance();  
  // TPC
  parameterDB->setXCutTpc(.6); // 6 mm
  parameterDB->setYCutTpc(.6); // 6 mm
  parameterDB->setZCutTpc(.6); // 6 mm
  parameterDB->setReqCommonHitsTpc(3); // Require 3 hits in common for tracks to be associated
  // FTPC
  parameterDB->setRCutFtpc(.3); // 3 mm
  parameterDB->setPhiCutFtpc(5*(3.1415927/180.0)); // 5 degrees
  parameterDB->setReqCommonHitsFtpc(3); // Require 3 hits in common for tracks to be associated
  // SVT
  parameterDB->setXCutSvt(.08); // 800 um
  parameterDB->setYCutSvt(.08); // 800 um
  parameterDB->setZCutSvt(.08); // 800 um
  parameterDB->setReqCommonHitsSvt(1); // Require 1 hits in common for tracks to be associated
    
    
  // now execute the chain member functions

  chain->PrintInfo();
  Int_t initStat = chain->Init(); // This should call the Init() method in ALL makers
  if (initStat) chain->Fatal(initStat, "during Init()");
    
  int istat=0,iev=1;
 EventLoop: if (iev<=nevents && istat!=2) {
   chain->Clear();
   cout << "---------------------- Processing Event : " << iev << " ----------------------" << endl;
   istat = chain->Make(iev); // This should call the Make() method in ALL makers
   if (istat == 2) { cout << "Last  Event Processed. Status = " << istat << endl; }
   if (istat == 3) { cout << "Error Event Processed. Status = " << istat << endl; }
   iev++; goto EventLoop;
 } // Event Loop
  examples->mAssociationCanvas = new TCanvas("mAssociationCanvas", "Histograms",200,10,600,600);
  TCanvas* myCanvas = examples->mAssociationCanvas;
  myCanvas->Divide(2,2);
  
//.........这里部分代码省略.........

void EfficiencyMaps(TString Signal, TString Analysis, TString Channel, bool entries = false) {

  int const nSRs = 12;
  TString SRname[nSRs] = {"SR1_NoJet", "SR1_NoTag", "SR1_Veto", "SR1_Tag", 
			  "SR2_NoJet", "SR2_NoTag", "SR2_Veto", "SR2_Tag", 
			  "SR3_Veto",  "SR3_Tag",   "SR3_Veto_isr", "SR3_Tag_isr"};
  TString SRtitle[nSRs] = {"SR1^{\\text{0Jet}}_{\\text{0Tag}}", "SR1_^{\\text{Jets}}_{\\text{0Tag}}", "SR1_{\\text{0Tag}}", "SR1_{\\text{Tags}}", 
			   "SR2^{\\text{0Jet}}_{\\text{0Tag}}", "SR2_^{\\text{Jets}}_{\\text{0Tag}}", "SR2_{\\text{0Tag}}", "SR2_{\\text{Tags}}", 
			   "SR3_{\\text{0Tag}}", "SR3_{\\text{Tags}}", "SR3^{\\text{ISR}}_{\\text{0Tag}}", "SR3^{\\text{ISR}}_{\\text{Tags}}"};
  bool UseSR[nSRs]; 
  if (Analysis=="_Stop") {
    UseSR[0] = 0; UseSR[1] = 0; UseSR [2] = 1; UseSR [3] = 1;
    UseSR[4] = 0; UseSR[5] = 0; UseSR [6] = 1; UseSR [7] = 1;
    UseSR[8] = 0; UseSR[9] = 0; UseSR[10] = 1; UseSR[11] = 1; 
  } else if (Analysis=="_Chargino") {
    UseSR[0] = 1; UseSR[1] = 1; UseSR [2] = 0; UseSR [3] = 0;
    UseSR[4] = 1; UseSR[5] = 1; UseSR [6] = 0; UseSR [7] = 0;
    UseSR[8] = 1; UseSR[9] = 0; UseSR[10] = 0; UseSR[11] = 0; 
  }

  int nSRsignal = 0;
  for (int sr = 0; sr<nSRs; sr++)
    if (UseSR[sr]) 
      nSRsignal++;

  TH2F *EffMap  = new TH2F("EffMap",  "", 7, 0., 140., nSRsignal, 0., nSRsignal);
  TH2F *StatMap = new TH2F("StatMap", "", 7, 0., 140., nSRsignal, 0., nSRsignal);

  EffMap->SetXTitle("M_{T2}(ll) [GeV]");
  StatMap->SetXTitle("M_{T2}(ll) [GeV]");
  
  TString Syst = entries ? "noweight" : "nominal"; 
  myfile0 = new TFile("../minitrees/rootfilesOct17/" + Syst + "/Stop/" + Signal + ".root", "read");

  int isr = 1;

  for (int sr = 0; sr<nSRs; sr++)
    if (UseSR[sr]) {

      EffMap->GetYaxis()->SetBinLabel(isr, SRtitle[sr]);
      StatMap->GetYaxis()->SetBinLabel(isr, SRtitle[sr]);

      TString HistoName = "Stop/02_" + SRname[sr] + "/h_MT2ll" + Channel; 
      TString HistoNameGen = HistoName; 
      HistoNameGen.ReplaceAll("SR1", "SR1gen");
      HistoNameGen.ReplaceAll("SR2", "SR2gen");
      HistoNameGen.ReplaceAll("SR3", "SR3gen");
      HistoNameGen.ReplaceAll("MT2ll", "MT2llgen");
      if (SRname[sr].Contains("_isr")) {
	HistoName.ReplaceAll("_isr", "");
	HistoName.ReplaceAll("MT2ll", "MT2llisr");
	HistoNameGen.ReplaceAll("_isr", "");
	HistoNameGen.ReplaceAll("MT2ll", "MT2llisr");
      }

      TH1F *ThisHisto = GetHisto(HistoName);
      TH1F *ThisHistoGen = GetHisto(HistoNameGen);

      for (int ib = 1; ib<=7; ib++) {
	
	//float XXX = ThisHisto->GetBinContent(ib);
	//float EEE = ThisHistoGen->GetBinContent(ib);
	float XXX = 35.896*(ThisHisto->GetBinContent(ib)+ThisHistoGen->GetBinContent(ib))/2.;
        float EEE = 35.896*(ThisHisto->GetBinError(ib)+ThisHistoGen->GetBinError(ib))/2.;
	float value = EEE/XXX;
	if (entries) value = ThisHisto->GetBinContent(ib);
	EffMap->SetBinContent(ib, isr, value);
	//EffMap->SetBinError(ib, isr, EEE);
	
      }	

      isr++;
      
    }
  
  gStyle->SetOptStat(0);
  //gStyle->SetPaintTextFormat("5.2f");

  TCanvas *CC = new TCanvas("CC", "", 600, 400);
  CC->Divide(1, 1); 

  EffMap->Draw("COLZtext");
  //EffMap->Draw("COLZtexte");
  
  TString channelFlag = (Channel=="_em") ? "_emu" : "_eemumu";
  TString PlotName = entries ? "SignalEntries" : "SignalStatError";
  CC->Print("../Plots/MapMySignal/" + PlotName + "Map_" + Signal + channelFlag + ".pdf");
  CC->Print("../Plots/MapMySignal/" + PlotName + "Map_" + Signal + channelFlag + ".png");

}

void pileup()
{
//  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HIData_Minbias_2760GeV/merged/HIData_Minbias_2760GeV_PPRereco_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep");
  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/HydjetReco_Minbias_2760GeV/merged/HydjetReco_Minbias_2760GeV_PeripheralNew_INCLMULT_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep");
//  TH3D* h = (TH3D*)GetHist("/net/hisrv0001/home/davidlw/scratch1/DiHadronCorrelations/outputs_312/PAData_Minbias_5TeV/merged/PAData_Minbias_5TeV_HM_PromptReco_INCLMULTAllPURej3Vtx_nmin-1_nmax-1_etatrg-2.4-2.4_etaass-2.4-2.4_centmin-1_centmax-1.root","multmaxvssecvszvtxsep");

  TF1* func[40];
  TGraph* xx_yy[40];
  double par0_arr[40];
  double par1_arr[40];
  double par2_arr[40];
  double dz_arr[40];
//  double binwidth=0.2;
  double binedges[16]={0.0,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.4,1.8,2.4,3.0,4.2,5.0,7.0,10.0};

  TCanvas* c16 = new TCanvas("c16","c16",800,600);
  c16->Divide(4,4);

  for(int i=0;i<15;i++)
  {
//    h->GetZaxis()->SetRange(h->GetZaxis()->FindBin(binwidth*i),h->GetZaxis()->FindBin(binwidth*i+binwidth));
    h->GetZaxis()->SetRange(h->GetZaxis()->FindBin(binedges[i]),h->GetZaxis()->FindBin(binedges[i+1]));
    h->GetYaxis()->SetRange(h->GetYaxis()->FindBin(2),h->GetYaxis()->FindBin(50));
    TH2D* h2D = (TH2D*)h->Project3DProfile(Form("yx%d",i));
    h2D->Rebin2D(10,1);
    TProfile* hprof_x = (TProfile*)h2D->ProfileX("hprof_x",-1,-1,"s");
    double xx[500];
    double yy[500];
    for(int ibin=1;ibin<=hprof_x->GetNbinsX();ibin++)
    {
      yy[ibin-1] = hprof_x->GetBinContent(ibin)+2*hprof_x->GetBinError(ibin);
      xx[ibin-1] = hprof_x->GetBinCenter(ibin);
    }
  
// Hydjet
//    TF1* func = new TF1("func","[0]+[1]*(x-[2])*(x-[2])",0,500);
//    func->SetParameters(20,-0.1,150);
  func[i] = new TF1(Form("func_%d",i),"[0]*x**[1]",0,500);
  func[i]->SetParameters(1,0.5);
  func[i]->SetParLimits(0,0,10000);
  func[i]->SetParLimits(1,0,10000);
  if(i==0) func[i]->FixParameter(0,8); //pPb
//  else if(i>11) func[i]->FixParameter(0,20); //pPb
  else func[i]->FixParameter(0,2.); //pPb

//  func->FixParameter(0,5); //PbPb

//    TF1* func = new TF1("func","[0]+[1]*x",40,500);
//    double par1 = (hprof_x->GetBinContent(6)+hprof_x->GetBinError(6)-hprof_x->GetBinContent(1)-hprof_x->GetBinError(1))/(hprof_x->GetBinCenter(6)-hprof_x->GetBinCenter(1));
//    double par0 = hprof_x->GetBinContent(6)+hprof_x->GetBinError(6)-par1*hprof_x->GetBinCenter(6);
//    func->SetParameters(par0,par1);

    xx_yy[i] = new TGraph(hprof_x->GetNbinsX(),xx,yy);
    xx_yy[i]->Fit(Form("func_%d",i),"NO","",0,300);

    c16->cd(i+1);
    xx_yy[i]->Draw("AP"); 
    func[i]->SetLineColor(kRed);
    func[i]->Draw("Lsame");   

    par2_arr[i] = func[i]->GetParameter(2);    
    par1_arr[i] = func[i]->GetParameter(1);
    par0_arr[i] = func[i]->GetParameter(0);
    dz_arr[i] = (binedges[i]+binedges[i+1])/2;
  }

  TGraph* gr_par0 = new TGraph(15,dz_arr,par0_arr);
  TGraph* gr_par1 = new TGraph(15,dz_arr,par1_arr);
  TGraph* gr_par2 = new TGraph(15,dz_arr,par2_arr);

  TCanvas* cc = new TCanvas("cc","",950,400);  
  cc->Divide(3,1);
  cc->cd(1);
  gr_par0->Draw("AP");
  TF1* func_par0 = new TF1("func_par0","[2]*exp(-x**2/[0])*x**[3]+[1]",0.0,10.0);
  func_par0->SetParameters(0.2,0,5,-1);
  func_par0->FixParameter(2,0);
  gr_par0->Fit("func_par0","NO","",0.3,4.0);
  func_par0->Draw("Lsame");
  cc->cd(2);
  gr_par1->Draw("AP");
//  TF1* func_par1 = new TF1("func_par1","[2]*exp(-x**2/[0])*x**[3]+[1]",0.,10.0);
//  func_par1->SetParameters(0.2,0,-0.0005,0);
//  TF1* func_par1 = new TF1("func_par1","[2]*exp(-x/[0])*x**[3]+[1]",0.,10.0);
  TF1* func_par1 = new TF1("func_par1","[1]*exp(-x**[2]/[0])*x**[3]",0.,10.0);
  func_par1->SetParameters(0.1,0.6,2,1);
//  func_par1->FixParameter(2,1);
//  func_par1->FixParameter(3,1);

  gr_par1->Fit("func_par1","NO","",0.3,4.0);
  func_par1->Draw("Lsame");
  cc->cd(3);
  gr_par2->Draw("AP"); 
  TF1* func_par2 = new TF1("func_par2","[2]*exp(-x**2/[0])*x**[3]+[1]",0.,10.0);
  func_par2->SetParameters(1.2,170,120,1);
  gr_par2->Fit("func_par2","NO","",0.2,2.0);
  func_par2->Draw("Lsame");

cout<<Form("par0=%.2f*exp(-x**2/%.2f)*(x**%.2f)+%.2f",func_par0->GetParameter(2),func_par0->GetParameter(0),func_par0->GetParameter(3),func_par0->GetParameter(1))<<endl;
cout<<Form("par1=%.6f*exp(-x/%.6f)*(x**%.6f)+%.2f",func_par1->GetParameter(2),func_par1->GetParameter(0),func_par1->GetParameter(3),func_par1->GetParameter(1))<<endl;
//.........这里部分代码省略.........

void Unfold2(int algo= 3,bool useSpectraFromFile=0, bool useMatrixFromFile=0, int doToy = 0, int isMC = 0,char *spectraFileName = (char*)"pbpb_spectra_akPu3PF.root",double recoJetPtCut = 60,double trackMaxPtCut = 0, int nBayesianIter = 4, int doBjets=0) // algo 2 =akpu2 ; 3 =akpu3 ; 4 =akpu4 ;1 = icpu5
{
  
  gStyle->SetErrorX(0.5);
  gStyle->SetPaintTextFormat("3.2f");
  gStyle->SetOptLogz(1);
  gStyle->SetPadRightMargin(0.13);	
  gStyle->SetOptTitle(0);

	const float pplumi=5.3e9;
        const bool SavePlot=kTRUE;

  TH1::SetDefaultSumw2();
  TH2::SetDefaultSumw2();
  
  // input files
  char *fileNamePP_data = (char*)"/net/hidsk0001/d00/scratch/maoyx/pPb/Btag/CodeMatt/NewFormatV4_bFractionMCTemplate_pppp1_SSVHEat2.0FixCL0_bin_0_40_eta_0_2.root";
  char *fileNamePbPb_data = (char*)"/net/hidsk0001/d00/scratch/maoyx/pPb/Btag/CodeMatt/AltBinningV6_bFractionMCTemplate_ppPbPb1_SSVHEat2.0FixCL0_bin_0_40_eta_0_2.root";
  char *fileNamePP_mc = NULL;
  if(doBjets) fileNamePP_mc = (char*)"/net/hidsk0001/d00/scratch/kjung/histos/ppMC_ppReco_ak3PF_BjetTrig_noIPupperCut.root";
  else fileNamePP_mc = (char*)"/net/hidsk0001/d00/scratch/kjung/histos//ppMC_ppReco_ak3PF_QCDjetTrig_noIPupperCut.root";
  char *fileNamePbPb_mc = NULL;
  if(doBjets)fileNamePbPb_mc = (char*) "/net/hisrv0001/home/mnguyen/scratch/bTaggingOutput/ntuples/PbPbBMC_pt30by3_ipHICalibCentWeight_noTrig.root";
  else fileNamePbPb_mc = (char*) "/net/hisrv0001/home/mnguyen/scratch/bTaggingOutput/ntuples/PbPbQCDMC_pt30by3_ipHICalibCentWeight_noTrig.root";

  // grab ntuples
  TFile *infPbPb_mc = new TFile(fileNamePbPb_mc);
  TFile *infPP_mc = new TFile(fileNamePP_mc);


  // Output file
  TFile *pbpb_Unfo;
  if (isMC) pbpb_Unfo = new TFile(Form("pbpb_Unfo_%s_MC.root",algoName[algo]),"RECREATE");
  else pbpb_Unfo  = new TFile(Form("pbpb_Unfo_%s_jtpt%.0f_trk%.0f.root",algoName[algo],recoJetPtCut,trackMaxPtCut),"RECREATE");
  // Histograms used by RooUnfold
  UnfoldingHistos *uhist[nbins_cent+1];
		
  // Initialize Histograms   
	
  for (int i=0;i<=nbins_cent;i++) uhist[i] = new UnfoldingHistos(i);
	
  // Initialize reweighting functions
  
  TCut dataSelection;
  TCut dataSelectionPP;
  TCut TriggerSelectionPP;
  TCut TriggerSelectionPbPb80;

  if(doBjets)dataSelection = "weight*(abs(refparton_flavorForB)==5&&abs(jteta)<2)";
  else dataSelection = "weight*(abs(jteta)<2)";
  

  if (isMC) cout<<"This is a MC closure test"<<endl;
  else cout<< "This is a data analysis"<<endl;    		     
	     	
  // Setup jet data branches, basically the jet tree branches are assigned to this object when we loop over the events
	
  JetDataPbPb *dataPbPb   = new JetDataPbPb(fileNamePbPb_mc,(char*)"nt"); // PbPb data	
  JetDataPP *dataPP = new JetDataPP(fileNamePP_mc,(char*)"nt");	// pp data
	
  TFile *fSpectra(0);		
  if (useSpectraFromFile||useMatrixFromFile){
    fSpectra = new TFile(spectraFileName,"read");
  }
  
  // Come back to the output file dir
  pbpb_Unfo->cd();

  // Get Jet spectra from data file
  cout <<"Reading data..."<<endl;

  // This doesn't seem to be relevant for the moment -Matt
  /*	

  TTree *tPbPbJet = (TTree*)infPbPb_mc->Get("nt");
  TTree *tPPJet  = (TTree*)infPP_mc->Get("nt");


  TCanvas * cInput = new TCanvas("cInput","Input",800,400);
  cInput->Divide(2,1);
		
  cout <<"Spectra..."<<endl;	
	
  for (int i=0;i<=nbins_cent;i++){
    cout <<nbins_cent<<endl;
    TCut centCut = Form("bin<%.0f&&bin>=%.0f",boundaries_cent[i+1],boundaries_cent[i]);
    if (useSpectraFromFile) {
      uhist[i]->hMeas = (TH1F*)fSpectra->Get(Form("hMeas_cent%d",i));
    } else {
      if (!isMC) {
	tPbPbJet->Project(Form("hMeas_cent%d",i),"jtptB", dataSelection&&centCut&&TriggerSelectionPbPb80);
      }   
    }
		
    if (useMatrixFromFile) {
      cout <<"Matrix"<<endl;
      uhist[i]->hMatrixFit = (TH2F*) fSpectra->Get(Form("hMatrixFit_cent%d",i));
      uhist[i]->hMeasMatch = (TH1F*)((TH2F*) fSpectra->Get(Form("hMatrixFit_cent%d",i)))->ProjectionY();
      uhist[i]->hMeasMatch->Divide(uhist[i]->hMeas);
    } else {
//.........这里部分代码省略.........

void Draw(const char* pairName, const char* canvasName, std::vector<CorrelationFunction> &plots, bool isIdentical)
{
	unsigned int i, k;
	int canvasNumber1 = -1, canvasNumber2 = 0;
	TLegend *hHKMlegend, *EPOSlegend;
	TCanvas *canv = new TCanvas(canvasName, canvasName, 1300, 2000);
	if(isIdentical)
		canv->Divide(2,4);
	else
		canv->Divide(2,5);

	// 
	// hHKM
	//
	hHKMlegend = new TLegend(legendX1, legendY1, legendX2, legendY2);

	k = 0;
	canv->cd(++++canvasNumber1);
	for(i = 0 ; i < plots.size(); ++i)
		if(plots[i].model.compare("hHKM") == 0)
		{
			DrawSinglePlot(plots[i].cf1D, (std::string(pairName) + std::string(" - hHKM - 1D;q_{inv} (GeV/c);C(q_{inv})")).c_str(), k);
			SetRanges(plots[i].cf1D, canvasName, "1D");
			hHKMlegend->AddEntry(plots[i].cf1D, plots[i].centrality.c_str(),"P");
		}
	if(plots.size() > 0)
	{
		hHKMlegend->Draw();
		hHKMlegend->SetFillColor(0);
	}
			
	k = 0;
	canv->cd(++++canvasNumber1);
	for(i = 0 ; i < plots.size(); ++i)
		if(plots[i].model.compare("hHKM") == 0)
		{
			DrawSinglePlot(plots[i].cfSH00, (std::string(pairName) + std::string(" - hHKM - Spherical Harmonics;k* (GeV/c);#RgothicC^{0}_{0}")).c_str(), k);
			SetRanges(plots[i].cfSH00, canvasName, "SH00");
		}

	if(!isIdentical)
	{
		k = 0;
		canv->cd(++++canvasNumber1);
		for(i = 0 ; i < plots.size(); ++i)
			if(plots[i].model.compare("hHKM") == 0)
			{
				DrawSinglePlot(plots[i].cfSH11, (std::string(pairName) + std::string(" - hHKM - Spherical Harmonics;k* (GeV/c);#RgothicC^{1}_{1}")).c_str(), k);
				SetRanges(plots[i].cfSH11, canvasName, "SH11");
			}
	}

	k = 0;
	canv->cd(++++canvasNumber1);
	for(i = 0 ; i < plots.size(); ++i)
		if(plots[i].model.compare("hHKM") == 0)
		{
			DrawSinglePlot(plots[i].cfSH20, (std::string(pairName) + std::string(" - hHKM - Spherical Harmonics;k* (GeV/c);#RgothicC^{0}_{2}")).c_str(), k);
			SetRanges(plots[i].cfSH20, canvasName, "SH20");
		}

	k = 0;
	canv->cd(++++canvasNumber1);
	for(i = 0 ; i < plots.size(); ++i)
		if(plots[i].model.compare("hHKM") == 0)
		{
			DrawSinglePlot(plots[i].cfSH22, (std::string(pairName) + std::string(" - hHKM - Spherical Harmonics;k* (GeV/c);#RgothicC^{2}_{2}")).c_str(), k);
			SetRanges(plots[i].cfSH22, canvasName, "SH22");
		}

	//
	// EPOS
	//
	EPOSlegend = new TLegend(legendX1, legendY1, legendX2, legendY2);

	k = 0;
	canv->cd(++++canvasNumber2);
	for(i = 0 ; i < plots.size(); ++i)
		if(plots[i].model.compare("EPOS") == 0)
		{
			DrawSinglePlot(plots[i].cf1D, (std::string(pairName) + std::string(" - EPOS - 1D;q_{inv} (GeV/c);C(q_{inv})")).c_str(), k);
			SetRanges(plots[i].cf1D, canvasName, "1D");
			EPOSlegend->AddEntry(plots[i].cf1D, plots[i].centrality.c_str(),"P");
		}
	if(plots.size() > 0)
	{
		EPOSlegend->Draw();
		EPOSlegend->SetFillColor(0);
	}
		
	k = 0;
	canv->cd(++++canvasNumber2);
	for(i = 0 ; i < plots.size(); ++i)
		if(plots[i].model.compare("EPOS") == 0)
		{
			DrawSinglePlot(plots[i].cfSH00, (std::string(pairName) + std::string(" - EPOS - Spherical Harmonics;k* (GeV/c);#RgothicC^{0}_{0}")).c_str(), k);
			SetRanges(plots[i].cfSH00, canvasName, "SH00");
		}

	if(!isIdentical)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      if (factor / background < 10) continue;

      // Get the total number of events in the gaussian 1 sigma
      Int_t sigRangeLow = hist->FindBin(mean - sigma);
      Int_t sigRangeHigh = hist->FindBin(mean + sigma);

      double sigRangeTotal = 0;
      for (Int_t bin = sigRangeLow; bin <= sigRangeHigh; bin++)
        sigRangeTotal += hist->GetBinContent(bin);

      // 2 * 1 sigma integral shoudl give ~ area under gaussian
      sigma /= sqrt(2 * sigRangeTotal);

      ptsX.push_back(hist->GetYaxis()->GetBinCenter(row));
      ptsY.push_back(mean);
      ptsErr.push_back(sigma);
    }

    if (ptsX.size() < 3) continue;

    std::vector<double> yvals = ptsY;
    std::sort(yvals.begin(), yvals.end());
    const double median = yvals[yvals.size()/2];
    double avgDeviation = 0;
    for (unsigned int i = 0; i < yvals.size(); i++)
      avgDeviation += fabs(yvals[i] - median);
    avgDeviation /= (double)yvals.size();

    std::vector<double> ptsXGood;
    std::vector<double> ptsYGood;
    std::vector<double> ptsErrGood;

    for (unsigned int i = 0; i < ptsX.size(); i++)
    {
      if (fabs(ptsY[i] - median) > 1.5*avgDeviation) continue;
      ptsXGood.push_back(ptsX[i]);
      ptsYGood.push_back(ptsY[i]);
      ptsErrGood.push_back(ptsErr[i]);
    }

    if (ptsXGood.size() < 3) continue;

    TGraphErrors* graph = new TGraphErrors(ptsXGood.size(),
                                           &(ptsXGood.at(0)),
                                           &(ptsYGood.at(0)), 0,
                                           &(ptsErrGood.at(0)));

    TF1* fitFunc = new TF1("f1", "1 ++ x");
    TF1* result = 0;

    graph->Fit(fitFunc, "Q0E").Get();
    result = graph->GetFunction(fitFunc->GetName());

    // Weight the angle by the slope uncertainty and the inverse of the chi2 normalized
    double weight = result->GetParError(1);
    const double chi2 = result->GetChisquare() / (double)result->GetNDF();
    fitChi2 += chi2;
    weight *= chi2;
    if (weight > 10 * DBL_MIN) weight = 1.0 / weight;
    else weight = 1.0;

    if (axis)
    {
      rotation -= weight * atan(result->GetParameter(1));
      offsetX = result->GetParameter(0);
    }
    else
    {
      rotation += weight * atan(result->GetParameter(1));
      offsetY = result->GetParameter(0);
    }

    angleWeights += weight;

    if (display)
    {
      TCanvas* can = new TCanvas("ResidualAlignment", "Residual Alignment", 900, 600);
      can->Divide(2);
      can->cd(1);
      hist->Draw("COLZ");
      can->cd(2);
      result->SetLineColor(46);
      result->SetLineWidth(2);
      graph->Draw("ap");
      result->Draw("SAME");
      can->Update();
      can->WaitPrimitive();
    }

    delete fitFunc;
    delete graph;
  }

  if (angleWeights > 10 * DBL_MIN)
    rotation /= angleWeights;
  std::cout << "relaxation: " << relaxation << std::endl;
  rotation *= relaxation;
  offsetX *= relaxation;
  offsetY *= relaxation;
}

//=============================================================================
// Standard constructor, initializes variables
//=============================================================================
void makeLongTracksonly2(  ) {
  gROOT->SetBatch(kTRUE);
  TString tuplelocation = "/afs/cern.ch/work/m/mwilkins/b_b-bar_cross-section/";
  TString data2011tuple = tuplelocation+"Strp20r1_SL_D0andDp_MD.root";
  TString MC2011tuple = tuplelocation+"MC_2011_Bu_D0Xmunu_cocktail_12873441_MC2011_S20r1_noPID_Tuples.root";
  TString data2015tuple = tuplelocation+"B2DMuNuX_tuples_05082015.root";
  TString MC2015tuple = tuplelocation+"MC_2015_Bu_D0Xmunu_cocktail_12873441_MC2015_S22_noPID_Tuples.root";
  TFile *f[4];
  f[0] = TFile::Open(data2011tuple);
  f[1] = TFile::Open(MC2011tuple);
  f[2] = TFile::Open(data2015tuple);
  f[3] = TFile::Open(MC2015tuple);
  cout<<"files opened"<<endl;
  
  const int nhpc = 2; //number of histograms per canvas
  const int ncanvases = 4;
  const int nstacks = ncanvases;
  THStack *hs[nstacks];
  TCanvas *c[ncanvases];
  TH1F *h[ncanvases][nhpc];
  TLegend *leg[ncanvases];
  TString outputlocation = "output/longtracksplots2/";//where the output will be saved
  TString plotfilename = "longtracksplots2.pdf";//name of file that will hold drawn plots

  TCanvas *cf = new TCanvas("cf","combined");//canvas to hold everything
  float sqnc = sqrt(ncanvases), sqncu = ceil(sqnc), sqncd = floor(sqnc);
  while(sqncu*sqncd<ncanvases) sqncu++;
  cf->Divide(sqncu,sqncd);
  TLegend *legf = new TLegend(0.84, 0.84, .99, .95);

  TString placeholder; //to avoid adding strings within functions; assign immediately before use
  TString placeholder2;
  
  cout<<"starting canvas loop..."<<endl;
  for(int ci =0; ci<ncanvases; ci++){ //loop over canvases
    //create indicators and strings:
    int year;
    if(ci<(ncanvases/2)){//2011 for 1st half canvases, 2015 for 2nd half
      year = 2011;
    }else year = 2015;
    TString yearstring = Form("%d",year);
    TString filetype;
    if(ci%2==0){//true every other canvas
      filetype = "data";
    }else filetype = "MC";
    int file_num=-1000;//assign file_num based on combinations of data and MC
    if(year==2011){
      if(filetype=="data") file_num=0;
      if(filetype=="MC") file_num=1;
    }
    if(year==2015){
      if(filetype=="data") file_num=2;
      if(filetype=="MC") file_num=3;
    }
    int B0;//these will be assigned within the histogram loop
    TString Bdecay;
    TString branch;
    TString  tracktype = "nLongTracks";
    TString cistring = Form("%d",ci);
    cout<<"loop "<<ci<<" indicators and strings created:"<<endl<<"year: "<<year<<endl;
    cout<<"file: "<<file_num<<endl;

    //create the plots
    placeholder = "c"+cistring;
    placeholder2=yearstring+filetype+tracktype;
    c[ci] = new TCanvas(placeholder,placeholder2,1200,800); //create the canvases
    c[ci]->cd();
    gStyle->SetOptStat("");
    leg[ci] = new TLegend(0.7, 0.7, .97, .93);//create legend
    placeholder = "hs"+cistring;
    hs[ci] = new THStack(placeholder,placeholder2); //create the stack to hold the histograms
    cout<<"starting histogram loop..."<<endl;
    for(int hi=0;hi<nhpc;hi++){ //loop over histograms in a canvas; B0=0 and B0=1
      if((int)floor(hi)%2 ==0){//every other histogram
        B0=0;
        Bdecay = "B^{-}->(D^{0}->K^{-} #pi^{+})#mu^{-}";
        if((filetype=="data")&&(year==2011)){
          branch="tupleb2D0Mu/tupleb2D0Mu";
        }else branch="Tuple_b2D0MuX/DecayTree";
      } else{
        B0=1;
        Bdecay = "B^{0}->(D^{-}->K^{+} #pi^{-} #pi^{-})#mu^{+}";
        if((filetype=="data")&&(year==2011)){
          branch="tupleb2DpMu/tupleb2DpMu";
        }else branch="Tuple_b2DpMuX/DecayTree";
      }
      //create convenient strings
      TString histring = Form("%d",hi);
      TString hname = "h"+cistring+histring;
      
      //create histograms
      if(tracktype=="nTracks") h[ci][hi] = new TH1F(hname,tracktype,104,0,1144);
      if(tracktype=="nLongTracks") h[ci][hi] = new TH1F(hname,tracktype,131,0,262);
      cout<<"histogram loop "<<hi<<" strings and histograms created"<<endl;

      //navigate files
      cout<<"navigating file..."<<endl;
//.........这里部分代码省略.........

void higgsAna::Plot(string histoName)
{
	string temp1 = "hee_" + histoName;
	string temp2 = "hmumu_" + histoName;
	string temp3 = histoName;
	char *eeName, *mumuName, *Name;
	eeName = new char [temp1.size()+1];
	mumuName = new char [temp2.size()+1];
	Name = new char [temp3.size()];
	strcpy (eeName, temp1.c_str());
	strcpy (mumuName, temp2.c_str());
	strcpy (Name, temp3.c_str());
	

	
	isData = 1;
	leptonType = 11;
	SetFiles();
	TFile *eeDataFile = new TFile(outputFileName, "READ");
	TH1F *hee_data = (TH1F*) eeDataFile->Get(eeName);
	hee_data->SetMarkerStyle(20);
	hee_data->SetMarkerColor(1);
	hee_data->SetLineColor(1);
	hee_data->SetMarkerSize(1);
	
	leptonType = 13;
	SetFiles();
	TFile *mumuDataFile = new TFile(outputFileName, "READ");
	TH1F *hmumu_data = (TH1F*) mumuDataFile->Get(mumuName);
	hmumu_data->SetMarkerStyle(20);
	hmumu_data->SetMarkerColor(1);
	hmumu_data->SetLineColor(1);
	hmumu_data->SetMarkerSize(1);

	TCanvas *can = new TCanvas(Name, Name, 0, 0, 840, 800);
	can->Divide(2,2);
	
	can->cd(1);
	hee_data->DrawCopy("e");
	can->cd(3);
	hmumu_data->DrawCopy("e");
	can->cd(2);
	can->cd(2)->SetLogy();
	hee_data->SetMinimum(0.0009);
	hee_data->DrawCopy("e");
	can->cd(4);
	can->cd(4)->SetLogy();
	hmumu_data->SetMinimum(0.0009);
	hmumu_data->DrawCopy("e");
	

	TFile *MCFile[20];
	TH1F *hee_MC[20];
	TH1F *hmumu_MC[20];

	isData = 0;
	
	cout << endl << "SCALING" << endl;
	for (int i = 0; i < 20; i++){
		SetFiles(i+11);
		MCFile[i] = new TFile(outputFileName, "READ");
		hee_MC[i] = (TH1F*) MCFile[i]->Get(eeName);
		hee_MC[i]->Scale(eeScaling[i]);	
		hmumu_MC[i] = (TH1F*) MCFile[i]->Get(mumuName);
		hmumu_MC[i]->Scale(mumuScaling[i]);
		
		
		//cout <<" ee:  " << eeScaling[i] << "   mumu: " << mumuScaling[i] << endl;
		
	
	}

	// add histograms

	Color_t colors[20] = {kGreen-10,0,0,0,0,0,0,0,0, kRed+2, kGreen-6, kGreen+1,0,0,0,0,0, kAzure-6, kBlue-6, kBlue-10}; 	
	//Int_t colors[20] = {30,0,0,0,0,0,0,0,0, 50, 72, 8,0,0,0,0,0, 54, 51, 38}; 
	//Int_t colors[8] = {50, 66, 91, 0, 0, 59, 51, 8}; 

	// DY + WJets + Top + WW + WZ + ZZ
	// WJets + Top + WW + WZ + ZZ
	// Top + WW + WZ + ZZ
	// WW + WZ + ZZ
	// WZ + ZZ
	// ZZ
	for (int j = 0; j < 20; j++){
		if ( (j >=1 && j<=8) || (j >= 12 && j <= 16) ) continue;
		for (int i = j+1; i < 20; i++){
			hee_MC[j]->Add(hee_MC[i]);
			hmumu_MC[j]->Add(hmumu_MC[i]);
		}
		can->cd(1);
		hee_MC[j]->SetFillColor(colors[j]);
		hee_MC[j]->DrawCopy("same");
		can->cd(2);
		hee_MC[j]->SetMinimum(0.0009);
		hee_MC[j]->DrawCopy("same");
		can->cd(3);
		hmumu_MC[j]->SetFillColor(colors[j]);
		hmumu_MC[j]->DrawCopy("same");
		can->cd(4);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	RooPlot* frame2 = mcstudy->plotError(Ntt,Bins(40)) ;
	RooPlot* frame3 = mcstudy->plotPull(Ntt,Bins(40),FitGauss(kTRUE)) ;
	RooPlot* frame1w = mcstudy->plotParam(Nbkg,Bins(40)) ;
	RooPlot* frame2w = mcstudy->plotError(Nbkg,Bins(40)) ;
	RooPlot* frame3w = mcstudy->plotPull(Nbkg,Bins(40),FitGauss(kTRUE)) ;
	RooPlot* frame1st = mcstudy->plotParam(NST,Bins(40)) ;
	RooPlot* frame2st = mcstudy->plotError(NST,Bins(40)) ;
	//RooPlot* frame3st = mcstudy->plotPull(NST,Bins(40),FitGauss(kTRUE)) ;
	
	// Plot distribution of minimized likelihood
	RooPlot* frame4 = mcstudy->plotNLL(Bins(40)) ;

	// Make some histograms from the parameter dataset
	TH1* hh_cor_ttbar_w = mcstudy->fitParDataSet().createHistogram("hh",Ntt,YVar(Nbkg)) ;

	// Access some of the saved fit results from individual toys
	//TH2* corrHist000 = mcstudy->fitResult(0)->correlationHist("c000") ;
	//TH2* corrHist127 = mcstudy->fitResult(127)->correlationHist("c127") ;
	//TH2* corrHist953 = mcstudy->fitResult(953)->correlationHist("c953") ;

	
	// Draw all plots on a canvas
	gStyle->SetPalette(1) ;
	gStyle->SetOptStat(0) ;

	TCanvas* cv = new TCanvas("cv","cv",600,600) ;
	hM3->SetFillColor(kRed);
	hWjets->SetFillColor(kGreen);
	hM3->Draw();
	hWjets->Draw("same");
	gPad->RedrawAxis();
	
	TCanvas* cva = new TCanvas("cva","cva",1800,600) ;
	cva->Divide(3);
	cva->cd(1) ;
	RooPlot *initialframe = mass.frame();
	//initial->SetMaximum(10);
	hpdf_ttbar.plotOn(initialframe,LineColor(kRed));
	hpdf_wjetsFast.plotOn(initialframe,LineColor(kGreen));
	hpdf_ST.plotOn(initialframe,LineColor(kYellow));
	initialframe->Draw();
	//initialframe->SetTitle();
	cva->cd(2);
	//retrieve data for only one PE
	
	int Npe = 10;
	RooPlot *genframe = mass.frame(Nbins);
	RooDataSet *gendata = mcstudy->genData(Npe);
	cout << " N events = " << gendata->numEntries() << endl;
	gendata->plotOn(genframe);
	//mcstudy->fitResult(Npe)->plotOn(genframe, model_histpdf);
	genframe->Draw();
	cva->cd(3);
	RooPlot *genframe2 = mass.frame(Nbins);
	mcstudy->fitResult(Npe)->Print("v");
	gendata->plotOn(genframe2);
	RooArgList arglist = mcstudy->fitResult(Npe)->floatParsFinal();

	
	//cout << "name of argument:" << arglist[2].GetName() << endl;
	//cout << "name of argument:" << arglist[1].GetName() << endl;
	//cout << "name of argument:" << arglist[0].GetName() << endl;
	
	RooAddPdf model_histpdf_fitted("modelfitted", "TTjets+Wjets", RooArgList(hpdf_ttbar,hpdf_wjetsFast,hpdf_ST),
							   RooArgList(arglist[2],arglist[1],arglist[0]) ) ;
	

void EcalRecHitsPlotCompare( TString currentfile = "EcalRecHitsValidation_new.root",
                             TString referencefile = "EcalRecHitsValidation_old.root")
{

    gROOT ->Reset();
    char*  rfilename = referencefile;
    char*  sfilename = currentfile;

    int rcolor = 2;
    int scolor = 4;

    delete gROOT->GetListOfFiles()->FindObject(rfilename);
    delete gROOT->GetListOfFiles()->FindObject(sfilename);

    TText* te = new TText();
    te->SetTextSize(0.1);
    TFile * rfile = new TFile(rfilename);
    TFile * sfile = new TFile(sfilename);

    rfile->cd("DQMData/EcalRecHitsV/EcalRecHitsTask");
    gDirectory->ls();

    sfile->cd("DQMData/EcalRecHitsV/EcalRecHitsTask");
    gDirectory->ls();

    Char_t histo[200];

    gStyle->SetOptStat("n");

    gROOT->ProcessLine(".x HistoCompare.C");
    HistoCompare * myPV = new HistoCompare();

//////////////////////////////////////////////////////////////

// General class: Particle Gun
    if (1)
    {
        TH1* meGunEnergy_;
        rfile->GetObject("DQMData/EcalRecHitsV/EcalRecHitsTask/EcalRecHitsTask Gun Momentum;1",meGunEnergy_);
        meGunEnergy_;
        meGunEnergy_->SetLineColor(rcolor);

        TH1* meGunEta_;
        rfile->GetObject("DQMData/EcalRecHitsV/EcalRecHitsTask/EcalRecHitsTask Gun Eta;1",meGunEta_);
        meGunEta_;
        meGunEta_->SetLineColor(rcolor);

        TH1* meGunPhi_;
        rfile->GetObject("DQMData/EcalRecHitsV/EcalRecHitsTask/EcalRecHitsTask Gun Phi;1",meGunPhi_);
        meGunPhi_;
        meGunPhi_->SetLineColor(rcolor);

        TH1* newmeGunEnergy_;
        sfile->GetObject("DQMData/EcalRecHitsV/EcalRecHitsTask/EcalRecHitsTask Gun Momentum;1",newmeGunEnergy_);
        newmeGunEnergy_;
        newmeGunEnergy_->SetLineColor(scolor);

        TH1* newmeGunEta_;
        sfile->GetObject("DQMData/EcalRecHitsV/EcalRecHitsTask/EcalRecHitsTask Gun Eta;1",newmeGunEta_);
        newmeGunEta_;
        newmeGunEta_->SetLineColor(scolor);

        TH1* newmeGunPhi_;
        sfile->GetObject("DQMData/EcalRecHitsV/EcalRecHitsTask/EcalRecHitsTask Gun Phi;1",newmeGunPhi_);
        newmeGunPhi_;
        newmeGunPhi_->SetLineColor(scolor);


        // --------------------------------
        TCanvas * Ecal = new TCanvas("Ecal","Ecal",800,1000);
        Ecal->Divide(1,3);

        Ecal->cd(1);
        if ( meGunEnergy_ && newmeGunEnergy_ ) {
            meGunEnergy_->Draw();
            newmeGunEnergy_->Draw("same");
            myPV->PVCompute( meGunEnergy_ , newmeGunEnergy_ , te );
        }

        Ecal->cd(2);
        if ( meGunEta_ && newmeGunEta_ ) {
            meGunEta_->Draw();
            newmeGunEta_->Draw("same");
            myPV->PVCompute( meGunEta_ , newmeGunEta_ , te );
        }

        Ecal->cd(3);
        if ( meGunPhi_ && newmeGunPhi_ ) {
            meGunPhi_->Draw();
            newmeGunPhi_->Draw("same");
            myPV->PVCompute( meGunPhi_ , newmeGunPhi_ , te );
        }

        Ecal->Print("ParticleGun_compare.eps");
        delete Ecal;
    }



// General class: sim/rec hit ratios
//.........这里部分代码省略.........

void FakePhotonSystematic(){
  
  TFile FMCRSSig("DsTolnu/SP1005Skim/DsReco/Final/HistosForFit.root");
  TFile FMCWSSig("DsTolnu/SP1005Skim/DsReco/FinalWS/HistosForFit.root");

  //TFile FMCRS("DsTolnu/MC/Final/HistosForFit.root");
  //TFile FMCWS("DsTolnu/MC/FinalWS/HistosForFit.root");
  TFile FMCRS("DsTolnu/SP1005Skim/DsReco/Final/HistosForFit.root");
  TFile FMCWS("DsTolnu/SP1005Skim/DsReco/FinalWS/HistosForFit.root");

  TFile FDataRS("DsTolnu/Data/Final/HistosForFit.root");
  TFile FDataWS("DsTolnu/Data/FinalWS/HistosForFit.root");

  

  TCanvas C;
  C.Clear();
  C.Print("FakePhotonSystematic.ps[");

  C.Clear();
  C.Divide(2,2);
  C.cd(1);  
//   TH1F* HMCRS=(TH1F*)FMCRS.Get("HEvtMultiplicity");
//   HMCRS->SetTitle(TString("MC RS =")+(long)HMCRS->GetMean()+"."+(long)(1000*(HMCRS->GetMean()-1.)));
//   HMCRS->Draw();
//   C.cd(2);
//   TH1F* HMCWS=(TH1F*)FMCWS.Get("HEvtMultiplicity");
//   HMCWS->SetTitle(TString("MC WS =")+(long)HMCWS->GetMean()+"."+(long)(1000*(HMCWS->GetMean()-1.)));
//   HMCWS->Draw();
  C.cd(3);  
  TH1F* HDataRS=(TH1F*)FDataRS.Get("HEvtMultiplicity");
  HDataRS->SetTitle(TString("Data RS =")+(long)HDataRS->GetMean()+"."+(long)(1000*(HDataRS->GetMean()-1.)));
  HDataRS->Draw();
  C.cd(4);
  TH1F* HDataWS=(TH1F*)FDataWS.Get("HEvtMultiplicity");
  HDataWS->SetTitle(TString("Data WS =")+(long)HDataWS->GetMean()+"."+(long)(1000*(HDataWS->GetMean()-1.)));
  HDataWS->Draw();
  C.Print("FakePhotonSystematic.ps");



  C.Clear();
  C.Divide(2,2);
  C.cd(1);  
  TH1F* HSigMCRS=(TH1F*)FMCRSSig.Get("HEvtMultiplicitySig");
  HSigMCRS->SetTitle(TString("MC Sig RS =")+(long)HSigMCRS->GetMean()+"."+(long)(1000*(HSigMCRS->GetMean()-1.)));
  HSigMCRS->Draw();
  C.cd(2);
  TH1F* HSigMCWS=(TH1F*)FMCWSSig.Get("HEvtMultiplicitySig");
  HSigMCWS->SetTitle(TString("MC Sig WS =")+(long)HSigMCWS->GetMean()+"."+(long)(1000*(HSigMCWS->GetMean()-1.)));
  HSigMCWS->Draw();
  C.cd(3);  
  TH1F* HBkgMCRS=(TH1F*)FMCRS.Get("HEvtMultiplicityBkg");
  HBkgMCRS->SetTitle(TString("MC Bkg RS =")+(long)HBkgMCRS->GetMean()+"."+(long)(1000*(HBkgMCRS->GetMean()-1.)));
  HBkgMCRS->Draw();
  C.cd(4);
  TH1F* HBkgMCWS=(TH1F*)FMCWS.Get("HEvtMultiplicityBkg");
  HBkgMCWS->SetTitle(TString("MC Bkg WS =")+(long)HBkgMCWS->GetMean()+"."+(long)(1000*(HBkgMCWS->GetMean()-1.)));
  HBkgMCWS->Draw();
  C.Print("FakePhotonSystematic.ps");
  
  C.Clear();
  THStack HMCStack("HMCStack","");
  HBkgMCRS->SetFillColor(1);HMCStack.Add(HBkgMCRS);
  HSigMCRS->SetFillColor(2);HSigMCRS->SetLineColor(2);HMCStack.Add(HSigMCRS);
  HMCStack.Draw();
  C.Print("FakePhotonSystematic.ps");

  //---------------------------------------------
  //subtract background from Data
  //---------------------------------------------
  TH1F*HDataRSBkg=(TH1F*)HDataWS->Clone("HDataRSBkg");
  TH1F*HBkgMCWSCl=(TH1F*)HBkgMCWS->Clone("HBkgMCWSCl");
  //need to fix the WS first by subtracting the signal component
  HSigMCWS->Scale(11350./HSigMCWS->Integral());//scale to number of WS signal events in Data
  HDataRSBkg->Add(HSigMCWS,-1);
  //construct correction ratio from MC to correct the WS distribution
  TH1F*HBkgMCRSCl=(TH1F*)HBkgMCRS->Clone("HBkgMCRSCl");
  HBkgMCWSCl->Scale(1./HBkgMCWSCl->Integral());
  HBkgMCRSCl->Scale(1./HBkgMCRSCl->Integral());
  TH1F*HBkgMCRatio=(TH1F*)HBkgMCRSCl->Clone("HBkgMCRatio");
  HBkgMCRatio->Divide(HBkgMCWSCl);
  C.Clear();
  delete HBkgMCRatio->FindObject("stats");
  HBkgMCRatio->SetFillColor(0);HBkgMCRatio->SetLineColor(1);
  HBkgMCRatio->SetStats(0);
  HBkgMCRatio->GetYaxis()->SetRangeUser(0,1.5);
  HBkgMCRatio->Draw("hist");
  C.Print("FakePhotonSystematic.ps");
  //correct the WS distribution
  HDataRSBkg->Multiply(HBkgMCRatio);
  C.Clear();
  HDataWS->Scale(1./HDataWS->Integral());
  HDataRSBkg->Scale(1./HDataRSBkg->Integral());
  delete HDataWS->FindObject("stats");
  HDataWS->SetTitle("");
  HDataWS->GetXaxis()->SetTitle("Reconstruction Multiplicity");
  HDataWS->Draw("p");
  HDataRSBkg->SetLineColor(0);
  HDataRSBkg->SetFillColor(4);
//.........这里部分代码省略.........