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

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


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

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

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

TH1D* GetFinalCorrFunc1D(int itrg, int jass)
{
                    TH1D* hbackground_phi = (TH1D*)GetRawBackground1D(itrg,jass);
		    TH1D* hcorrphi = (TH1D*)GetJetCorrFunc1D_ZYAM(itrg,jass);
                    hcorrphi->Multiply(hbackground_phi);
                    float max = hcorrphi->GetBinContent(hcorrphi->GetMaximumBin());
		    hcorrphi->SetAxisRange(ymin,max*1.5,"Y");
                    hcorrphi->SetYTitle("#frac{1}{N^{trg}}#frac{d^{2}N}{d#Delta#phi}");
                    return hcorrphi;
}

void ReFit(int lo, int hi) {
	TF1* floo = (TF1*) gROOT->FindObject("floo");
	if (floo == 0) return;
	TH1D* hCe = (TH1D*) gROOT->FindObject("hCe");
	if (hCe == 0) return;
	hCe->GetXaxis()->SetRangeUser(lo,hi);
	double max = hCe->GetMaximum();
	int binmax = hCe->GetMaximumBin();
	double xmax = hCe->GetXaxis()->GetBinCenter(binmax);
	floo->SetParameters(max,xmax,10,max*0.1,50,20,5);
	hCe->Fit(floo,"Q","",lo,hi);
}

void timeWindow( TFile *infile, double w0, double w1, double w2, double w3)
{
	// w0 -> w1 = time window around peak of time-expt
	// w2 -> w0 = window to left of w01
	// w1 -> w3 = window to right of w01

	// grab trees from infile
	TTree *dirc = infile->Get("dirc");
	TTree *reco = infile->Get("reco");

	// get lens, bar angle, and beam momentum
	int lens(0);
	double alpha(0),beam(0);
	dirc->SetBranchAddress("lens",&lens);
	dirc->SetBranchAddress("theta",&alpha);
	dirc->GetEntry(0);

	// define histograms of diff, outside diff, and counters
	TH1D *diff = new TH1D("diff","diff",1000,-50,50);
	TH1D *dout = new TH1D("dout","dout",1000,-50,50);
	TH1D *tof1 = new TH1D("time-tof1","time-tof1",1000,0,1000);
	TH1D *tof2 = new TH1D("time-tof2","time-tof2",1000,0,1000);
	TH1D *trig = new TH1D("time-trig","time-trig",1000,0,1000);

	// project to histograms
	reco->Project("diff","diff");
	reco->Project("dout","diff");
	reco->Project("time-tof1","time-tof1");
	reco->Project("time-tof2","time-tof2");
	reco->Project("time-trig","time-trig");

	// get diff value at peak and define initial cut
	double diffpeak = diff->GetMaximumBin();
	diffpeak = (diffpeak-500)/10;
	char *diffcut = Form("%f<diff-%f && diff-%f<%f",w0,diffpeak,diffpeak,w1);
	cout << "Initial cut\t" << diffcut << endl;
	
	
}

//Fit a peak in an energy spectrum
void FitPeak(int run, int decay, int shaping, int top) {
	TH1D* hCe = (TH1D*) gROOT->FindObject("hCe");
	if (hCe != 0) delete hCe;
	hCe = new TH1D("hCe","hCe",3000,0,6000);

	trapf.SetPath("Files/Dec15/Fixed");
	if (!trapf.Open(run,decay,shaping,top)) return;
	for (int i = 0; i < trapf.GetNumEvents(); i++) {
		trapf.GetEvent(i);
		if (trapf.Trap_event.ch == 9 && trapf.Trap_event.AveE > 0) {
			hCe->Fill(trapf.Trap_event.AveE);
		}
	}
	gROOT->cd();
	hCe->GetXaxis()->SetRangeUser(abslo,abshi); //avoid low energy stuffs
	hCe->Draw();
	double max = hCe->GetMaximum();
	int binmax = hCe->GetMaximumBin();
	double xmax = hCe->GetXaxis()->GetBinCenter(binmax);
	TF1* floo = (TF1*) gROOT->FindObject("floo");
	if (floo != 0) delete floo;
	floo = new TF1("floo", floof, 0, 6000, 7);
	floo->SetParameters(max,xmax,10,max*0.1,50,20,5);
	floo->SetParLimits(0,0,2*max);
	floo->SetParLimits(3,0,2*max);
	floo->SetParLimits(4,0,500);
	floo->SetParName(2,"sigma");
	hCe->Fit(floo,"Q","",1000,xmax+60);
	double mean, fwhm;
	double lo, hi;
	GetFWHM(floo,mean,fwhm);
	lo = mean-8.5*fwhm; hi = mean+3*fwhm;
	if (lo < abslo) lo = abslo;
	ReFit(lo,hi);
	GetFWHM(floo,mean,fwhm);
	lo = mean-8.5*fwhm; hi = mean+3*fwhm;
	if (lo < abslo) lo = abslo;
	ReFit(lo,hi);
	GetFWHM(floo,mean,fwhm);
	lo = mean-8.5*fwhm; hi = mean+3*fwhm;
	if (lo < abslo) lo = abslo;
	ReFit(lo,hi);
	GetFWHM(floo,mean,fwhm);
	//draw parts
	TF1* part1 = (TF1*) gROOT->FindObject("part1");
	if (part1 != 0) delete part1;
	part1 = new TF1("part1", gaussian, 0, 6000, 3);
	for (int i=0;i<3;i++)
		part1->SetParameter(i,floo->GetParameter(i));
	part1->SetLineColor(kViolet);
	part1->Draw("same");
	TF1* part2 = (TF1*) gROOT->FindObject("part2");
	if (part2 != 0) delete part2;
	part2 = new TF1("part2", hypermet, 0, 6000, 4);
	part2->SetParameter(0,floo->GetParameter(3));
	part2->SetParameter(1,floo->GetParameter(1));
	part2->SetParameter(2,floo->GetParameter(2));
	part2->SetParameter(3,floo->GetParameter(4));
	part2->SetLineColor(kBlue);
	part2->Draw("same");
	TF1* part3 = (TF1*) gROOT->FindObject("part3");
	if (part3!= 0) delete part3;
	part3 = new TF1("part3", stepbkgd, 0, 6000, 4);
	part3->SetParameter(0,floo->GetParameter(5));
	part3->SetParameter(1,floo->GetParameter(6));
	part3->SetParameter(2,floo->GetParameter(1));
	part3->SetParameter(3,floo->GetParameter(2));
	part3->SetLineColor(kGreen);
	part3->Draw("same");
	fgt->SetLineColor(kBlack);
	fgt->Draw("same");
}

//.........这里部分代码省略.........
				newBaseHist->SetBinContent(iBin,baseHist->GetBinContent(iBin));
				newBaseHist->SetBinError(iBin,baseHist->GetBinError(iBin));
				newDMHist->SetBinContent(iBin,profileHist1[title_id][0]->GetBinContent(iBin));
				newDMHist->SetBinError(iBin,profileHist1[title_id][0]->GetBinError(iBin));
				newOverHist->SetBinContent(iBin,overHist->GetBinContent(iBin));
				newOverHist->SetBinError(iBin,overHist->GetBinError(iBin));
				newDMHist2->SetBinContent(iBin,profileHist1[title_id][1]->GetBinContent(iBin));
				newDMHist2->SetBinError(iBin,profileHist1[title_id][1]->GetBinError(iBin));	
				
			 }
			 newBaseHist->Multiply(newDMHist);
			 newOverHist->Multiply(newDMHist2);
			 baseHist=newBaseHist;
			 overHist=newOverHist;
			 delete newDMHist;
			 delete newDMHist2;
		 }
	
		 if(title_id==0) baseHist->GetXaxis()->SetRangeUser(0,100);
		 baseHist->GetXaxis()->SetLabelOffset(0.1);		
		 baseHist->GetXaxis()->SetTitle(XTitle[title_id]);
		 baseHist->GetYaxis()->SetTitle("Fake rate");
		 baseHist->SetMarkerStyle(20);
		 baseHist->SetMarkerSize(2);
		 baseHist->SetLineWidth(2);
		 
		 overHist->SetMarkerColor(kRed);
		 overHist->SetMarkerStyle(25);
		 overHist->SetMarkerSize(2);
		 overHist->SetLineWidth(2);
		 overHist->SetLineStyle(2);
		 overHist->SetLineColor(kRed);
		
		 int binmax = baseHist->GetMaximumBin();
	     double max = 1.2*(baseHist->GetBinContent(binmax)+baseHist->GetBinError(binmax));
	     int binmax2 = overHist->GetMaximumBin();
	     double max2 = 1.2*(overHist->GetBinContent(binmax2)+overHist->GetBinError(binmax2));
	     max = max2 > max ? max2 : max;
	 
	     baseHist->SetMaximum(10.0);
	     baseHist->SetMinimum(1e-3);
	     
		 
		 histPad->cd();	
			 histPad->SetLogy();
		 baseHist->Draw();
		 overHist->Draw("same");	
		 
		TLegend* leg = new TLegend(0.8,0.80,0.9,0.9,NULL,"brNDC");
		leg->SetFillColor(0);
		leg->SetTextSize(0.035);
		leg->SetBorderSize(0);
		
		leg->AddEntry(baseHist, leg1, "lp");
		leg->AddEntry(overHist, leg2, "lp");
		
		leg->Draw();
		 //~ 
		 
    TString lumist="";
	TPaveText *ll = new TPaveText(0.15, 0.95, 0.95, 0.99, "NDC");
	ll->SetTextSize(0.03);
	ll->SetTextFont(62);
	ll->SetFillColor(0);
	ll->SetBorderSize(0);
	ll->SetMargin(0.01);

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

for(int q =1; q<ttgamma->GetNbinsX(); q++){
mcEff->GetXaxis()->SetBinLabel(q, step[q]);
dataEff->SetBinContent(q, data->GetBinContent(q+1)/data->GetBinContent(q));
dataEff->SetBinError(q, dataEff->GetBinContent(q)*sqrt(pow(data->GetBinError(q+1)/data->GetBinContent(q+1),2)+pow(data->GetBinError(q)/data->GetBinContent(q),2)));
mcEff->SetBinContent(q, allMC->GetBinContent(q+1)/allMC->GetBinContent(q));

}  

	
	TCanvas *c2 = new TCanvas("cutflow eff","cutflow eff",600, 500);
  	
	mcEff->SetLineColor(kRed);
	mcEff->Draw();
	dataEff->Draw("Esame");
	
	TLegend *tleg3;
	tleg3 = new TLegend(0.7,0.7,0.8,0.9);
	tleg3->SetTextSize(0.04);
	tleg3->SetBorderSize(0);
	tleg3->SetFillColor(10);
	tleg3->AddEntry(dataEff , "data", "l");
	tleg3->AddEntry(mcEff , "mc", "l");
	
	tleg3->Draw("same");
	
	c2->SaveAs("plots/cutFlow/cutEff_"+Variable+".pdf");
	c2->SaveAs("plots/cutFlow/cutEff_"+Variable+".png");
	delete c2;

  //draw histos to files
  TCanvas *c1 = new TCanvas("cutflow","cutflow",900, 600);

  hs->SetMaximum(data->GetBinContent(data->GetMaximumBin())*1.3);
  
  if(logPlot ==true)
  hs->SetMinimum(10.);

  hs->Draw("hist");
  
  for(int q =0; q<ttgamma->GetNbinsX(); q++)
  hs->GetXaxis()->SetBinLabel(q+1, step[q]);
  
  data->Draw("E same");
  data->SetMarkerStyle(20);
  
//  hs->GetXaxis()->SetLimits(MinX, MaxX);
  hs->GetXaxis()->SetTitle(Xtitle); hs->GetXaxis()->SetTitleSize(0.05);
  hs->GetYaxis()->SetTitle("Number of Events");hs->GetYaxis()->SetTitleSize(0.05);
  
  
  	TLegend *tleg2;
	tleg2 = new TLegend(0.7,0.7,0.8,0.9);
	tleg2->SetTextSize(0.04);
	tleg2->SetBorderSize(0);
	tleg2->SetFillColor(10);
	tleg2->AddEntry(data , "2012 data", "lpe");
	tleg2->AddEntry(ttgamma , "t#bar{t}+gamma", "lf");
	tleg2->AddEntry(tt , "t#bar{t}", "lf");
	tleg2->AddEntry(T_tW, "single top", "lf");
	tleg2->AddEntry(DY1 , "Drell-Yan", "lf");
	tleg2->AddEntry(ZZ , "Di-boson", "lf");
	tleg2->AddEntry(wjets , "w+jets", "lf");
	tleg2->AddEntry(QCD_all , "QCD", "lf");
	
	//tleg2->AddEntry(singtEff, "single-t"      , "l");

//.........这里部分代码省略.........
  qcd1->Add(qcd5);
  qcd1->Add(qcd6);
  qcd1->Add(qcd7);
  qcd1->Add(qcd8);
  qcd1->Add(qcd9);
  qcd1->Add(qcd10);
  qcd1->Add(qcd11);
  }
  
TH1D* qcd_data = getQCD(Obj, Variable, rebinFact);
qcd_data->Scale(qcd1->Integral());

THStack *hs = new THStack("hs","test");
  if(inclQ == true){
  hs->Add(qcd_data);
  }else{
  hs->Add(qcd_data);
  }
    
  zjets->Add(z2jets);
  zjets->Add(z3jets);
  zjets->Add(z4jets);  

  wjets->Add(w2jets);
  wjets->Add(w3jets);
  wjets->Add(w4jets);
  
  hs->Add(zjets);
  hs->Add(wjets);
      
  hs->Add(top_t);
  hs->Add(top_tw);
  hs->Add(top_s);
  hs->Add(tbar_t);
  hs->Add(tbar_tw);
  hs->Add(tbar_s);
  
  hs->Add(tt);
  
  TH1D* allMC = (TH1D*)tt->Clone("allMC");
  allMC->Add(sing_top); allMC->Add(wjets); allMC->Add(zjets); allMC->Add(qcd1);
  
  cout << "tt: " << tt->Integral()  << " wjets: " << wjets->Integral() << " zjets: " << zjets->Integral() << " single top: " << sing_top->Integral() << " qcd: " << qcd1->Integral() << " all: " << allMC->Integral() << " data: " << data->Integral() << endl;
  	
  //draw histos to files
  TCanvas *c1 = new TCanvas("Plot","Plot",900, 600);

  hs->SetMaximum(data->GetBinContent(data->GetMaximumBin())*1.3);

  hs->Draw();
  data->Draw("E same");
  data->SetMarkerStyle(20);
  
  hs->GetXaxis()->SetLimits(MinX, MaxX);
  hs->GetXaxis()->SetTitle(Xtitle); hs->GetXaxis()->SetTitleSize(0.05);
  hs->GetYaxis()->SetTitle("Number of Events");hs->GetYaxis()->SetTitleSize(0.05);
  
  
  	TLegend *tleg2;
	tleg2 = new TLegend(0.6,0.6,0.8,0.9);
	tleg2->SetTextSize(0.04);
	tleg2->SetBorderSize(0);
	tleg2->SetFillColor(10);
	tleg2->AddEntry(data , "2012 data", "lpe");
	tleg2->AddEntry(tt , "t#bar{t} l+jets", "lf");
	//tleg2->AddEntry(tt_lep , "t#bar{t} di-lepton", "lf");
	tleg2->AddEntry(top_t, "single top", "lf");
	tleg2->AddEntry(wjets , "w+jets", "lf");
	tleg2->AddEntry(zjets , "z+jets", "lf");
	tleg2->AddEntry(qcd , "QCD", "lf");

	//tleg2->AddEntry(singtEff, "single-t"      , "l");
	//tleg2->AddEntry(singtwEff, "single-tW"      , "l");
 	tleg2->Draw("same");	

	TText* textPrelim = doPrelim(0.12,0.96);
	textPrelim->Draw();

  if(logPlot ==true){
  c1->SetLogy();
  }	
  
  TString plotName("plots/Control/Muon/");
  
  if(logPlot ==true){
    plotName += Variable+"Log";
    plotName += Nbtags+".pdf";
    
  }else{
    plotName += Variable;  
    plotName += Nbtags+".pdf";
  }
 
 
  c1->SaveAs(plotName);
  delete c1;
  
  }
  	
}

void DoTDeriMax1090Correction(TString SpectrumFileInput = "/lustre/miifs05/scratch/him-specf/hyp/steinen/COSYBeamtestAna/COSYnewMogon/June2014/COSYJune2014Dataset11_200,100,0,5339_SR1.root", TString FitFileInput = "/lustre/miifs05/scratch/him-specf/hyp/steinen/COSYBeamtestAna/COSYnewMogon/Fit/FitCOSYJune2014Dataset11_200,100,0,5339_SR1.root")
{
	TH2D *hSpectrumTDeriMax1090_EnergyChannel;
	TH2D *hSpectrumTDeriMax1090Rel_EnergyChannel;
	
	TH2D *hSpectrumT1090_EnergyChannelCorr1;
	TNtuple *DataNTuple;
	TFile *SpectrumInput = new TFile(SpectrumFileInput.Data());
	
	hSpectrumTDeriMax1090_EnergyChannel = (TH2D*) SpectrumInput->Get("Histograms/Energy_DeriMaxT90/Energy_DeriMaxT90_01");
	hSpectrumTDeriMax1090_EnergyChannel->SetDirectory(0);
	hSpectrumTDeriMax1090Rel_EnergyChannel = (TH2D*) SpectrumInput->Get("Histograms/Energy_DeriMaxT90Rel/Energy_DeriMaxT90Rel_01");
	hSpectrumTDeriMax1090Rel_EnergyChannel->SetDirectory(0);
	hSpectrumT1090_EnergyChannelCorr1 = (TH2D*) SpectrumInput->Get("Histograms/EnergyRt1090/EnergyRt1090CorrectionRt_01");
	hSpectrumT1090_EnergyChannelCorr1->SetDirectory(0);
	SpectrumInput->Close();
	//hSpectrumTDeriMax1090_EnergyChannel->Draw("colz");
	
	
	TFile *FitInput = new TFile(FitFileInput.Data(),"Update");
	DataNTuple = (TNtuple*)FitInput->Get("DataNTuple");
	DataNTuple->Scan();
	
	Int_t entries = (Int_t)DataNTuple->GetEntries();
	cout<<"Number of Entries: "<<entries<< endl;
	const int entriesArrayValue =entries;
	TF1 *FitFunc[entriesArrayValue];
	float Energy=0;
	float ChannelPeakPos=0;
	float ChannelRangeMin=0;
	float ChannelRangeMax=0;
	DataNTuple->SetBranchAddress("Energy",&Energy);
	DataNTuple->SetBranchAddress("ChannelPeakPos",&ChannelPeakPos);
	DataNTuple->SetBranchAddress("ChannelRangeMin",&ChannelRangeMin);
	DataNTuple->SetBranchAddress("ChannelRangeMax",&ChannelRangeMax);
	
	TCanvas* can=new TCanvas();
	gPad->SetLogz();
	hSpectrumTDeriMax1090_EnergyChannel->GetXaxis()->SetRangeUser(-30,250);
	hSpectrumTDeriMax1090Rel_EnergyChannel->GetXaxis()->SetRangeUser(-0.1,0.95);
	
	for (Int_t ki=0;ki<entries;ki++)
	{
		DataNTuple->GetEntry(ki);
		//if (int(Energy) == 1332)
		//if (int(Energy) == 510)
		//if (ki == entries-1)
		{
			cout << ChannelRangeMin << " " << ChannelRangeMax << endl;
			
			
			//first correction via TDeriMaxT90Rel
			hSpectrumTDeriMax1090_EnergyChannel->GetYaxis()->SetRangeUser(ChannelRangeMin,ChannelRangeMax);
			hSpectrumTDeriMax1090Rel_EnergyChannel->GetYaxis()->SetRangeUser(ChannelRangeMin,ChannelRangeMax);
			//TF1* FitFuncSlices = new TF1("FitFuncSlices","gaus(0)+pol0(3)",ChannelRangeMin,ChannelRangeMax);
			//FitFuncSlices->SetParameters(1000,ChannelPeakPos-10,4,10);
			//FitFuncSlices->SetParLimits(1,ChannelRangeMin,ChannelRangeMax);
			//FitFuncSlices->SetParLimits(2,0,10);
			//FitFuncSlices->SetParLimits(3,0,100);
	
			//gDirectory->ls();
			
			TH1D *hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually=new TH1D("hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually","",hSpectrumTDeriMax1090Rel_EnergyChannel->GetNbinsX(),-0.3,1.3);
			//cout <<hSpectrumTDeriMax1090_EnergyChannel_MaxPos->GetEntries()<< endl;
			for(int binX = hSpectrumTDeriMax1090Rel_EnergyChannel->GetXaxis()->FindBin(-0.1);binX <= hSpectrumTDeriMax1090Rel_EnergyChannel->GetXaxis()->FindBin(0.90);binX++)
			{
				cout << "binx " << binX << endl;
				TH1D *hProfileY =hSpectrumTDeriMax1090Rel_EnergyChannel->ProjectionY("_py",binX,binX);	
				double MaxValue=hProfileY->GetBinCenter(hProfileY->GetMaximumBin());
				TF1* FitFuncSlices = new TF1("FitFuncSlices","gaus(0)+[3]",MaxValue-20,MaxValue+20);
				TF1* FitFuncGausSlices = new TF1("FitFuncGausSlices","gaus(0)",MaxValue-20,MaxValue+20);
				FitFuncGausSlices->SetParameters(hProfileY->GetBinContent(hProfileY->GetMaximumBin()),MaxValue,4);
				
				hProfileY->Fit(FitFuncGausSlices,"RNQ");
				
				FitFuncSlices->SetParameters(FitFuncGausSlices->GetParameter(0),FitFuncGausSlices->GetParameter(1),FitFuncGausSlices->GetParameter(2),10);
				
				FitFuncSlices->SetParLimits(0,0,10000);
				FitFuncSlices->SetParLimits(1,MaxValue-10,MaxValue+10);
				FitFuncSlices->SetParLimits(2,0,10);
				FitFuncSlices->SetParLimits(3,0,100);
				hProfileY->Fit(FitFuncSlices,"RNQ");
				
				cout <<MaxValue<<"  " << FitFuncSlices->GetParameter(1) << "   " << FitFuncSlices->GetParError(1) <<endl;
				hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->SetBinContent(binX, FitFuncSlices->GetParameter(1));
				hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->SetBinError(binX, FitFuncSlices->GetParError(1));
				
				
			}
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->GetYaxis()->SetRangeUser(ChannelPeakPos-100,ChannelPeakPos+50);
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->GetXaxis()->SetRangeUser(-0.05,0.9);
			TF1 * funcCorrConst = new TF1("funcCorrConst","pol0",-0.03,0.03);
			funcCorrConst->SetLineColor(kRed);
			TF1 * funcCorrPol = new TF1("funcCorrPol",PolPiecewise,-0.05,0.9,6);
			
			funcCorrPol->SetLineColor(kBlue);
			funcCorrPol->SetParameter(0,0.04);
			funcCorrPol->SetParameter(1,ChannelPeakPos);
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->Fit(funcCorrConst,"R");
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->Fit(funcCorrPol,"R+");
//.........这里部分代码省略.........

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

  hs->Add(ttgamma);
  

std::cout << "Total number of events" << std::endl;
std::cout << "Data: " << data->Integral() << std::endl;
std::cout << "ttbar + gamma: " << ttgamma->Integral() << std::endl;
std::cout << "Inclusive ttbar: " << tt->Integral() << std::endl;
std::cout << "W+Jets: " << wjets->Integral() << std::endl;
std::cout << "Drell-Yan1: " << DY1->Integral() << std::endl;
std::cout << "Drell-Yan2: " << DY2->Integral() << std::endl;
std::cout << "Single top: " << T_tW->Integral() << std::endl;
std::cout << "Anti-Single top: " << Tbar_tW->Integral() << std::endl;
std::cout << "ZZ: " << ZZ->Integral() << std::endl;
std::cout << "WW: " << WW->Integral() << std::endl;
std::cout << "WZ: " << WZ->Integral() << std::endl;
std::cout << "QCD: " << QCD_all->Integral() << std::endl;


  //draw histos to files
  TCanvas *c1 = new TCanvas("Plot","Plot",900, 600);
  TPad *pad1 = new TPad("pad1","pad1",0,0.3,1,1);
  pad1->SetBottomMargin(0);
  pad1->Draw();
  pad1->cd();	  

  data->Draw("e x0");
  data->SetAxisRange(MinX, MaxX, "X");

  hs->Draw("hist");
  allMC->Draw("same e2");
	
  hs->SetMaximum(data->GetBinContent(data->GetMaximumBin())*1.3);

  data->Draw("E same");
  data->SetMarkerStyle(20);
  
  hs->GetXaxis()->SetLimits(MinX, MaxX);
  hs->GetXaxis()->SetTitle(Xtitle); hs->GetXaxis()->SetTitleSize(0.05);
  hs->GetYaxis()->SetTitle("Number of Events");hs->GetYaxis()->SetTitleSize(0.05);
   
  
  	TLegend *tleg2;
	tleg2 = new TLegend(0.7, 0.6, 0.8, 0.9);
	tleg2->SetTextSize(0.04);
	tleg2->SetBorderSize(0);
	tleg2->SetFillColor(10);
	tleg2->AddEntry(data , "2012 data", "lpe");
	tleg2->AddEntry(ttgamma , "t#bar{t}#gamma", "lf");

	if(inclTop == false){
        tleg2->AddEntry(tt_sig , "t#bar{t}#gamma (MG)", "lf");
        tleg2->AddEntry(tt_lep , "t#bar{t} l #rightarrow #gamma (MG)", "lf");
        tleg2->AddEntry(tt , "t#bar{t} fake (MG)", "lf");
        }else{
        tleg2->AddEntry(tt , "t#bar{t}", "lf");
        }
	tleg2->AddEntry(T_tW, "Single Top"      , "lf");
//	tleg2->AddEntry(Tbar_tW, "anti-single-tW"      , "lf");	
	tleg2->AddEntry(DY1 , "Z+Jets", "lf");
//	tleg2->AddEntry(DY2 , "DYJetsToLL", "lf");
	tleg2->AddEntry(ZZ, "Diboson", "lf");
//	tleg2->AddEntry(WW, "WW", "lf");
//	tleg2->AddEntry(WZ , "WZ", "lf");
	tleg2->AddEntry(wjets , "W+Jets", "lf");

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


  cdEdxPure->cd(1); 
  f2dHistdEdxSignalPionFromLambda->Draw("colz"); 
  cdEdxPure->cd(2); 
  f2dHistResponsePionFromLambda->Draw("colz"); 
  cdEdxPure->cd(3); 
  fLowPtPionResponse->Draw();
  cdEdxPure->cd(4); 
  fMidPtPionResponse->Draw();

  cdEdxPure->cd(5); 
  f2dHistdEdxSignalProtonFromLambda->Draw("colz"); 
  cdEdxPure->cd(6); 
  f2dHistResponseProtonFromLambda->Draw("colz"); 
  cdEdxPure->cd(7); 
  fLowPtProtonResponse->Draw();
  cdEdxPure->cd(8); 
  fMidPtProtonResponse->Draw();



  //Write explanations on canvases
  cdEdxPure->cd(1); Tl.DrawLatex(.25, .9277, "#pi^{-} from #Lambda: TPC Signal");
  cdEdxPure->cd(2); Tl.DrawLatex(.15, .9277, "#pi^{-} from #Lambda: AliPIDResponse Value");
  cdEdxPure->cd(3); Tl.DrawLatex(.21, .9277, "#pi^{-} N#sigma, 0.5 < p_{T} (GeV/c) < 1.0");
  cdEdxPure->cd(4); Tl.DrawLatex(.21, .9277, "#pi^{-} N#sigma, 1.5 < p_{T} (GeV/c) < 2.5");

  cdEdxPure->cd(5); Tl.DrawLatex(.25, .9277, "p from #Lambda: TPC Signal");
  cdEdxPure->cd(6); Tl.DrawLatex(.15, .9277, "p from #Lambda: AliPIDResponse Value");
  cdEdxPure->cd(7); Tl.DrawLatex(.21, .9277, "p N#sigma, 0.5 < p_{T} (GeV/c) < 1.0");
  cdEdxPure->cd(8); Tl.DrawLatex(.21, .9277, "p N#sigma, 1.5 < p_{T} (GeV/c) < 2.5");

  Double_t lLowPtPeakPion         = fLowPtPionResponse->GetBinCenter( fLowPtPionResponse->GetMaximumBin() );
  Double_t lMidPtPeakPion         = fMidPtPionResponse->GetBinCenter( fMidPtPionResponse->GetMaximumBin() );
  Double_t lLowPtPeakProton       = fLowPtProtonResponse->GetBinCenter( fLowPtProtonResponse->GetMaximumBin() );
  Double_t lMidPtPeakProton       = fMidPtProtonResponse->GetBinCenter( fMidPtProtonResponse->GetMaximumBin() );

  //List Maximal Values
  cout<<"Maximal Value for pion from Lambda at low pt...............: " <<lLowPtPeakPion<<endl;
  cout<<"Maximal Value for pion from Lambda at mid pt...............: " <<lMidPtPeakPion<<endl;
  cout<<"Maximal Value for proton from Lambda at low pt.............: " <<lLowPtPeakProton<<endl;
  cout<<"Maximal Value for proton from Lambda at mid pt.............: " <<lMidPtPeakProton<<endl;

  if( TMath::Abs( lLowPtPeakPion ) > 0.3)   cout<<"*** WARNING: Check Low Pt pion PID Response! ***"<<endl;
  if( TMath::Abs( lMidPtPeakPion ) > 0.3)   cout<<"*** WARNING: Check Mid Pt pion PID Response! ***"<<endl;
  if( TMath::Abs( lLowPtPeakProton ) > 0.3) cout<<"*** WARNING: Check Low Pt proton PID Response! ***"<<endl;
  if( TMath::Abs( lMidPtPeakProton ) > 0.3) cout<<"*** WARNING: Check Mid Pt proton PID Response! ***"<<endl;

  TLatex Tlq;
  Tlq.SetNDC();
  Tlq.SetTextSize(0.06);

  //Draw Arrows to be sure!
  Double_t lFractionHeight = 0.33;
  cdEdxPure->cd(3);  
  TArrow *ar1 = new TArrow(lLowPtPeakPion,lFractionHeight * fLowPtPionResponse->GetMaximum(),lLowPtPeakPion,0.0, 0.008 ,"|>");
  ar1->SetLineWidth(2);
  ar1->Draw();
  if( TMath::Abs( lLowPtPeakPion ) < 0.3) {
    Tlq.SetTextColor(8); 
    Tlq.DrawLatex ( 0.15,0.8, "OK!");
  } else { 
    Tlq.SetTextColor(kRed); 
    Tlq.DrawLatex ( 0.15,0.8, "Prob.! Check dEd/x!");
  }

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

wjets->Add(w2jets);
wjets->Add(w3jets);
wjets->Add(w4jets);

//make combined top and single top template
TH1D* qcd_all = (TH1D*)qcd->Clone("qcd_all");
qcd_all->Add(qcd2);
qcd_all->Add(qcd3);
qcd_all->Add(qcd4);
qcd_all->Add(qcd5);
qcd_all->Add(qcd6);
qcd_all->Add(qcd7);
qcd_all->Add(qcd8);
qcd_all->Add(qcd9);
qcd_all->Add(qcd10);
qcd_all->Add(qcd11);
  
THStack *hs = new THStack("hs","test");

qcd_all->SetLineColor(kBlack);
zjets->SetLineColor(kBlack);
wjets->SetLineColor(kBlack);	
sing_top->SetLineColor(kBlack);
tt->SetLineColor(kBlack);	

  hs->Add(qcd_all);
  hs->Add(zjets);
  hs->Add(wjets);
  hs->Add(sing_top);
  hs->Add(tt);

  //draw histos to files
  TCanvas *c1 = new TCanvas("Plot","Plot",900, 600);
		
  hs->SetMaximum(data->GetBinContent(data->GetMaximumBin())*1.2);
  hs->SetMinimum(1.);
  
  hs->Draw();
  data->Draw("E same");
  data->SetMarkerStyle(20);

//events:
cout << "ttbar: " << tt->Integral() << endl;
cout << "data: " << data->Integral() << endl;
  
  if(logPlot == true){
  hs->GetXaxis()->SetLimits(MinX, MaxX);
  }else{
  hs->GetXaxis()->SetLimits(MinX, 4.5);
  }
  
  hs->GetXaxis()->SetTitle(Xtitle); hs->GetXaxis()->SetTitleSize(0.05);
  hs->GetYaxis()->SetTitle("Number of Events");hs->GetYaxis()->SetTitleSize(0.05);
  
  
  	TLegend *tleg2;
	tleg2 = new TLegend(0.7,0.7,0.8,0.9);
	tleg2->SetTextSize(0.04);
	tleg2->SetBorderSize(0);
	tleg2->SetFillColor(10);
	tleg2->AddEntry(data , "2012 data", "lpe");
	tleg2->AddEntry(tt , "t#bar{t}", "f");
	tleg2->AddEntry(sing_top, "single top", "f");
	tleg2->AddEntry(wjets , "w+jets", "f");
	tleg2->AddEntry(zjets , "z+jets", "f");
	tleg2->AddEntry(qcd_all , "QCD", "f");
	
	//tleg2->AddEntry(singtEff, "single-t"      , "l");
	//tleg2->AddEntry(singtwEff, "single-tW"      , "l");
 	tleg2->Draw("same");	
	
	if(logPlot == true){
  	TText* textPrelim = doPrelim(0.12,0.96, "#geq 0 btags");
	textPrelim->Draw();
  	}else{
	TText* textPrelim = doPrelim(0.2,0.96, "");
	textPrelim->Draw();
	}
	
  if(logPlot ==true){
  c1->SetLogy();
  }	
  
  TString plotName("Plots/TTbarPlusVarAnalysis/Nbtags/");
  
  if(logPlot ==true){
    plotName += Variable+"_Log";  
  }else{
    plotName += Variable;  
  }
 
 
  c1->SaveAs(plotName+".pdf");
  c1->SaveAs(plotName+".png");
  delete c1;
  
  }
  	
}

//.........这里部分代码省略.........
 		
dataEff->Sumw2();
for(int i =1; i<11; i++){
mcEff->GetXaxis()->SetBinLabel(i, step[i]);
dataEff->SetBinContent(i, data->GetBinContent(i+1)/data->GetBinContent(i));
dataEff->SetBinError(i, sqrt(pow(data->GetBinContent(i+1),2)+pow(data->GetBinContent(i),2)));
mcEff->SetBinContent(i, allMC->GetBinContent(i+1)/allMC->GetBinContent(i));

}  

	
	TCanvas *c2 = new TCanvas("cutflow eff","cutflow eff",600, 500);
  	
	mcEff->SetLineColor(kRed);
	mcEff->Draw();
	dataEff->Draw("Esame");
	
	TLegend *tleg3;
	tleg3 = new TLegend(0.7,0.7,0.8,0.9);
	tleg3->SetTextSize(0.04);
	tleg3->SetBorderSize(0);
	tleg3->SetFillColor(10);
	tleg3->AddEntry(dataEff , "data", "l");
	tleg3->AddEntry(mcEff , "mc", "l");
	
	tleg3->Draw("same");
	
	c2->SaveAs("plots/cutFlow/cutEff_data.png");
	delete c2;

  //draw histos to files
  TCanvas *c1 = new TCanvas("cutflow","cutflow",900, 600);

  hs->SetMaximum(data->GetBinContent(data->GetMaximumBin())*1.3);
  
  if(logPlot ==true)
  hs->SetMinimum(10000.);

  hs->Draw();
  for(int i =0; i<tt->GetNbinsX(); i++)
  hs->GetXaxis()->SetBinLabel(i+1, step[i]);
  
  data->Draw("E same");
  data->SetMarkerStyle(20);
  
//  hs->GetXaxis()->SetLimits(MinX, MaxX);
  hs->GetXaxis()->SetTitle(Xtitle); hs->GetXaxis()->SetTitleSize(0.05);
  hs->GetYaxis()->SetTitle("Number of Events");hs->GetYaxis()->SetTitleSize(0.05);
  
  
  	TLegend *tleg2;
	tleg2 = new TLegend(0.7,0.7,0.8,0.9);
	tleg2->SetTextSize(0.04);
	tleg2->SetBorderSize(0);
	tleg2->SetFillColor(10);
	tleg2->AddEntry(data , "2012 data", "lpe");
	tleg2->AddEntry(tt , "t#bar{t}", "lf");
	tleg2->AddEntry(top_t, "single top", "lf");
	tleg2->AddEntry(wjets , "w+jets", "lf");
	tleg2->AddEntry(zjets , "z+jets", "lf");
	tleg2->AddEntry(qcd , "QCD", "lf");
	
	//tleg2->AddEntry(singtEff, "single-t"      , "l");
	//tleg2->AddEntry(singtwEff, "single-tW"      , "l");
 	tleg2->Draw("same");	
	

void reconstruction(TString infile = "../build/reco_spr.root", 
					int num_events = 10000, 
					char *runtype = "data", 
					double tdiff = 1.5)
{
	gErrorIgnoreLevel = kWarning; // ignore 'Info in...' messages
	gStyle->SetOptFit(1111); // show fit parameters
	gStyle->SetOptStat(0);   // remove stats box
	TGaxis::SetMaxDigits(4); // put axis in scientific notation

	TFile *f = TFile::Open(infile); 
	TTree *t = (TTree*)f->Get("reco"); //Lee's Tree
	TTree *d = (TTree*)f->Get("dirc"); //Roman's Tree

	// set value of lens and alpha from dirc tree
	int lens(0);
	double angle(0),beam(0),yield(0),nph(0);
	d->SetBranchAddress("lens",&lens);
	d->SetBranchAddress("theta",&angle);
	d->SetBranchAddress("beam",&beam);
	d->SetBranchAddress("yield",&yield);
	d->SetBranchAddress("nph",&nph);
	d->GetEntry(0);

	// get time-expt peak value
	TH1D *diff = new TH1D("diff","diff",100,-5,5);
	t->Project("diff","diff");
	double diffmax = diff->GetBinCenter(diff->GetMaximumBin());

	// define cuts and save path
	char *sim_cut  = Form("abs(diff-%f)<%f",diffmax,tdiff);
	char *data_cut = Form("abs(diff-%f)<%f",diffmax,tdiff);
	char *savepath = "../macro/timecuts";

	// define cut and project from tree to histogram w/ proper cuts
	theta = new TH1D("theta","theta",100,0.6,1.0);
	if(runtype=="data") // do all the data things
	{
		// make counter histograms for cut
		TH1D *tof1 = new TH1D("time-tof1","time-tof1",1000,0,1000);
		TH1D *tof2 = new TH1D("time-tof2","time-tof2",1000,0,1000);
		TH1D *trig = new TH1D("time-trig","time-trig",1000,0,1000);

		// project onto counter histograms
		t->Project("time-tof1","time-tof1");
		t->Project("time-tof2","time-tof2");
		t->Project("time-trig","time-trig");
	
		// make cuts w/ counters
		cntrCut(tof1, data_cut, "time-tof1");
		cntrCut(tof2, data_cut, "time-tof2");
		cntrCut(trig, data_cut, "time-trig");
		//cout << "data cut\t" << data_cut << endl;

		/*
		// make counter hists pretty (and distinguishable)
		tof1->SetTitle(data_cut);
		tof1->SetTitleSize(0.5);
		tof1->SetLineColor(kBlue);
		tof2->SetLineColor(kRed);
		trig->SetLineColor(kGreen);

		// draw histograms and save canvas
		char *cname = Form("%s/reco_%s_%d_%d_ccuts.png",
					  savepath,runtype,lens,(int)angle);
		TCanvas *ccanv = new TCanvas();
		gPad->SetLogy();
		tof1->Draw();
		tof2->Draw("same");
		trig->Draw("same");
		ccanv->Print(cname);
		*/

		t->Project("theta","theta",data_cut);
	}
	else
		t->Project("theta","theta",sim_cut);

	// preliminary paramters for fit of theta
	double height = theta->GetMaximum();
	double center = theta->GetBinCenter(theta->GetMaximumBin());
	double sigma  = 0.015;
	double slope  = 1.0;
	double shift  = height/3;

	// define custom fit, set parameters, and fit to data
	TF1 *gaus0 = new TF1("gaus0", "gaus+[3]*x+[4]", 
						 center-0.05, center+0.05);
	gaus0->SetParameters(height,center,sigma,slope,shift);
	gaus0->SetParNames("height","mean","sigma","slope","const");
	theta->Fit("gaus0","QR");
	double truesig = 1000*gaus0->GetParameter("sigma");
	
	// define strings for picture and root file names
	char *namt = Form("theta_{c} lens %d angle %d diff %1.2f",
					  lens,(int)angle,tdiff);
	char *pict = Form("%s/reco_%s_%d_%d_%1.2f.png",
					  savepath,runtype,lens,(int)angle,tdiff);

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

//.........这里部分代码省略.........
     ttbarX->SetFillColor(kGray);
     Stack->Add(ttbarX);
   }
   if(WZ_Sherpa){
     WZ_Sherpa->SetLineColor(kBlue-9);
     WZ_Sherpa->SetFillColor(kBlue-9);
     Stack->Add(WZ_Sherpa);
   }
   if(otherBkgs){
     otherBkgs->SetLineColor(kOrange-2);
     otherBkgs->SetFillColor(kOrange-2);
     Stack->Add(otherBkgs);
   }

   TH1D* backgroundsForRatio = (TH1D*)allBackgrounds->Clone("bkgds_for_ratio");

   // create main pad
   const double mainPad_ylow  = 0.3;
   const double mainPad_yhigh = 0.95;
   const double mainPad_xlow = 0.;
   const double mainPad_xhigh = 0.95;
   const double pad_margin_left = 0.2;
   const double pad_margin_right  = 0.02;
   double main_y_max = -99;
   double main_y_min = -99;

   TPad* mainPad = new TPad("main","main",mainPad_xlow,mainPad_ylow,mainPad_xhigh,mainPad_yhigh);
   mainPad->SetMargin(pad_margin_left,pad_margin_right,0.,.05); // left, right, bottom, top
   mainPad->Draw();
   mainPad->cd();

   // find max y
   if(allBackgrounds){
     int maxBinBkg = allBackgrounds->GetMaximumBin();
     double bkgYmax = allBackgrounds->GetBinContent(maxBinBkg);
     main_y_max = bkgYmax;
   }
   if(data){
     int maxBinData = data->GetMaximumBin();
     double dataYmax = data->GetBinContent(maxBinData);
     if(dataYmax > main_y_max) main_y_max = dataYmax;
   }
   if(main_y_max < .002) main_y_max = .1;

   if(doLogAxis){
     mainPad->SetLogy();
     main_y_min = 0.07;
     main_y_max*=60;
   }else{
     main_y_min = 0.001;
     main_y_max*=1.6;
   }

   // draw axis
   allBackgrounds->SetMaximum(main_y_max);
   allBackgrounds->SetMinimum(main_y_min);
   std::string title = std::string(";") + axisName + ";Events";
   std::stringstream binw; 
   
   binw<<allBackgrounds->GetBinWidth(1);
   std::string width = binw.str();
   if (axisName.find("GeV") != std::string::npos){
     width += " GeV";
     title += " / "+width;
   }
   allBackgrounds->SetTitle(title.c_str());