C++ TCut::GetTitle方法代码示例

void Muti_Plot_centScale(std::string mu_title, const char **lineKind_Arr, TCut selectionC , std::string Var ,  const int ngr, const int nVarBins, double * VarBins_array , double **yArr2, double **yErrArr2)
{ 

	char filename[]="JetESR_status.txt";
	fstream myfile;
	myfile.open (filename,ios::out|ios::app);

	double VarBinsCenter[nVarBins];
	double VarBinsErr[nVarBins];
	for (int ibin=0; ibin<nVarBins; ibin++)
	{
		VarBinsCenter[ibin]= (VarBins_array[ibin] + VarBins_array[ibin+1]) /2;
		//		VarBinsErr[ibin]=0;
		VarBinsErr[ibin]=(VarBins_array[ibin+1]- VarBins_array[ibin])/2;
	}

	cout<<"lineKind_Arr 1 = "<<lineKind_Arr[0]<<" ,lineKind_Arr2 = "<<lineKind_Arr[1]<<endl;
	cout<<"mu_title = "<<mu_title<<" selectionC = "<< selectionC<<endl;

	Can_result[counter1] = new TCanvas(Form("Can_result_%i",counter1));
	mutiGR[counter1] = new TMultiGraph("name",Form("%s",mu_title.c_str()));
	TGraphErrors *gr_mean[ngr];
	legend[counter1] = new TLegend(0.68,0.7,0.88,0.85);

	for(int igr=0; igr<ngr; igr++)
	{
		gr_mean[igr] = new TGraphErrors(nVarBins, VarBinsCenter, yArr2[igr], VarBinsErr , yErrArr2[igr]);
		gr_mean[igr]->SetTitle(Form("some title_%i", igr));
		gr_mean[igr]->SetFillStyle(0);
		gr_mean[igr]->SetMarkerStyle(1);
		gr_mean[igr]->SetLineColor(igr+1);	
		gr_mean[igr]->SetMarkerColor(igr+1);
		mutiGR[counter1]->Add(gr_mean[igr]);
		legend[counter1]->AddEntry( gr_mean[igr] ,lineKind_Arr[igr],"l");  // not work for unknow reason

	}

	myfile<<"hello, Var = "<<Var<<"Var.compare(Centrality) = "<<Var.compare("Centrality")<<endl;

	mutiGR[counter1]->Draw("AP");
	mutiGR[counter1]->GetXaxis()->SetTitle("Centrality"); // must after draw to create a vitual histogram like object to set title.
	mutiGR[counter1]->GetYaxis()->SetTitle("#mu_{Reco./Gen.}");
	mutiGR[counter1]->Draw("AP"); // draw again to let all the settting on it.
	//	Can_result[counter1]->BuildLegend(0.7,0.7,0.9,0.9);
	legend[counter1]->SetBorderSize(0);
	legend[counter1]->Draw();

	TLatex *tex = new TLatex(0.25,0.7,"#bf{#color[2]{|#eta_{jet}|<2.0 &75<refpt<120}}");
	tex->SetNDC();
	tex->SetTextFont(43);
	tex->SetTextSize(24);
	tex->SetLineWidth(2);
	tex->Draw();

	Can_result[counter1]->SaveAs(Form("./Plots/result_%s_%s_%s.pdf",mu_title.c_str(),selectionC.GetTitle() , Var.c_str())); 
	Can_result[counter1]->SaveAs(Form("./Plots/result_%s_%s_%s.C",mu_title.c_str(),selectionC.GetTitle() , Var.c_str()));
	counter1++;
	cout<<"finish mutigraph"<<endl;
}

void make_psi2s_mupt(float ptmin=0.0, float ptmax=30.0, float ymin=0.0, float ymax=2.4, bool absRapidity=true, bool saveFile=false)
{
  TCanvas *c1 = new TCanvas("c1","c1");
  c1->SetLogy();
  TCut defaultCut = "Reco_QQ_sign==0&&Reco_QQ_4mom.M()>3.35&&Reco_QQ_4mom.M()<4.0&&Reco_QQ_ctauTrue>-10";
  TCut ptCut = Form("Reco_QQ_4mom.Pt()>%4.1f&&Reco_QQ_4mom.Pt()<%4.1f",ptmin,ptmax);
  TCut rapCut;
  if (absRapidity)
    rapCut = Form("abs(Reco_QQ_4mom.Rapidity())>%3.1f&&abs(Reco_QQ_4mom.Rapidity())<%3.1f",ymin,ymax);
  else
    rapCut = Form("Reco_QQ_4mom.Rapidity()>%3.1f&&Reco_QQ_4mom.Rapidity()<%3.1f",ymin,ymax);

  unsigned int trigBit=2; // DoubleMu0_HighQ
  TCut trigCut = Form("(HLTriggers&%u)==%u&&(Reco_QQ_trig&%u)==%u",trigBit,trigBit,trigBit,trigBit);

  TString fname;
  fname = Form("20140324/MC_psi2s_pp_mupt_Rap_%3.1f-%3.1f_Pt_%3.1f-%3.1f.root",ymin,ymax,ptmin,ptmax);

  std::cout << fname << std::endl;


  std::cout << "default: " << defaultCut.GetTitle() << std::endl;
  std::cout << "pt cut: " << ptCut.GetTitle() << std::endl;
  std::cout << "rapidity cut: " << rapCut.GetTitle() << std::endl;
  std::cout << "trigger bit: " << trigCut.GetTitle() << std::endl;

  TH1F *hMuPlPtRec = new TH1F("hMuPlPtRec","hMuPlPtRec;p_{T} (#mu^{+}) (GeV/c);Events",200,0,20);
  TH1F *hMuMiPtRec = new TH1F("hMuMiPtRec","hMuMiPtRec;p_{T} (#mu^{-}) (GeV/c);Events",200,0,20);
  TH1F *hMuPtRec = new TH1F("hMuPtRec","hMuPtRec;p_{T} (#mu^{#pm}) (GeV/c);Events",200,0,20);
  hMuPlPtRec->Sumw2();
  hMuMiPtRec->Sumw2();
  hMuPtRec->Sumw2();

  hMuPlPtRec->SetMarkerColor(kRed);
  hMuMiPtRec->SetMarkerColor(kBlue);

  TChain *myTree = new TChain("myTree");
  myTree->Add("../root_files/PRpsi2SMC_Histos_2013pp_GlbGlb_STARTHI53_V28-v1_muLessPV.root");

  myTree->Draw("Reco_QQ_mupl_4mom.Pt()>>hMuPlPtRec",defaultCut&&ptCut&&rapCut&&trigCut,"e");
  myTree->Draw("Reco_QQ_mumi_4mom.Pt()>>hMuMiPtRec",defaultCut&&ptCut&&rapCut&&trigCut,"e");

  hMuPtRec->Add(hMuPlPtRec,hMuMiPtRec);

  hMuPtRec->Draw();
  hMuPlPtRec->Draw("same");
  hMuMiPtRec->Draw("same");

  TFile *outf = NULL;
  if (saveFile) {
    outf = new TFile(fname,"RECREATE");
    hMuPtRec->Write();
    outf->Close();
  }

  return;
}

void zeeAnaDiff(TH1D* hsee, TCut addCut,TH1D* hseeMc) {
  
  TFile *f1=new TFile(fname1.Data());
  TTree *zp = (TTree*)f1->Get("tz");
  zp->AddFriend("yEvt=yongsunHiEvt"       ,fname1.Data());
  zp->AddFriend("ySkim=yongsunSkimTree"   ,fname1.Data());
  zp->AddFriend("yHlt=yongsunHltTree"     ,fname1.Data());
  zp->AddFriend("tgj",                     fname1.Data());
  
  TH1D* hmass = new TH1D("hmass",";inv. mass (GeV);event fraction",34,30,200);
  hmass->Sumw2();

  zp->Draw("invm>>hmass",addCut && "isEle1==1 && isEle2==1");
  TH1D* htemp = (TH1D*)hsee->Clone("htemp");   htemp->Reset();
  zp->Draw("see1>>htemp",addCut && "isEle1==1 && isEle2==1" && "invm>80 && invm<100" );
  hsee->Reset();
  hsee->Add(htemp);
  handsomeTH1(hmass,2);
  scaleInt(hmass);
  
  TH1D* htempMC;
  TFile *f1mc;
  TTree *zpmc ;
  f1mc=new TFile(fname1MC.Data());
  zpmc = (TTree*)f1mc->Get("tz");
  zpmc->AddFriend("yEvt=yongsunHiEvt"       ,fname1MC.Data());
  zpmc->AddFriend("ySkim=yongsunSkimTree"   ,fname1MC.Data());
  zpmc->AddFriend("yHlt=yongsunHltTree"     ,fname1MC.Data());
  zpmc->AddFriend("yPho=yongsunPhotonTree"     ,fname1MC.Data());
  zpmc->AddFriend("tgj",                     fname1MC.Data());
  hseeMc->Sumw2();
  htempMC = (TH1D*)hseeMc->Clone("htempMc");   htempMC->Reset();
  TH1D* hmassMC = new TH1D("hmassMC11",";inv. mass (GeV);event fraction",34,30,200);
  hmassMC->Sumw2();
  zpmc->Draw("invm>>hmassMC11",Form("(%s) *yPho.ncoll",addCut.GetTitle()));
  zpmc->Draw("see1>>htempMc",Form("(%s) *yPho.ncoll", (addCut && "invm>80 && invm<100").GetTitle()) );
  hseeMc->Reset();
  hseeMc->Add(htempMC);
  handsomeTH1(hmassMC,1);
  scaleInt(hmassMC);
  
  
  hmassMC->SetAxisRange(0,0.6,"Y");
  hmassMC->DrawCopy();
  hmass->DrawCopy("same");                                                                                             
  TLegend* leg1 =  new TLegend(0.5580645,0.5381356,0.8923387,0.7648305,NULL,"brNDC");
  easyLeg(leg1,"Z->ee mass");
  leg1->AddEntry(hmass,"DATA","pl");
  leg1->AddEntry(hmassMC,"MC","pl");
  leg1->Draw();
  
}

void getseePhoMC(TH1D* hsee, TCut addCut){

  TFile *f1=new TFile(fPho.Data());
  TTree *photon1 = (TTree*)f1->Get("yongsunPhotonTree");
  photon1->AddFriend("yEvt=yongsunHiEvt"       ,fPho.Data());
  photon1->AddFriend("ySkim=yongsunSkimTree"   ,fPho.Data());
  photon1->AddFriend("yHlt=yongsunHltTree"     ,fPho.Data());
  photon1->AddFriend("tgj",                     fPho.Data());
 
  TH1D* htemp = (TH1D*)hsee->Clone("htemp");  htemp->Reset();
  photon1->Draw("sigmaIetaIeta>>htemp",Form("(%s) *ncoll",addCut.GetTitle()));
  hsee->Reset();
  hsee->Add(htemp);
}

void makePlot(vector<TTree*> sigTree,vector<double> sigWeight,
	      vector<int> ptHatLo, vector<int> ptHatHi,
              std::string var,TCut cut,TH1F* h,bool norm)
{
   TH1F *hRes = (TH1F*)h->Clone();
   hRes->SetName("hRes");
   hRes->Sumw2();
  
   char tmp[300];
   for (unsigned int i=0; i<sigTree.size(); i++)
     {
       // first determine the pthat cut
       
       sprintf(tmp, "ptHat >= %d && ptHat <= %d",ptHatLo[i], ptHatHi[i]);
       TCut ptHatCut = tmp;
       TCut allCut = cut + ptHatCut;
       cout << "Current cut = " << allCut.GetTitle() << endl;
       TH1F *htmp = (TH1F*)h->Clone();
       htmp->SetName("htmp");
       sigTree[i]->Draw(Form("%s>>htmp",var.data()),allCut);
       htmp->Sumw2();
       htmp->Scale(sigWeight[i]);
       cout << "scale = " << sigWeight[i] << endl;
       cout << "After scaling htmp -> entries() " << htmp->GetEntries() << endl;
       cout << "After scaling htmp -> Integral() " << htmp->Integral() << endl;
       cout << "After scaling htmp -> GetMean()  " << htmp->GetMean() << endl;
       cout << "After scaling htmp -> GetRMS()  " << htmp->GetRMS() << endl;
       hRes->Add(htmp);
       delete htmp;
   }
   h->Sumw2();
   h->Add(hRes);
   if(norm)h->Scale(1.0/(double)h->Integral(0,1000));
   cout << "After scaling h-> entries() " << h->GetEntries() << endl;
   cout << "After scaling h-> Integral() " << h->Integral() << endl;
   cout << "After scaling h -> GetMean()  " << h->GetMean() << endl;
   cout << "After scaling h -> GetRMS()  " << h->GetRMS() << endl;
   
   delete hRes;
}

// ---------------------------------------------------------------------------------------------
TString weightCut (TCut cut, bool needMCweight, int treeVersion = 2, bool pTweighted = false)
{
  TString sWMC (treeVersion <= 1 ? "xsec" : "w");
  TString sPT (treeVersion <= 2 ? "jtpt" : "jtrawpt");
  TString title( stripStart( stripEnd( cut.GetTitle() ) ) );
  
  TString title_in_parenthesis( title );
  if( ! allInBrackets( title ) ) title_in_parenthesis = "("+title+")";
  if( needMCweight || pTweighted ) { // need some sort of weight
    TString sW;
    if( needMCweight && !pTweighted ) sW = sWMC.Data();
    if( pTweighted && !needMCweight ) sW = sPT.Data();
    if( needMCweight && pTweighted )  sW = Form("(%s*%s)", sWMC.Data(), sPT.Data());
    if( title.IsWhitespace() ) {
      return sW;
    } else {
      return title_in_parenthesis + "*" + sW;
    }
  } else { // no weights needed. 
    // But this is a cut, not a weigh, so add parenthesis and a test to prevent weighting
    if( title.IsWhitespace() ) return "";
    return title_in_parenthesis+">0";
  }
}

//.........这里部分代码省略.........
  tnpcut += os;
  tnpcut += tag_pt30;

  tnpcut += met30;
  tnpcut += nbl0;

  TCut  lepid;
  TCut  lepiso;
  char* lepchar = "";

  if( leptype == 0 ){
    tnpcut += tag_eta14;
    tnpcut += probe_eta14;
    tnpcut += eltnp;
    tnpcut += eltnptrig;
    lepid   = TCut(elid);
    lepiso  = TCut(eliso);
    lepchar = "el";
  }
  else if( leptype == 1 ){
    tnpcut += tag_eta21;
    tnpcut += probe_eta21;
    //tnpcut += mutrk;
    //tnpcut += mud0;
    //tnpcut += mudz;
    tnpcut += mutnp;
    tnpcut += mutnptrig;
    lepid   = TCut(muid);
    lepiso  = TCut(muiso);
    lepchar = "mu";
  }
  
  //tnpcut += njets2;
  cout << "Selection  : " << tnpcut.GetTitle()          << endl;
  cout << "Ndata      : " << chdata->GetEntries(tnpcut) << endl;
  cout << "NMC        : " << chmc->GetEntries(tnpcut)   << endl;
  cout << "ID cut     : " << lepid.GetTitle()           << endl;
  cout << "iso cut    : " << lepiso.GetTitle()          << endl;

  chmc->  Draw("abs(probe->eta()):probe->pt()>>hmcid_deno", 	tnpcut+lepiso,	      	"goff");
  chmc->  Draw("abs(probe->eta()):probe->pt()>>hmcid_num", 	tnpcut+lepiso+lepid,	"goff");
  chmc->  Draw("abs(probe->eta()):probe->pt()>>hmciso_deno", 	tnpcut+lepid,	      	"goff");
  chmc->  Draw("abs(probe->eta()):probe->pt()>>hmciso_num", 	tnpcut+lepid+lepiso,	"goff");
  chdata->Draw("abs(probe->eta()):probe->pt()>>hdataid_deno", 	tnpcut+lepiso,		"goff");
  chdata->Draw("abs(probe->eta()):probe->pt()>>hdataid_num", 	tnpcut+lepiso+lepid,	"goff");
  chdata->Draw("abs(probe->eta()):probe->pt()>>hdataiso_deno", 	tnpcut+lepid,	       	"goff");
  chdata->Draw("abs(probe->eta()):probe->pt()>>hdataiso_num", 	tnpcut+lepid+lepiso,	"goff");

  // get efficiencies 
  hmcid->Divide(hmcid_num,hmcid_deno,1,1,"B");
  hmciso->Divide(hmciso_num,hmciso_deno,1,1,"B");
  hdataid->Divide(hdataid_num,hdataid_deno,1,1,"B");
  hdataiso->Divide(hdataiso_num,hdataiso_deno,1,1,"B");

  // hmcid->Divide(hmcid_num,hmcid_deno,1,1);
  // hmciso->Divide(hmciso_num,hmciso_deno,1,1);
  // hdataid->Divide(hdataid_num,hdataid_deno,1,1);
  // hdataiso->Divide(hdataiso_num,hdataiso_deno,1,1);
	
  // get scale factors
  hsfid->Divide(hdataid, hmcid, 1, 1);
  hsfiso->Divide(hdataiso, hmciso, 1, 1);

  // Draw histograms	
  //hmcid->Draw("text");

void Muti_Plot_JER_ce30to100(std::string mu_title, const char **lineKind_Arr, TCut selectionC , std::string Var ,  const int ngr, const int nVarBins, double * VarBins_array , double **yArr2, double **yErrArr2)
{ 

	char filename[]="JetESR_status.txt";
	fstream myfile;
	myfile.open (filename,ios::out|ios::app);

	// Double_t Xptbin_chris_mean_ce0to30[20] = {7.999, 9.680157, 11.71464, 14.17672, 17.15625, 20.762, 25.12557, 30.40623, 36.79673, 44.53033, 53.88931, 65.21529, 78.92163, 95.50865, 115.5818, 139.8737, 169.2711, 204.8469, 247.8998, 300.001};
  Double_t Xptbin_chris_mean_ce30to100[21] = {-2,-1.8,-1.6,-1.4,-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8,2};

	Double_t Xpt_chris_mean_ce30to100[20];
	Double_t XErrpt_chris_mean_ce30to100[20];
	for (int ibin =0; ibin<20;ibin++){
		Xpt_chris_mean_ce30to100[ibin]= (Xptbin_chris_mean_ce30to100[ibin+1] +Xptbin_chris_mean_ce30to100[ibin])/2;
		XErrpt_chris_mean_ce30to100[ibin]= (Xptbin_chris_mean_ce30to100[ibin+1] -Xptbin_chris_mean_ce30to100[ibin])/2;
	}
//	double ypt_chris_mean_ce0to30[19]={2.614117,2.282987,2.04587,1.74664,1.548742,1.396857,1.275476,1.194122,1.140162,1.120749,1.075858,1.05422,1.009865,0.9785294,0.945,0.99,0.9000001,1.163599,1.15};
//	double yerrpt_chris_mean_ce0to30[19]={0.03932904,0.02605798,0.0244768,0.02502056,0.02119035,0.02074938,0.01284258,0.01424918,0.01133906,0.01246782,0.01578119,0.0227899,0.03189043,0.02962023,0.0438062,0.05302978,0.01194279,0.009490237,0.01854081};
  double ypt_chris_mean_ce30to100[20]={0.2311591,0.2477474,0.2655955,0.2576075,0.2220231,0.206768,0.2064983,0.2020811,0.19917,0.2115871,0.1952014,0.1965805,0.2158946,0.2003286,0.213619,0.213695,0.2466632,0.2597152,0.2553281,0.2344014};
  double yerrpt_chris_mean_ce30to100[20]={0.009894431,0.01000488,0.01052403,0.009913006,0.008527985,0.007570562,0.007309829,0.007320613,0.00699798,0.00776784,0.007146697,0.00713059,0.007530691,0.007364154,0.007896077,0.007855752,0.009496669,0.01034782,0.01046585,0.009252961};

	TGraphErrors *gr_chris_JERpt_ce30to100= new TGraphErrors(20,Xpt_chris_mean_ce30to100,ypt_chris_mean_ce30to100,XErrpt_chris_mean_ce30to100,yerrpt_chris_mean_ce30to100);

	gr_chris_JERpt_ce30to100->SetFillStyle(0);
	gr_chris_JERpt_ce30to100->SetMarkerStyle(1);
	gr_chris_JERpt_ce30to100->SetLineColor(ngr+1);
	gr_chris_JERpt_ce30to100->SetMarkerColor(ngr+1);


	double VarBinsCenter[nVarBins];
	double VarBinsErr[nVarBins];
	for (int ibin=0; ibin<nVarBins; ibin++)
	{
		VarBinsCenter[ibin]= (VarBins_array[ibin] + VarBins_array[ibin+1]) /2;
		//		VarBinsErr[ibin]=0;
		VarBinsErr[ibin]=(VarBins_array[ibin+1]- VarBins_array[ibin])/2;
	}

	cout<<"lineKind_Arr 1 = "<<lineKind_Arr[0]<<" ,lineKind_Arr2 = "<<lineKind_Arr[1]<<endl;
	cout<<"mu_title = "<<mu_title<<" selectionC = "<< selectionC<<endl;

	Can_result[counter1] = new TCanvas(Form("Can_result_%i",counter1));
	mutiGR[counter1] = new TMultiGraph("name",Form("%s",mu_title.c_str()));
	TGraphErrors *gr_mean[ngr];
	legend[counter1] = new TLegend(0.62,0.7,0.82,0.85);

	for(int igr=0; igr<ngr; igr++)
	{
		gr_mean[igr] = new TGraphErrors(nVarBins, VarBinsCenter, yArr2[igr], VarBinsErr , yErrArr2[igr]);
		//    char *name_temp = {Form("gr_mean_%i", (int)igr)};
		//		gr_mean[igr]->SetName(Form("gr_mean_%i", (int)igr) );
		//		gr_mean[igr]->SetName(name_temp);
		gr_mean[igr]->SetTitle(Form("some title_%i", igr));
		gr_mean[igr]->SetFillStyle(0);
		gr_mean[igr]->SetMarkerStyle(1);
		gr_mean[igr]->SetLineColor(igr+1);	
		gr_mean[igr]->SetMarkerColor(igr+1);
		mutiGR[counter1]->Add(gr_mean[igr]);
		legend[counter1]->AddEntry( gr_mean[igr] ,lineKind_Arr[igr],"l");  // not work for unknow reason

	}

	mutiGR[counter1]->Add(gr_chris_JERpt_ce30to100);
	legend[counter1]->AddEntry(gr_chris_JERpt_ce30to100,"from chris","l");

	myfile<<"hello, Var = "<<Var<<"Var.compare(Centrality) = "<<Var.compare("Centrality")<<endl;

	mutiGR[counter1]->Draw("AP");
	mutiGR[counter1]->GetXaxis()->SetTitle(Form("%s",Var.c_str())); // must after draw to create a vitual histogram like object to set title.
	mutiGR[counter1]->GetYaxis()->SetTitle("#sigma_{Reco./Gen.}");
	mutiGR[counter1]->Draw("AP"); // draw again to let all the settting on it.
	//	Can_result[counter1]->BuildLegend(0.7,0.7,0.9,0.9);
	legend[counter1]->SetBorderSize(0);
	legend[counter1]->Draw();

	TLatex *tex = new TLatex(0.35,0.8,"#bf{#color[2]{|#eta_{jet}|<2.0} & refpt>30}");
	tex->SetNDC();
	tex->SetTextFont(43);
	tex->SetTextSize(26);
	tex->SetLineWidth(2);
	tex->Draw();
  tex = new TLatex(0.17,0.25,"Centrality:30-100%");
	tex->SetNDC();
  tex->SetTextFont(43);
  tex->SetTextSize(26);
  tex->SetLineWidth(2);
  tex->Draw();


	Can_result[counter1]->SaveAs(Form("./Plots/result_%s_%s_%s.pdf",mu_title.c_str(),selectionC.GetTitle() , Var.c_str())); 
	Can_result[counter1]->SaveAs(Form("./Plots/result_%s_%s_%s.C",mu_title.c_str(),selectionC.GetTitle() , Var.c_str()));
	counter1++;
	cout<<"finish mutigraph"<<endl;
}

//TH1F* doDYestimate( TChain *ch , TCut sel ){
void doDYestimate( TChain *ch , TCut sel , TH1F* hyield ){

  //-------------------
  // set parameters
  //-------------------

  bool  verbose = true;
  float R       = 0.13;
  float R_err   = 0.07;
  float k       = 1.11;

  if( verbose ){
    cout << "R : " << R << " +/- " << R_err << endl;
    cout << "k : " << k << endl;
  }

  //-------------------
  // invert Z-veto
  //-------------------

  TString selstring(sel.GetTitle());
  selstring.ReplaceAll("passz == 0","dilmass>76&&dilmass<106");
  TCut newsel = TCut(selstring);

  if( verbose ){
    cout << "Pre  : " << sel.GetTitle() << endl;
    cout << "Post : " << newsel.GetTitle() << endl;
  }

  //-------------------
  // get data yields
  //-------------------

  float ndata_ee_in = ch->GetEntries(newsel+"leptype==0");
  float ndata_mm_in = ch->GetEntries(newsel+"leptype==1");
  float ndata_em_in = ch->GetEntries(newsel+"leptype==2");

  if( verbose ){
    cout << "nee     " << ndata_ee_in << endl;
    cout << "nmm     " << ndata_mm_in << endl;
    cout << "nem     " << ndata_em_in << endl;
  }

  //-------------------
  // do DY estimate
  //-------------------

  float neepred     = R     * ( ndata_ee_in - (0.5/k)  * ndata_em_in );
  float nmmpred     = R     * ( ndata_mm_in - k/2.     * ndata_em_in );

  float neeprederr  = 0.;
  float nmmprederr  = 0.;

  neeprederr  += pow( R_err * ( ndata_ee_in - (0.5/k) * ndata_em_in ) , 2 );
  nmmprederr  += pow( R_err * ( ndata_mm_in - k/2.    * ndata_em_in ) , 2 );

  neeprederr  += pow( R * sqrt( ndata_ee_in + pow(0.5/k,2)  * ndata_em_in ) , 2 );
  nmmprederr  += pow( R * sqrt( ndata_mm_in + pow(k/2.,2)   * ndata_em_in ) , 2 );

  neeprederr = sqrt(neeprederr);
  nmmprederr = sqrt(nmmprederr);

  float ntotpred    = neepred + nmmpred;
  float ntotprederr = sqrt( pow(neeprederr,2) + pow(nmmprederr,2) );

  if( verbose ){
    cout << "nee  pred " << neepred  << " +/- " << neeprederr  << endl;
    cout << "nmm  pred " << nmmpred  << " +/- " << nmmprederr  << endl;
    cout << "ntot pred " << ntotpred << " +/- " << ntotprederr << endl;
  }

  //-------------------
  // set hist contents
  //-------------------

  hyield->SetBinContent( 1 , neepred );
  hyield->SetBinContent( 2 , nmmpred );
  hyield->SetBinContent( 3 , 0 );
  hyield->SetBinContent( 4 , 0 );

  hyield->SetBinError( 1 , neeprederr );
  hyield->SetBinError( 2 , nmmprederr );
  hyield->SetBinError( 3 , 0 );
  hyield->SetBinError( 4 , 0 );

}

void printYields( vector<TChain*> chmc , vector<char*> labels , TChain* chdata , TCut sel , TCut weight , bool latex ){

  initSymbols( latex );

  TCanvas *ctemp = new TCanvas();

  printLine(latex);
  printHeader();
  printLine(latex);

  TH1F* hyield = new TH1F("hyield","yield",4,0,4);
  TH1F* hmctot = new TH1F("hmctot","hmctot",4,0,4);
  hyield->Sumw2();
  hmctot->Sumw2();

  //----------------------
  // print SM MC samples
  //----------------------

  for(unsigned int imc = 0 ; imc < chmc.size() ; imc++){

    bool correlatedError = false;

    if( TString(labels[imc]).Contains("T2") ) continue;

    // data-driven DY estimate
    if( strcmp(labels[imc],"DYdata")   == 0 ){
      //hyield = doDYestimate( chmc[imc] , sel );
      doDYestimate( chmc[imc] , sel , hyield );
    }

    // fake estimate
    else if( strcmp(labels[imc],"single fakes")   == 0 || strcmp(labels[imc],"double fakes") == 0){

      //correlatedError = true;

      TString weightstring(weight.GetTitle());
      weightstring.ReplaceAll("ndavtxweight","1");
      TCut newweight = TCut(weightstring);

      chmc[imc]->Draw("leptype>>hyield",sel*newweight);

      // SF --> SF - 2 X DF
      if( strcmp(labels[imc],"single fakes")   == 0 ){
	
	TH1F *hyielddf = new TH1F("hyielddf","hyielddf",4,0,4);

	chmc[imc+1]->Draw("leptype>>hyielddf",sel*newweight);
	hyield->Add(hyielddf,-2);
      }      

      hyield->SetBinError(1,0.5*hyield->GetBinContent(1));
      hyield->SetBinError(2,0.5*hyield->GetBinContent(2));
      
    }

    //vanilla MC
    else{
      chmc[imc]->Draw("leptype>>hyield",sel*weight);

      //do efficiency correction
      //hyield->SetBinContent  ( 2 , hyield->GetBinContent(2) * 0.90);
      //hyield->SetBinContent  ( 3 , hyield->GetBinContent(3) * 0.95);
      //hyield->SetBinError    ( 2 , hyield->GetBinError(2)   * 0.90);
      //hyield->SetBinError    ( 3 , hyield->GetBinError(3)   * 0.95);

    }

    if( imc == 0 ) hmctot = (TH1F*) hyield->Clone();
    else           hmctot->Add(hyield);
    
    print( hyield , labels[imc] , correlatedError );

    //hyield->Reset();
  }

  printLine(latex);

  //-------------------------------
  // print sum of SM MC samples
  //-------------------------------

  print( hmctot , "total SM MC" );

  printLine(latex);
 
  chdata->Draw("leptype>>hyield",sel);

  print( hyield , "data" );
    
  printLine(latex);

  //----------------------
  // print T2tt MC samples
  //----------------------

  for(unsigned int imc = 0 ; imc < chmc.size() ; imc++){

    if( !TString(labels[imc]).Contains("T2") ) continue;

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

void gammaTrkHistProducer(sampleType collision = kPPDATA, float photonPtThr=20, float photonPtThrUp=9999, int icent =7){
  TH1::SetDefaultSumw2();
  
  TString stringSampleType = getSampleName(collision); "";
  
  TDatime* date = new TDatime();
  TString  outName=  Form("photonTrackCorr_%s_output_photonPtThr%d_to_%d_%d.root",stringSampleType.Data(),(int)photonPtThr, (int)photonPtThrUp,  date->GetDate());
  delete date;
  
 
  int lowerCent(0),  upperCent(0); 
  TCut centCut  = "";
  if ( (collision ==kHIDATA) || (collision ==kHIMC) ) {
     lowerCent = ((icent/100)%100) *2 ;
     upperCent =  (icent%100)*2 -1 ; 
     centCut = Form("cBin >= %d && cBin<= %d",lowerCent,upperCent);
  }
  else if (  (collision ==kPPDATA) || (collision==kPPMC)  ){  // if it's pp 
    centCut = "(1==1)";
    //    icent = 7;   // for pp, centrality is set as the smearing 
  }
  else { // pPb
    centCut = Form( "hf4Sum > %f && hf4Sum <= %f", (float)centBinPa[icent-1], (float)centBinPa[icent]);
  }
  
  cout <<" centrality : " << centCut.GetTitle() << endl;
      
  ///////// Photon cut //////////////////////////////////////////////////////////////////////////////
  
  cout <<" photon pt >" << photonPtThr << " GeV" << endl;
  TCut ptPhoCut  = Form("photonEt>%.1f && photonEt<%.1f", (float)photonPtThr, (float)photonPtThrUp  );
  TCut caloIso;
  
  if ( (collision==kPPMC) || (collision==kPPDATA) ) 
    caloIso = "(ecalIso < 4.2  &&  hcalIso < 2.2  &&  trackIso < 2) && hovere<0.1";
  else if ( (collision==kHIMC) || (collision==kHIDATA) )
    caloIso = "(sumIso<5) && hovere<0.1";
  else {
    caloIso = "ecalIso < 4.2  &&  hcalIso < 2.2  &&  trackIso < 2 && hovere<0.1";
  }
  
  TCut sbIso   = "(sumIso>5) && (sumIso<20) && hovere<0.1";
  //  if ( (collision==kPPMC) || (collision==kPPDATA) || (collision==kPAMC) || (collision==kPADATA)  )
  //  sbIso   = "ecalIso < 4.2  &&  hcalIso < 2.2 && trackIso > 2 && trackIso < 5 && hovere<0.1";

  TCut basicPhoCut = centCut && ptPhoCut && caloIso ;
  TCut sbPhoCut    = centCut && ptPhoCut && sbIso   ;
  TCut evtSeltCut = basicPhoCut;
  TCut sbSeltCut  = sbPhoCut;

  TCut phoCandCut   = "sigmaIetaIeta<0.010";
  TCut phoDecayCut  = "(sigmaIetaIeta>0.011) && (sigmaIetaIeta<0.017)";
  if  ( ( collision == kHIMC ) || (collision == kPPMC) || (collision == kPAMC))  
    phoCandCut = phoCandCut && "genIso<5 && abs(genMomId)<=22";
  
  
  TString fname = "";
  if ( collision == kHIDATA)      fname = fnameHIDATA; //
  else if ( collision == kPADATA) fname = fnamePADATA;
  else if ( collision == kPPDATA) {
    if  ( icent == 7 ) fname = fnamePPDATA;
    else if ( icent == 10010 ) fname = fnamePPDATA0010;
    else if ( icent == 11030 ) fname = fnamePPDATA1030;
    else if ( icent == 13050 ) fname = fnamePPDATA3050;
    else if ( icent == 15099 ) fname = fnamePPDATA5099;
    else if ( icent == 10030 ) fname = fnamePPDATA0030;
    else if ( icent == 13099 ) fname = fnamePPDATA30100;
  }  
  else fname = "";
  
  multiTreeUtil* tgj = new multiTreeUtil();
  tgj->addFile(fname,  "tgj",  evtSeltCut,  1);
  tgj->AddFriend("yTrk");
 
  float purity(0);
  
  TString canvasName = Form("gifs/purity_%s_output_icent%d_photonPtThr%d-%d", stringSampleType.Data(),  (int)icent, (int)photonPtThr, (int)photonPtThrUp);
  
  if ( collision == kPPDATA) {  
    purity = 0.85;
  }
  else {
    fitResult fitr = getPurity(fname, collision, evtSeltCut, sbSeltCut, canvasName, photonPtThr, photonPtThrUp);
    purity = fitr.purity010;
  }
  GjSpectra* gSpec = new GjSpectra();
  gSpec->init(Form("icent%d",(int)icent) );
  tgj->Draw2(gSpec->hPtPhoCand,  "photonEt", phoCandCut, "");
  tgj->Draw2(gSpec->hPtPhoDecay, "photonEt", phoDecayCut, "");
  
  // Obtain background subtracted spectra
  
  float candInt = gSpec->hPtPhoCand->Integral();
  float decayInt = gSpec->hPtPhoDecay->Integral();
  gSpec->hPtPhoSig->Reset();
  gSpec->hPtPhoSig->Add(gSpec->hPtPhoCand);
  gSpec->hPtPhoSig->Add(gSpec->hPtPhoDecay, -(1. - purity) * candInt / decayInt);
  gSpec->hPtPhoSig->Scale(1./purity ) ;
  
  TFile outf = TFile(Form("ffFiles/%s",outName.Data()),"update");
//.........这里部分代码省略.........

void fit_combine(TChain *tc, std::string Var, const double * VarBins_array , const int nVarBins , TCut Gcut,  double *mean, double *meanErr, double *sigma, double *sigmaErr ){

	std::string jetTitle=tc->GetTitle();

	cout<<"tc title = "<<tc->GetTitle()<<" ,Var = "<<Var<<" ,nVarBins = "<<nVarBins<<", Gcut = "<<Gcut<<endl;
	Can_Temp[counter] = new TCanvas(Form("Can_Temp_%i",counter));
	if(nVarBins<=9){Can_Temp[counter]->Divide(3,3);}
	else {Can_Temp[counter]->Divide(3,4);}


	TH1D *h_Ratio[nVarBins];
	TF1 *f_Ratio[nVarBins];
	double nRatiobins= 100;
	double minRatio=0;
	double maxRatio=3;


	for(int ibin =0; ibin<nVarBins; ibin++)
	{
		Can_Temp[counter]->cd(ibin+1);
		double VarMin = VarBins_array[ibin];
		double VarMax = VarBins_array[ibin+1];

		h_Ratio[ibin] = new TH1D( Form("h_Ratio[%d]",ibin),Form("%s_%s_%i_%i",jetTitle.c_str(),Var.c_str(),(int)VarMin,(int)VarMax) ,nRatiobins,minRatio,maxRatio);  // h_Ration_%d ??
		h_Ratio[ibin]->Sumw2();

		TCut Varcut = Form("%s>%f && %s<%f", Var.c_str(), VarMin, Var.c_str(), VarMax);
		TCut allcut = Gcut && Varcut;
		cout<<"Varcut = "<<Varcut<<" ,allcut = "<<allcut<<endl;

		tc->Draw( Form("jtpt/refpt>>h_Ratio[%d]",ibin ) , allcut*"weight" );

		h_Ratio[ibin]->GetXaxis()->SetTitle(Form("%s, jtpt/refpt",Gcut.GetTitle()));

		f_Ratio[ibin] = new TF1(Form("f_Ratio_%d",ibin), "gaus",0.5,1.5 );
		f_Ratio[ibin]->SetParameter(0,h_Ratio[ibin]->GetEntries()); // get number for normalization
		if(ibin==0){
		f_Ratio[ibin]->SetParameter(1,h_Ratio[ibin]->GetMean()) ;   // for mean
		f_Ratio[ibin]->SetParameter(2,h_Ratio[ibin]->GetStdDev()) ; // for sigma   
		}
		if(ibin>=1){
    f_Ratio[ibin]->SetParameter(1,mean[ibin-1]) ;   // for mean
    f_Ratio[ibin]->SetParameter(2,sigma[ibin-1]) ; // for sigma   
		} 
		h_Ratio[ibin]->Fit(Form("f_Ratio_%d",ibin),"MR");
    h_Ratio[ibin]->Fit(Form("f_Ratio_%d",ibin),"MRL");
    h_Ratio[ibin]->Fit(Form("f_Ratio_%d",ibin),"MR");
		h_Ratio[ibin]->Fit(Form("f_Ratio_%d",ibin),"MR");


		mean[ibin] = f_Ratio[ibin]->GetParameter(1);
		meanErr[ibin] = f_Ratio[ibin]->GetParError(1);
		sigma[ibin] = f_Ratio[ibin]->GetParameter(2);
		sigmaErr[ibin] = f_Ratio[ibin]->GetParError(2);

	} // end   for(int ibin =0; ibin<nVarBins; ibin++)

	Can_Temp[counter]->SaveAs(Form("./fitPlots/fit_%s_%s_%s.pdf",Var.c_str(),jetTitle.c_str(),Gcut.GetTitle() ));

	counter++; 
}

void printHisto( TCanvas *can , TChain *data , TChain *mc , TCut num , TCut denom , char* var , int nbins , float xmin , float xmax , char* xtitle , char* ytitle){

  can->cd();

  TPad *plotpad = new TPad("plotpad","plotpad",0.0,0.0,1.0,0.8);
  plotpad->Draw();
  plotpad->cd();

  float ptbin[] = {30., 40., 60. , 80. , 100. , 120. , 150. , 200 , 300 };
  int   nptbin  = 8;

  //TH1F* hpass   = new TH1F(Form("hpass_%i",iplot),Form("hpass_%i",iplot),nptbin,ptbin);
  //TH1F* hall    = new TH1F(Form("hall_%i" ,iplot),Form("hall_%i" ,iplot),nptbin,ptbin);

  // TH1F* hpass_data = new TH1F(Form("hpass_data_%i",iplot),Form("hpass_data_%i",iplot),nbins,xmin,xmax);
  // TH1F* hall_data  = new TH1F(Form("hall_data_%i" ,iplot),Form("hall_data_%i" ,iplot),nbins,xmin,xmax);
  // TH1F* hpass_mc   = new TH1F(Form("hpass_mc_%i"  ,iplot),Form("hpass_mc_%i"  ,iplot),nbins,xmin,xmax);
  // TH1F* hall_mc    = new TH1F(Form("hall_mc_%i"   ,iplot),Form("hall_mc_%i"   ,iplot),nbins,xmin,xmax);

  TH1F* hpass_data   = new TH1F(Form("hpass_data_%i"  ,iplot),Form("hpass_data_%i"  ,iplot),nptbin,ptbin);
  TH1F* hall_data    = new TH1F(Form("hall_data_%i"   ,iplot),Form("hall_data_%i"   ,iplot),nptbin,ptbin);
  TH1F* hratio_data  = new TH1F(Form("hratio_data_%i" ,iplot),Form("hratio_data_%i" ,iplot),nptbin,ptbin);
  TH1F* hpass_mc     = new TH1F(Form("hpass_mc_%i"    ,iplot),Form("hpass_mc_%i"    ,iplot),nptbin,ptbin);
  TH1F* hall_mc      = new TH1F(Form("hall_mc_%i"     ,iplot),Form("hall_mc_%i"     ,iplot),nptbin,ptbin);
  TH1F* hratio_mc    = new TH1F(Form("hratio_mc_%i"   ,iplot),Form("hratio_mc_%i"   ,iplot),nptbin,ptbin);
  TH1F* hsf          = new TH1F(Form("hsf_%i"         ,iplot),Form("hsf_%i"         ,iplot),nptbin,ptbin);

  hpass_data->Sumw2();
  hall_data->Sumw2();
  hratio_data->Sumw2();
  hpass_mc->Sumw2();
  hall_mc->Sumw2();
  hratio_mc->Sumw2();
  hsf->Sumw2();

  //TCanvas *can = new TCanvas(Form("can_%i",iplot),Form("can_%i",iplot),600,600);
  //can->cd();
  
  data->Draw(Form("min(%s,%f)>>hpass_data_%i"  , var,xmax-0.0001,iplot),denom+num);
  data->Draw(Form("min(%s,%f)>>hall_data_%i"   , var,xmax-0.0001,iplot),denom);
  mc->Draw  (Form("min(%s,%f)>>hpass_mc_%i"    , var,xmax-0.0001,iplot),denom+num);
  mc->Draw  (Form("min(%s,%f)>>hall_mc_%i"     , var,xmax-0.0001,iplot),denom);

  TGraphAsymmErrors *grdata = new TGraphAsymmErrors();
  grdata->BayesDivide(hpass_data,hall_data);

  TGraphAsymmErrors *grmc = new TGraphAsymmErrors();
  grmc->BayesDivide(hpass_mc,hall_mc);

  // cout << "data all  " << hall_data->GetBinContent(8) << endl;
  // cout << "data pass " << hpass_data->GetBinContent(8) << endl;
  // cout << "data eff  " << hpass_data->GetBinContent(8) / hall_data->GetBinContent(8) << endl;

  // Double_t x;
  // Double_t y;
  // grdata->GetPoint(7,x,y);
  // cout << "data eff2 " << y << endl;

  gPad->SetGridx();
  gPad->SetGridy();

  grdata->SetMarkerColor(2);
  grdata->SetLineColor(2);
  grmc->SetMarkerColor(4);
  grmc->SetLineColor(4);
  grmc->SetMarkerStyle(25);

  grdata->GetXaxis()->SetRangeUser(30,300);
  grdata->GetYaxis()->SetRangeUser(0.6,1.0);
  grdata->GetXaxis()->SetTitle(xtitle);
  grdata->GetYaxis()->SetTitle(ytitle);
  grdata->Draw("AP");
  grmc->Draw("sameP");

  TLegend *leg = new TLegend(0.5,0.2,0.7,0.4);
  leg->AddEntry(grdata ,"data","lp");
  leg->AddEntry(grmc   ,"mc"  ,"lp");
  leg->SetBorderSize(1);
  leg->SetFillColor(0);
  //leg->Draw();

  TLatex *t = new TLatex();
  t->SetNDC();

  if( TString(denom.GetTitle()).Contains("njets==0") ) t->DrawLatex(0.2,0.2,"n_{jets}=0");
  if( TString(denom.GetTitle()).Contains("njets==1") ) t->DrawLatex(0.2,0.2,"n_{jets}=1");
  if( TString(denom.GetTitle()).Contains("njets==2") ) t->DrawLatex(0.2,0.2,"n_{jets}=2");
  if( TString(denom.GetTitle()).Contains("njets==3") ) t->DrawLatex(0.2,0.2,"n_{jets}=3");
  if( TString(denom.GetTitle()).Contains("njets>=4") ) t->DrawLatex(0.2,0.2,"n_{jets}#geq4");

  can->cd();

  TPad *respad = new TPad("respad","respad",0.0,0.8,1.0,1.0);
  respad->Draw();
  respad->cd();
  respad->SetGridy();

  // TGraphAsymmErrors* gr_ratio = (TGraphAsymmErrors*) grdata->Clone("gr_ratio");
  // gr_ratio->Divide(grmc);
  // gr_ratio->Draw();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            }
            else {
                std::cout << Form("entries[%d] = ", iInFileArg) << entriesTmp << std::endl;
            }

            output->cd();
            for (int i=0; i<nHistos; ++i) {

                int treeIndex = 0;
                if (nHistosInput == nTrees)  treeIndex = i%nTrees;
                // std::cout << "treePath = " << treePaths.at(treeIndex).c_str() << ", ";

                int iInFileArg = 0;
                if (mode == INPUT_MODE::k_comparison) {
                    iInFileArg = i%nInputFileArguments;
                    std::cout << "iInFileArg = " << iInFileArg << ", ";
                }

                std::string formula = formulas.at(0).c_str();
                std::string selection = selections.at(0).c_str();
                std::string weight = weights.at(0).c_str();
                if (nHistosInput == nFormulas)  formula = formulas.at(i%nFormulas).c_str();
                if (nHistosInput == nSelections)  selection = selections.at(i%nSelections).c_str();
                if (nHistosInput == nWeights)  weight = weights.at(i%nWeights).c_str();

                std::string selectionSplit = "";
                if (nSelectionSplitter > 1)  selectionSplit = selectionSplitter.at(i/ (nHistos/nSelectionSplitter)).c_str();

                // std::cout << "drawing histogram i = " << i << ", ";

                TCut selectionFinal = selectionBase.c_str();
                selectionFinal = selectionFinal && selection.c_str();
                if (selectionSplit.size() > 0)  selectionFinal = selectionFinal && selectionSplit.c_str();
                Long64_t entriesSelectedTmp = trees[treeIndex][iInFileArg]->GetEntries(selectionFinal.GetTitle());
                // std::cout << "entriesSelected in file = " << entriesSelectedTmp << std::endl;
                entriesSelected[i] += entriesSelectedTmp;

                TCut weight_AND_selection = Form("(%s)*(%s)", weight.c_str(), selectionFinal.GetTitle());
                trees[treeIndex][iInFileArg]->Draw(Form("%s >>+ %s", formula.c_str(), h[i]->GetName()), weight_AND_selection.GetTitle(), "goff");
            }
            fileTmp->Close();
        }
    }
    std::cout << "TTree::Draw() ENDED" <<std::endl;
    for (int i = 0; i < nInputFileArguments; ++i) {

        if (nInputFileArguments == 1)  {
            std::cout << "entries = " << entries[0] << std::endl;
        }
        else {
            std::cout << Form("entries[%d] = ", i) << entries[i] << std::endl;
        }
    }
    std::cout << "### selected entries" << std::endl;
    for (int i = 0; i < nHistos; ++i) {

        std::cout << "TH1D i = " << i << ", ";
        int treeIndex = 0;
        if (nHistosInput == nTrees)  treeIndex = i%nTrees;
        std::cout << "treePath = " << treePaths.at(treeIndex).c_str() << ", ";

        std::cout << "entriesSelected = " << entriesSelected[i] << std::endl;
    }
    std::cout << "###" << std::endl;

    // print info about histograms

void gammaJetHistProducer_jetEnergyScaledMinus2percent(sampleType collision = kPADATA, float photonPtThr=60, float photonPtThrUp=9999, float jetPtThr=30, int icent =1) {
  TH1::SetDefaultSumw2();
  
  TString stringSampleType = getSampleName(collision); "";
  
  TDatime* date = new TDatime();
  TString  outName=  Form("photonTrackCorr_%s_output_photonPtThr%d_to_%d_jetPtThr%d_%d.root",stringSampleType.Data(),(int)photonPtThr, (int)photonPtThrUp, (int)jetPtThr,  date->GetDate());
  delete date;
  
  int lowerCent(0),  upperCent(0); 
  TCut centCut  = "";
  if ( (collision ==kHIDATA) || (collision==kHIMC) )   {
    lowerCent = centBin1[icent-1];
    upperCent = centBin1[icent]-1;
    if ( icent > 9999) {
      lowerCent = ((icent/100)%100)/2.5;
      upperCent =  (icent%100)/2.5 -1;
    }   
    centCut = Form("cBin >= %d && cBin<= %d",lowerCent,upperCent);
  }
  else if (  (collision ==kPPDATA) || (collision==kPPMC)  ){  // if it's pp 
    centCut = "(1==1)";
    //    icent = 7;   // for pp, centrality is set as the smearing 
  }
  else { // pPb
    centCut = Form( "hf4Sum > %f && hf4Sum <= %f", (float)centBinPa[icent-1], (float)centBinPa[icent]);
  }
  
  cout <<" centrality : " << centCut.GetTitle() << endl;
      
  ///////// Photon cut //////////////////////////////////////////////////////////////////////////////
  
  cout <<" photon pt >" << photonPtThr << " GeV" << endl;
  TCut ptPhoCut  = Form("photonEt>%.1f && photonEt<%.1f", (float)photonPtThr, (float)photonPtThrUp  );
  TCut caloIso;
  
  if ( (collision==kPPMC) || (collision==kPPDATA) ) 
    caloIso = "(ecalIso < 4.2  &&  hcalIso < 2.2  &&  trackIso < 2) && hovere<0.1";
  else if ( (collision==kHIMC) || (collision==kHIDATA) )
    caloIso = "(sumIso<1) && hovere<0.1";
  else {
    caloIso = "ecalIso < 4.2  &&  hcalIso < 2.2  &&  trackIso < 2 && hovere<0.1";
  }
  
  TCut sbIso   = "(sumIso>10) && (sumIso<20) && hovere<0.1";
  //  if ( (collision==kPPMC) || (collision==kPPDATA) || (collision==kPAMC) || (collision==kPADATA)  )
  //  sbIso   = "ecalIso < 4.2  &&  hcalIso < 2.2 && trackIso > 2 && trackIso < 5 && hovere<0.1";

  TCut basicPhoCut = centCut && ptPhoCut && caloIso ;
  TCut sbPhoCut    = centCut && ptPhoCut && sbIso   ;
  TCut evtSeltCut = basicPhoCut;
  TCut sbSeltCut  = sbPhoCut;

  TCut phoCandCut   = "sigmaIetaIeta<0.010";
  TCut phoDecayCut  = "(sigmaIetaIeta>0.011) && (sigmaIetaIeta<0.017)";
  if  ( ( collision == kHIMC ) || (collision == kPPMC) || (collision == kPAMC))  
    phoCandCut = phoCandCut && "genIso<5 && abs(genMomId)<=22";
  
  
  TString fname = "";
  if ( collision == kHIDATA)      fname = fnameHIDATA_Minus2percentScaled;
  else if ( collision == kPADATA) fname = fnamePADATA_Minus2percentScaled;
  else if ( collision == kPPDATA) {
    if ( icent == 7 ) fname = fnamePPDATA_Minus2percentScaled;
  }  
  else fname = "";
  
  multiTreeUtil* tgj = new multiTreeUtil();
  multiTreeUtil* tgjMC = new multiTreeUtil();
  if (  ( collision == kHIDATA)   || ( collision==kPADATA) || ( collision == kPPDATA) ) {
    tgj->addFile(fname,  "tgj",  evtSeltCut,  1);
  }
  else if ( collision == kPPMC ) {
    tgj->addFile(fnamePPMC_AllQcdPho30to50,    "tgj", evtSeltCut, wPPMC_AllQcdPho30to50 );
    tgj->addFile(fnamePPMC_AllQcdPho50to80,    "tgj", evtSeltCut, wPPMC_AllQcdPho50to80 );
    tgj->addFile(fnamePPMC_AllQcdPho80to120,   "tgj", evtSeltCut, wPPMC_AllQcdPho80to120 );
    tgj->addFile(fnamePPMC_AllQcdPho120to9999, "tgj", evtSeltCut, wPPMC_AllQcdPho120to9999 );
  }
  else if ( collision == kPAMC ) {
    tgj->addFile(fnamePAMC_AllQcdPho30to50,    "tgj", evtSeltCut, wPAMC_AllQcdPho30to50 );
    tgj->addFile(fnamePAMC_AllQcdPho50to80,    "tgj", evtSeltCut, wPAMC_AllQcdPho50to80 );
    tgj->addFile(fnamePAMC_AllQcdPho80to120,   "tgj", evtSeltCut, wPAMC_AllQcdPho80to120 );
    tgj->addFile(fnamePAMC_AllQcdPho120to9999, "tgj", evtSeltCut, wPAMC_AllQcdPho120to9999 );
  }
  else  {    // kHIMC
    tgj->addFile(fnameHIMC_AllQcdPho30to50,    "tgj", evtSeltCut, wHIMC_AllQcdPho30to50 );
    tgj->addFile(fnameHIMC_AllQcdPho50to80,    "tgj", evtSeltCut, wHIMC_AllQcdPho50to80 );
    tgj->addFile(fnameHIMC_AllQcdPho80to9999,  "tgj", evtSeltCut, wHIMC_AllQcdPho80to9999 );
  }
  tgj->AddFriend("yJet");
 
  // get purity with the current jet cut ! 
  float purity(0);
  
  TString canvasName = Form("gifs/purity_%s_output_icent%d_photonPtThr%d-%d_jetPtThr%d", stringSampleType.Data(),  (int)icent, (int)photonPtThr, (int)photonPtThrUp, (int)jetPtThr);
  
  //  if ( (collision==kPPDATA) && (photonPtThr < 50 ) ) {
  //   purity = 0.86 ;   
  //   cout << " !!!!!!!" << endl << endl << " purity is set as  0.86 for this bin because we don't have pp MC low pt sample " << endl;
  //   cout << endl << endl << endl << " !!!!!!" << endl;
//.........这里部分代码省略.........