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

int RampCorrellation_single(std::string file_name, TGraphErrors &graph, TGraphErrors &graph2, vector<double> &vec_errors)
{

	TTree *t = new TTree();
	const std::string file_path = "/home/jlab/github/MagCloak_Analysis/calibration/data-calib/DATA_MegaVIEW/";
//	std::string file_name = "DataFile_2016-12-08_06-59-11.csv";
	std::string file = file_path+file_name;

	t->ReadFile(file.c_str());
//	t->Print();

	TCanvas *c = new TCanvas();
	t->Draw("TMath::Abs(B1/B2):time");
	TGraph *gh = (TGraph*)c->GetListOfPrimitives()->FindObject("Graph");
	double ratio_mean = gh->GetMean(2);
	double ratio_std = gh->GetRMS(2);

	t->Draw("Bnom:time");
	TGraph *gh1 = (TGraph*)c->GetListOfPrimitives()->FindObject("Graph");
	double nom_mean = gh1->GetMean(2);
	double nom_std = 0;

//	cout << "At " << nom_mean << " mT, B1/B2 is: " << ratio_mean << " +/- " << ratio_std << endl;
	
	c->Close();

	int n = graph.GetN();
	graph.SetPoint(n,nom_mean,ratio_mean);
	graph.SetPointError(n,nom_std,ratio_std);

//-------------------------------------------------------------------------------

	TCanvas *c2 = new TCanvas();
	t->Draw("B3:time");
	TGraph *gh2 = (TGraph*)c2->GetListOfPrimitives()->FindObject("Graph");
	double B3_mean = -1*(gh2->GetMean(2));
	double B3_std = gh2->GetRMS(2);
	
	int n2 = graph2.GetN();
	graph2.SetPoint(n2, nom_mean, B3_mean);
	graph2.SetPointError(n2, nom_std, B3_std);
	
	c2->Close();

	if(nom_mean < 500) vec_errors.push_back(B3_std);

	return 0;
}

void CalcQCDNormFactor() {
  //TFile *f = TFile::Open("results/Plotter_out_2016_05_29_22h19m32.root"); // 76X Silver JSON
  TFile *f = TFile::Open("results/Plotter_out_2016_06_21_15h27m59.root"); // 76X Golden JSON
  TCanvas *c = (TCanvas*)f->Get("NJet/PlotSamples_Pass5Cuts_PassHLT");
  THStack *s = (THStack*)c->GetListOfPrimitives()->At(1);
  TH1D *data = (TH1D*)c->GetListOfPrimitives()->At(3);
  double MC_integral = 0;
  double QCD_count = 0;
  for (int i=s->GetHists()->GetEntries()-1; i>=0; --i) {
    TH1D* h = (TH1D*)s->GetHists()->At(i);
    if (i==s->GetHists()->GetEntries()-1) QCD_count = h->Integral();
    std::cout<<h->GetName()<<" "<<h->Integral()<<std::endl;
    MC_integral += h->Integral();
  }
  double data_QCD_estimate = data->Integral()- (MC_integral-QCD_count);
  double QCD_scale = data_QCD_estimate/QCD_count;
  std::cout<<"MC:                  "<<MC_integral<<std::endl;
  std::cout<<"Data:                "<<data->Integral()<<std::endl;
  std::cout<<"MC   (QCD only):     "<<QCD_count<<std::endl;
  std::cout<<"Data (QCD only est): "<<data_QCD_estimate<<std::endl;
  std::cout<<"QCD Scale: "<<QCD_scale<<std::endl;

  TH1D* qcd = (TH1D*)s->GetHists()->At(s->GetHists()->GetEntries()-1);
  qcd->Scale(QCD_scale);
  c->Draw();
  
  //TCanvas *c = (TCanvas*)f->Get("NJet/PlotSamples_Pass5Cuts_PassHLT_Ratio");
  //
  //TH1D* ratio = (TH1D*)((TVirtualPad*)c->cd(2))->GetListOfPrimitives()->At(0);
  //TF1* fit = new TF1("fit","pol1", 400, 2000);
  //ratio->Fit("fit","RBQ0");
  //fit->SetLineColor(2);
  //fit->SetLineWidth(1);
  //TF1* fit_up   = (TF1*)fit->Clone("fit_up");
  //TF1* fit_down = (TF1*)fit->Clone("fit_down");
  //fit_up  ->SetParameter(0,fit->GetParameter(0)+fit->GetParError(0)*1);
  //fit_down->SetParameter(0,fit->GetParameter(0)-fit->GetParError(0)*1);
  //fit_up  ->SetParameter(1,fit->GetParameter(1)+fit->GetParError(1)*1);
  //fit_down->SetParameter(1,fit->GetParameter(1)-fit->GetParError(1)*1);
  //fit_up  ->SetLineColor(4); fit_up  ->Draw("SAME");
  //fit_down->SetLineColor(4); fit_down->Draw("SAME");
  //fit->Draw("SAME");
  //
  //std::cout<<"Fit result:"<<std::endl;
  //std::cout<<"p0: "<<fit->GetParameter(0)<<" +- "<<fit->GetParError(0)*1<<std::endl;
  //std::cout<<"p1: "<<fit->GetParameter(1)<<" +- "<<fit->GetParError(1)*1<<std::endl;
  //f->Close();  
}

TH1D* getSimplePlot(TString INPUTDIR_PREFIX, TString SCENARIO, TString dataMC, TString vartype, TString SCEN_TRIG, TString RR) {

   TString DENOM = "_Glb_pass_&_Tight2012_pass"; //"_&_Glb_pass_&_Tight2012_pass";
   if (SCENARIO == "Glb_Tight2012") DENOM = "";
   else if (SCENARIO == "Glb_Tight2012_IsolPFRelCombNoEGammaR03PU_"+SCEN_TRIG) DENOM = "_Glb_pass_&_IsolPFRelCombNoEGammaR03PU_pass_&_Tight2012_pass";
   else if (SCENARIO == "Glb_Tight2012_IsolPFRelCombNoEGammaR03_"+SCEN_TRIG) DENOM = "_Glb_pass_&_IsolPFRelCombNoEGammaR03_pass_&_Tight2012_pass";

   TString POSTFIX = "";
   if (vartype == "vtx") POSTFIX = ""; //"_vtx";
   else if (vartype == "run") POSTFIX = ""; //"_rrr";
   else if (vartype == "rrr2") POSTFIX = "_rrr2";
   else if (vartype == "rrr3") POSTFIX = "_rrr3";

   if (dataMC == "datalike_mc") {
     RR = "";
   } else {
     RR = "_"+RR;
   }
   TFile *thisf = new TFile(INPUTDIR_PREFIX+"TnP_2011_MuonID_item_"+dataMC+"_"+SCENARIO+POSTFIX+"_"+vartype+RR+".root");
   cout << "HERE " << INPUTDIR_PREFIX+"TnP_2011_MuonID_item_"+dataMC+"_"+SCENARIO+POSTFIX+"_"+vartype+RR+".root" << endl;
   thisf->cd();
   gDirectory->cd("tpTree/"+SCENARIO+"_"+vartype+"/fit_eff_plots");
   cout << "tpTree/"+SCENARIO+"_"+vartype+"/fit_eff_plots" << endl;

   TCanvas* c = (TCanvas*)gDirectory->Get(getFolder(vartype)+DENOM);
   cout << getFolder(vartype)+DENOM << endl;
   c->GetListOfPrimitives();
   RooHist* hpt = (RooHist*)c->GetPrimitive("hxy_fit_eff");

   TH1D* test = rooHist_to_TH1D(hpt,vartype,SCENARIO);
   for (int ibin = 0; ibin < test->GetXaxis()->GetNbins();ibin++) {
       cout << "AFTER " << test->GetBinContent(ibin+1) << endl;
    }
   return rooHist_to_TH1D(hpt,vartype,SCENARIO);
}

  double AutoSetYRange(TCanvas& canv, double maxScale) {
    TList* list      = canv.GetListOfPrimitives();
    double maximum   = 0;
    int    firstHist = -1;
    //    int isCanvasLogY = canv.GetLogy();
    for (int iPrims = 0; iPrims <= list->LastIndex(); ++iPrims) {
      TH1* hist = dynamic_cast<TH1*>(list->At(iPrims));
      if (hist) {
        //Remember histo to set maximum of, which is the first one drawn
        if (firstHist == -1) {
          firstHist = iPrims;
        }
        if (hist->GetMaximum() > maximum) {
          maximum = hist->GetMaximum();
        }
      }
    }

    if (firstHist != -1) {
      dynamic_cast<TH1*>(list->At(firstHist))->SetMaximum(maximum * maxScale);
      return maximum * maxScale;
    } else {
      std::cout << __func__ << " No Histograms found" << std::endl;
      return -1;
    }
  }

//----------  Retrieve data histo  -----------------
TH1F* GetData(TFile* fin, string region, string varname)
{
    string cname = CHANNEL_NAME+string("/")+region+"/"+varname;
    TCanvas* c = (TCanvas*) fin->Get(cname.c_str());
    TList* l = c->GetListOfPrimitives();
    TPad* pad = (TPad*) l->At(0);
    string hname = "v:"+varname+"|r:"+region+string("|c:")+CHANNEL_NAME+string("|t:1DSumData");
    TH1F* h = (TH1F*) pad->GetPrimitive(hname.c_str());
    return (TH1F*) h->Clone();
}

TH1F* GetHisto(TFile* fin, string region, string process, string varname, float& norm, bool do_norm, float input_norm)
{
    string cname = CHANNEL_NAME+string("/")+region+"/"+varname;
    TCanvas* c = (TCanvas*) fin->Get(cname.c_str());
    string hname = "v:"+varname+"|p:"+process+"|r:"+region+string("|c:")+CHANNEL_NAME+string("|t:1DEntries");
    TH1F* h = 0;
    if(VERBOSE>0){
 	 cerr<<"cname :"<<cname<<endl;
   	 cerr<<"histo name: "<<hname<<endl;
   	 cerr<<"pointer: "<<c<<endl;
    } 
    TList* l = c->GetListOfPrimitives();
    TPad* pad = (TPad*) l->At(0);
    THStack* stack = (THStack*) pad->GetPrimitive("");
    h = (TH1F*) stack->GetHists()->FindObject(hname.c_str());
    if(do_norm) h->Scale(input_norm/h->Integral());
    norm = h->Integral();
    return (TH1F*) h->Clone();
}

OptRes Optimization(string var = "topness", string signal = "T2tt_850_100", string region = "baseline", string channel = "allChannels", string filename = "plots/plotsProducer/1DStack.root"){

	OptRes res;
	
	//TFile* fplots = new TFile("plots/plotsProducer/1DStack.root","READ");
	TFile* fplots = new TFile(filename.c_str(),"READ");
	//string channel = "allChannels";
	//string region = "baseline";
	//string var = "topness";
	//string var = "topness_m5";
	//string signal = "T2tt_850_100";
	
	
	//retrieve canvas
	string cname = channel+"/"+region+"/"+var;
	TCanvas* c = fplots->Get(cname.c_str());
	//c->Draw();
	
	//signal
	string splot_name = "v:"+var+"|p:"+signal+"|r:"+region+"|c:"+channel+"|t:1DEntries";
	TH1F* h_sig = (TH1F*) c->GetPrimitive(splot_name.c_str());
	//sig->Draw();
	
	//////////////////////////////////////////////
	// Search the Stack that contains all bkg
	//////////////////////////////////////////////
	vector<TH1F*> h_bkgs;
	//c->GetListOfPrimitives()->Print();
	TIter next(c->GetListOfPrimitives());
	TObject* h;
	while ((h = (TH1F*) next())){
	//obj->Draw(next.GetOption());
		if(h->InheritsFrom("THStack")){
			THStack* stack = (THStack*) h;
			for(int i=0;i<stack->GetStack()->GetEntries();i++){
				TH1F* hist = (TH1F*) stack->GetStack()->At(i);
				//cout<<hist->IsA()->GetName()<<endl;
				//cout<<hist->GetName()<<" "<<hist->Integral()<<endl;
				//search tt_2l
				//if(string(hist->GetName()).find("W+jets")!=string::npos) continue;
				if(string(hist->GetName()).find("ttbar_2l")!=string::npos){
					//cout<<"mean = "<<hist->GetMean()<<" rms = "<<hist->GetRMS()<<" kurtosis = "<<hist->GetKurtosis()<<" skewness = "<<hist->GetSkewness()<<endl;
				}
				//else continue;
				//conly save the last one because plots are cumulative !!!
				// that's cheat ....
				if(i==stack->GetStack()->GetEntries()-1) h_bkgs.push_back(hist);
			}
			//cout<<((THStack*)h)->GetStack()->GetEntries()<<endl;
			//TH1* h_bkg = ((TH1*)(h->GetStack()->Last()));
		}
	}
	
	
	////////////////////////////////////////////////
	//     Loop over all bins and compute eff/sig
	////////////////////////////////////////////////
	
	res.h_eff_sig = (TH1F*) h_sig->Clone();	
	res.h_eff_bkg = (TH1F*) h_sig->Clone();	
	res.h_Zbi = (TH1F*) h_sig->Clone();	
	TH1F* h_ROC = new TH1F("h_ROC","ROC curve",10,0,1);

	//Compute the integral
	float integ_sig = res.h_eff_sig->Integral(0,res.h_eff_sig->GetNbinsX()+1);
	float integ_bkg = 0;
	for(unsigned ih=0;ih<h_bkgs.size();ih++){
		integ_bkg+=h_bkgs[ih]->Integral(0,h_bkgs[ih]->GetNbinsX()+1);
	}
	///////////////////////

	//-- efficiencoy for signal
	for(int i=1;i<res.h_eff_sig->GetNbinsX()+1;i++){
		float eff = 0;
		if(integ_sig!=0)  eff = res.h_eff_sig->Integral(i,res.h_eff_sig->GetNbinsX()+1)/integ_sig;
		res.h_eff_sig->SetBinContent(i,eff);
		res.h_eff_sig->SetBinError(i,0);
	}
	//-- efficiency for bkg
	for(int i=1;i<res.h_eff_sig->GetNbinsX()+1;i++){
		float eff = 0;
		for(unsigned ih=0;ih<h_bkgs.size();ih++){
			eff += h_bkgs[ih]->Integral(i,h_bkgs[ih]->GetNbinsX()+1);
		}
		if(integ_bkg!=0)  eff /= integ_bkg;
		res.h_eff_bkg->SetBinContent(i,eff);
		res.h_eff_bkg->SetBinError(i,0);
	}
	//-- compute significance
	//vector<pair<float,float> > roc_v;
	float* roc_sig = new float[res.h_eff_sig->GetNbinsX()];
	float* roc_bkg = new float[res.h_eff_sig->GetNbinsX()];
	int imax = 0;
	for(int i=1;i<res.h_eff_sig->GetNbinsX()+1;i++){
		float Ns = res.h_eff_sig->GetBinContent(i)*integ_sig;
		//float Ns = res.h_eff_sig->GetBinContent(i);
		//float Nb = res.h_eff_bkg->GetBinContent(i);
		float Nb = res.h_eff_bkg->GetBinContent(i)*integ_bkg;
		float signif = 0;
		if(Nb!=0) signif = Ns/sqrt(Nb);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  Float_t legendX2 = .955; //.70;
  Float_t legendY1 = .53; //.65;
  Float_t legendY2 = .78;

  TLegend* catLeg1 = new TLegend(legendX1,legendY1,legendX2,legendY2);
  //catLeg1->SetTextSize(0.060);
  catLeg1->SetTextSize(0.044);
  catLeg1->SetTextFont(42);
  catLeg1->SetFillColor(0);
  catLeg1->SetLineColor(1);
  catLeg1->SetBorderSize(1);

  //
  // Define canvas
  //
  TCanvas *canvas = new TCanvas("canvas","canvas",10,10,W,H);

  canvas->SetFillColor(0);
  canvas->SetBorderMode(0);
  canvas->SetFrameFillStyle(0);
  canvas->SetFrameBorderMode(0);
  canvas->SetLeftMargin( L/W );
  canvas->SetRightMargin( R/W );
  canvas->SetTopMargin( T/H );
  canvas->SetBottomMargin( B/H );
  canvas->SetTickx(0);
  canvas->SetTicky(0);

  canvas->Divide(1, 2);
  
  //
  // Define pads
  //
  TPad* canvas_up = (TPad*) canvas->GetListOfPrimitives()->FindObject("canvas_1");
  TPad* canvas_dw = (TPad*) canvas->GetListOfPrimitives()->FindObject("canvas_2");

  //
  // define the size
  double up_height     = 0.8;  // please tune so that the upper figures size will meet your requirement
  double dw_correction = 1.30; // please tune so that the smaller canvas size will work in your environment
  double font_size_dw  = 0.1;  // please tune the font size parameter for bottom figure
  double dw_height     = (1. - up_height) * dw_correction;
  double dw_height_offset = 0.04; // KH, added to put the bottom one closer to the top panel

  //
  // set pad size
  canvas_up->SetPad(0., 1 - up_height,    1., 1.00);
  canvas_dw->SetPad(0., 0.,               1., dw_height+dw_height_offset);
  //
  canvas_up->SetFrameFillColor(0);
  canvas_up->SetFillColor(0);
  canvas_up->SetTopMargin(0.12);
  canvas_up->SetLeftMargin(0.1);
  //
  canvas_dw->SetFillColor(0);
  canvas_dw->SetFrameFillColor(0);
  canvas_dw->SetBottomMargin(0.35);
  canvas_dw->SetTopMargin(0);
  canvas_dw->SetLeftMargin(0.1);
  
  //
  // draw top figure
  canvas_up->cd();

  TH1D * GenHist, * EstHist,* thist;
  TH1D * GenHistTemp, * EstHistTemp;

void recupThePlots(){
//	TDirectory *theDr = (TDirectory*) myFile->Get("eleIDdir");///denom_pt/fit_eff_plots");
	//theDr->ls();
	cout << "coucou" << theOutFileName << endl;
	
	TFile *myOutFile = new TFile(theOutFileName,"RECREATE");
	
	TSystemDirectory dir(thePath, thePath);
	TSystemFile *file;
	TString fname;
	TIter next(dir.GetListOfFiles());
	while (((file=(TSystemFile*)next()))) {
		fname = file->GetName();
		if ((fname.BeginsWith("TnP"))&&fname.Contains("data")) {
                        cout <<  "--------------------"<< "\n";
                        cout << fname << "\n";
			TFile *myFile = new TFile(fname);
			TIter nextkey(myFile->GetListOfKeys());                   
			TKey *key;
			while ((key = (TKey*)nextkey())) {
                                // cout << key << "\n";
				TString theTypeClasse = key->GetClassName();
				TString theNomClasse = key->GetTitle();
                                   cout << "theTypeClasse:  "<< theTypeClasse << " ,  " << theNomClasse << "\n";
				if ( theTypeClasse == "TDirectoryFile" ){
                                        //    cout << "we are here 1" << "\n";
					TDirectory *theDr = (TDirectory*) myFile->Get(theNomClasse);
					TIter nextkey2(theDr->GetListOfKeys());
					TKey *key2;
					while ((key2 = (TKey*)nextkey2())) {
						TString theTypeClasse2 = key2->GetClassName();
						TString theNomClasse2 = key2->GetTitle();	
						if ( theTypeClasse == "TDirectoryFile" || theTypeClasse == "TGraphAsymmErrors" ){
							TDirectory *theDr2 = (TDirectory*) myFile->Get(theNomClasse+"/"+theNomClasse2+"/fit_eff_plots");
							TIter nextkey3(theDr2->GetListOfKeys());
							TKey *key3;
							while ((key3 = (TKey*)nextkey3())) {
								TString theTypeClasse3 = key3->GetClassName();
								TString theNomClasse3 = key3->GetName();	
				 				cout << "type = " << theTypeClasse3 << " nom = " << theNomClasse3 << endl;
						
								TCanvas *theCanvas = (TCanvas*) myFile->Get(theNomClasse+"/"+theNomClasse2+"/fit_eff_plots/"+theNomClasse3);
								TIter nextObject(theCanvas->GetListOfPrimitives());
								TObject *obj;
								while ((obj = (TObject*)nextObject())) {
									if (obj->InheritsFrom("TGraphAsymmErrors")) {
					        				cout << "histo: " << obj->GetName() << endl;
										myOutFile->cd();
										obj->Write(theNomClasse2+"_"+theNomClasse3);
										myFile->cd();
									}
									if (obj->InheritsFrom("TH2F")) {
										cout << "the TH2F = " << obj->GetName() << endl;
										myOutFile->cd();
										obj->Write(theNomClasse2+"_"+theNomClasse3);
										myFile->cd();
									}
								}
							}
						}

					}
			
				}
			}
			delete myFile;
		}
	
	}
	myOutFile->Close();

}

void drawSingleFit( DrawBase* db, const std::string& varName, int ptMin, int ptMax, int rhoBin ) {

    TFile* fitFile = db->get_mcFile(0).file;

    char canvasName[400];
    sprintf( canvasName, "%s_pt_%d_%d_rho%d", varName.c_str(), ptMin, ptMax, rhoBin );

    char gluonHistoName[400];
    sprintf( gluonHistoName, "%s_gluon_pt%d_%d_rho%d", varName.c_str(), ptMin, ptMax, rhoBin );
    char quarkHistoName[400];
    sprintf( quarkHistoName, "%s_quark_pt%d_%d_rho%d", varName.c_str(), ptMin, ptMax, rhoBin );

    TCanvas* thisCanvas = (TCanvas*)fitFile->Get(canvasName);
    thisCanvas->cd();


    TList* list = thisCanvas->GetListOfPrimitives();

    TH1D* h1_gluon = (TH1D*)list->FindObject(gluonHistoName);
    TH1D* h1_quark = (TH1D*)list->FindObject(quarkHistoName);

    int quarkColor = 46;
    int gluonColor = 38;

    h1_gluon->SetMarkerStyle(21);
    h1_gluon->SetMarkerColor(gluonColor);
    h1_gluon->SetLineColor(gluonColor);
    h1_gluon->SetLineWidth(2);
    h1_gluon->SetFillColor(gluonColor);
    h1_gluon->SetFillStyle(3004);

    h1_quark->SetMarkerStyle(20);
    h1_quark->SetLineColor(quarkColor);
    h1_quark->SetMarkerColor(quarkColor);
    h1_quark->SetLineWidth(2);
    h1_quark->SetFillColor(quarkColor);
    h1_quark->SetFillStyle(3005);


    std::string functionName = (varName=="ptD") ? "functionPtD" : "gamma";
    std::string functionName_gluon = functionName + "_gluon";
    std::string functionName_quark = functionName + "_quark";

    TF1* f1_gluon = (TF1*)list->FindObject(functionName_gluon.c_str());
    TF1* f1_quark = (TF1*)list->FindObject(functionName_quark.c_str());

    f1_gluon->SetLineColor(gluonColor);
    f1_quark->SetLineColor(quarkColor);


    std::string axisName;
    if( varName=="nCharged" )
        axisName = "Jet Charged Multiplicity";
    else if( varName=="nNeutral" )
        axisName = "Jet Neutral Multiplicity";
    else if( varName=="ptD" )
        axisName = "Jet p_{T}D";


    float xMin = h1_quark->GetXaxis()->GetXmin();
    float xMax = (varName=="ptD") ? h1_quark->GetXaxis()->GetXmax() : 50.;
    float yMax_gluon = h1_gluon->GetMaximum();
    float yMax_quark = h1_quark->GetMaximum();
    float yMax = (yMax_gluon>yMax_quark) ? yMax_gluon : yMax_quark;
    yMax *= 1.3;

    std::string yTitle = "Normalized to Unity";

    TH2D* h2_axes = new TH2D("axes", "", 10, xMin, xMax, 10, 0., yMax);
    h2_axes->SetXTitle(axisName.c_str());
    h2_axes->SetYTitle(yTitle.c_str());

    TPaveText* label_cms = db->get_labelCMS();
    TPaveText* label_sqrt = db->get_labelSqrt();

    int rhoMin = rhoBin;
    int rhoMax = rhoBin+1;

    char ptText[400];
    sprintf( ptText, "%d < p_{T} < %d GeV", ptMin, ptMax );
    char rhoText[400];
    sprintf( rhoText, "%d < #rho_{PF} < %d GeV", rhoMin, rhoMax );

    TPaveText* pt_rho_label = new TPaveText( 0.25, 0.79, 0.45, 0.91, "brNDC" );
    pt_rho_label->SetTextSize(0.035);
    pt_rho_label->SetFillColor(0);
    pt_rho_label->AddText( ptText );
    pt_rho_label->AddText( rhoText );

    char legendTitle[500];
    sprintf( legendTitle, "%d < p_{T} < %d, %d < #rho < %d", ptMin, ptMax, rhoMin, rhoMax );

    //TLegend* legend = new TLegend( 0.55, 0.65, 0.9, 0.9, legendTitle );
    TLegend* legend = new TLegend( 0.65, 0.75, 0.9, 0.9 );
    legend->SetFillColor(0);
    legend->SetTextSize(0.035);
    legend->AddEntry( h1_gluon, "Gluons", "PL" );
    legend->AddEntry( h1_quark, "Quarks", "PL" );

    TCanvas* c1 = new TCanvas("c1", "", 600, 600);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  }
  TFile *f = TFile::Open(filename.c_str());
  for (size_t iMethod = 0; iMethod<4; ++iMethod) if (!(iMethod==0&&samples[0].size()==0)){
    // Print Top row for each method
    if (latex) {
      if (iMethod==0) {
        printf("\\hline\n");
        printf("Method 2 & A & B & C & D = B*C/A & D obs. & Ratio pred./obs. & Pull & KS test\\\\\n");
      }	else if (iMethod==1){
        printf("\\hline\n");
        printf("Method 1 & A & B & C & D = B*C/A & D obs. & Ratio pred./obs. & Pull & KS test\\\\\n");
      }
      printf("\\hline\n");
    } else {
      std::stringstream r_sb_cut;
      if (R_CUT_LOW==0) r_sb_cut<<"R<"<<R_CUT;
      else r_sb_cut<<R_CUT_LOW<<"<R<"<<R_CUT;
      const char* prime = ABCD_prime ? "'" : "";
      //if (iMethod==0) printf("| *Sample* | *A (DPhi>2.8, SB)* | *B (DPhi<2.8, SB)* | *C (DPhi>2.8, Sig.B.)* | - | *D = B*C/A pred.* | *D (DPhi<2.8, Sig.B.) obs.* | *Ratio pred./obs.* |\n");
      //else if (iMethod==1) printf("| *Sample* | *A (R<%1.1f, SB)* | *B (R>%1.1f, SB)* | *C (R<%1.1f, Sig.B.)* | *D = B (R fit, SB) * C/A pred.* | *D = B*C/A pred.* | *D (R>%1.1f, Sig.B.) obs.* | *Ratio pred./obs.* | \n", R_CUT, R_CUT, R_CUT, R_CUT);
      if (iMethod==0) printf("| *Sample* | *A (DPhi>%1.1f, %s)* | *B (DPhi<%1.1f, %s)* | *C (DPhi>%1.1f, R>0.4)* | *D = B*C/A pred.* | *D (DPhi<%1.1f, R>0.4) obs.* | *Ratio pred./obs.* | *Pull (pred-obs)/error* | *KS test* |\n", DPHI_CUT, r_sb_cut.str().c_str(), DPHI_CUT, r_sb_cut.str().c_str(), DPHI_CUT, DPHI_CUT);
      else if (iMethod==1) printf("| *Sample* | *A%s (%s, <2 tag)* | *B%s (R>%1.1f, <2 tag)* | *C%s (%s, 2 tag)* | *D%s = B%s*C%s/A%s pred.* | *D%s (R>%1.1f, 2 tag) obs.* | *Ratio pred./obs.* | *Pull (pred-obs)/error* | *KS test* |\n", prime, r_sb_cut.str().c_str(), prime, R_CUT, prime, r_sb_cut.str().c_str(), prime, prime, prime, prime, prime, R_CUT);
    }
    TH1D *h_side_sum, *h_signal_sum;
    for (size_t iSample = 0; iSample<samples[iMethod].size(); ++iSample)  {
      std::string canname = 
	iMethod==0 ? std::string("DPhiBins")+(ABCD_prime ? "/RBins_0To1HadTop_" : "/RBins_2HadTop_")+samples[iMethod][iSample] :
	iMethod==1 ? std::string("RFine/Tau32Cuts_")+(ABCD_prime ? "Fail" : "Pass")+"DPhiCut_"+samples[iMethod][iSample] :
	iMethod==2 ? std::string("RFine/Tau32Cuts_")+(ABCD_prime ? "Fail" : "Pass")+"DPhiCut_Background" :
	iMethod==3 ? std::string("RFine/Tau32Cuts_")+(ABCD_prime ? "Fail" : "Pass")+"DPhiCut_"+samples[iMethod][iSample] : "";
      TCanvas *can = (TCanvas*)(f->Get(canname.c_str()));
      can = (TCanvas*)can->Clone();
      can->Draw();
      TH1D *h_side = (TH1D*)can->GetListOfPrimitives()->At(i_h_side[iMethod]);
      TH1D *h_signal = (TH1D*)can->GetListOfPrimitives()->At(i_h_signal[iMethod]);
      // Simulate different cross section by scaling a certain background
      if (iMethod==1) {
        TH1D *h_side_temp_scaled = (TH1D*)h_side->Clone(); h_side_temp_scaled->Scale(scale_factors[iSample]);
        TH1D *h_signal_temp_scaled = (TH1D*)h_signal->Clone(); h_signal_temp_scaled->Scale(scale_factors[iSample]);
        if (iSample==0) {
          h_side_sum = h_side_temp_scaled;
          h_signal_sum = h_signal_temp_scaled;
        } else {
          h_side_sum->Add(h_side_temp_scaled);
          h_signal_sum->Add(h_signal_temp_scaled);
        }
      } else if (iMethod==2) {
        h_side = h_side_sum;
        h_signal = h_signal_sum;
      }
      TH1D *h_pred =(TH1D*)h_side->Clone();
      if (iMethod!=0&&rebin>1) { h_side->Rebin(rebin); h_signal->Rebin(rebin); h_pred->Rebin(rebin); }
      TLegend *leg = (TLegend*)can->GetListOfPrimitives()->At(can->GetListOfPrimitives()->GetEntries()-1);
      leg->SetX1(0.35); leg->SetX2(0.65); leg->SetY1(0.6);
      
      // Add ratio plot
      int y1 = 350;
      int y2 = 150;
      int mid2 = 10;
      can->Divide(1,2);
      // Pad 1 (80+500+20 x 40+500)
      TVirtualPad* p = can->cd(1);
      p->SetPad(0,(y2+60+mid2)/(y1+y2+100.0+mid2),1,1);
      p->SetTopMargin(40.0/(y1+40));
      p->SetBottomMargin(0);
      p->SetRightMargin(0.05);

Plot_searchBin_Data(string sample="stacked",string histname="searchH",int choice=1){

  ///////////////////////////////////////////////////////////////////////////////////////////
  ////Some cosmetic work for official documents.
  //
  // Set basic style
  //
  gROOT->LoadMacro("tdrstyle.C");
  setTDRStyle();
  gStyle->SetPalette(1) ; // for better color output
  gROOT->LoadMacro("CMS_lumi_v2.C");

  int W = 600;
  int H = 600;
  int H_ref = 600;
  int W_ref = 800;
  float T = 0.08*H_ref;
  float B = 0.12*H_ref;
  float L = 0.12*W_ref;
  float R = 0.08*W_ref;

  ///////////////////////////////////////////////////////////////////////////////////////////
  //
  // More specific style set, opening input files etc

  gStyle->SetOptStat(0);  ///to avoid the stat. on the plots
  char tempname[200];
  char xtitlename[200];
  char ytitlename[200];

  if(sample.find("stack")==string::npos)sprintf(tempname,"TauHad/GenInfo_HadTauEstimation_%s.root",sample.c_str());
  else sprintf(tempname,"TauHad/Stack/GenInfo_HadTauEstimation_%s.root",sample.c_str());
//cout << "warning:\n Warning \n \n  using elog195 for pre and  exp \n \n ";
  TFile * GenFile = new TFile(tempname,"R");
  sprintf(tempname,"TauHad2/Elog378_HadTauEstimation_data_SingleMuon_v15cd_.root");
  TFile * EstFile = new TFile(tempname,"R");

  //
  // Define legend
  //
  Float_t legendX1 = .55; //.50;
  Float_t legendX2 = .85; //.70;
  Float_t legendY1 = .75; //.65;
  Float_t legendY2 = .90;

  TLegend* catLeg1 = new TLegend(legendX1,legendY1,legendX2,legendY2);
  catLeg1->SetTextSize(0.032);
  catLeg1->SetTextFont(42);

  catLeg1->SetTextSize(0.042);
  catLeg1->SetTextFont(42);
  catLeg1->SetFillColor(0);
  catLeg1->SetLineColor(0);
  catLeg1->SetBorderSize(0);

  //
  // Define canvas
  //
  TCanvas *canvas = new TCanvas("canvas","canvas",10,10,W,H);

  canvas->SetFillColor(0);
  canvas->SetBorderMode(0);
  canvas->SetFrameFillStyle(0);
  canvas->SetFrameBorderMode(0);
  canvas->SetLeftMargin( L/W );
  canvas->SetRightMargin( R/W );
  canvas->SetRightMargin( 0.2 );
  canvas->SetTopMargin( T/H );
  canvas->SetBottomMargin( B/H );
  canvas->SetTickx(0);
  canvas->SetTicky(0);

  canvas->Divide(1, 2);

  //
  // Define pads
  //
  TPad* canvas_up = (TPad*) canvas->GetListOfPrimitives()->FindObject("canvas_1");
  TPad* canvas_dw = (TPad*) canvas->GetListOfPrimitives()->FindObject("canvas_2");
 
  // define the size
  double up_height     = 0.8;  // please tune so that the upper figures size will meet your requirement
  double dw_correction = 1.30; // please tune so that the smaller canvas size will work in your environment
  double font_size_dw  = 0.1;  // please tune the font size parameter for bottom figure
  double dw_height    = (1. - up_height) * dw_correction;
  double dw_height_offset = 0.040; // KH, added to put the bottom one closer to the top panel
 
  // set pad size
  canvas_up->SetPad(0., 1 - up_height, 0.97, 1.);
  canvas_dw->SetPad(0., dw_height_offset, 0.97, dw_height+dw_height_offset);
  canvas_up->SetFrameFillColor(0);
  canvas_up->SetFillColor(0);
  canvas_dw->SetFillColor(0);
  canvas_dw->SetFrameFillColor(0);
  canvas_dw->SetBottomMargin(0.25);
  
  // set top margin 0 for bottom figure
  canvas_dw->SetTopMargin(0);
  
  // draw top figure
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	legDta->AddEntry(gr,f->label,"l");
      }
    }
    f++;
  }
    

  mgrDta->Draw("p");
  //  legDta->Draw();   don't draw the data legend
  
  TMultiGraph* mgrThry = new TMultiGraph("Theory","G_{E}^{n}");
  TLegend *legThry = new TLegend(.54,.6,.875,.9,"","brNDC");

  wgr = mgrThry;
  wlg = legThry;

  // the theory
  wlg->SetBorderSize(0); // turn off border
  wlg->SetFillStyle(0);
  
  f = theoryfiles1;
  gr=0;
  while ( f && f->filename ) {
    gr=OneGraph(f);
    if (gr) {
      TGraphAsymmErrors *egr = dynamic_cast<TGraphAsymmErrors*>(gr);
      if (egr && egr->GetN()>1 && egr->GetEYhigh() && egr->GetEYhigh()[1]>0) {
	gr = toerror_band(egr);
	gr->SetFillStyle(3000+f->style);
      }
      if (f->lnpt) {
	wgr->Add(gr,f->lnpt);
	wlg->AddEntry(gr,f->label,f->lnpt);
      }
      else if (gr->GetMarkerStyle()>=20) {
	wgr->Add(gr,"p");
	wlg->AddEntry(gr,f->label,"p");
      }	
      else {
	wgr->Add(gr,"l");
	wlg->AddEntry(gr,f->label,"l");
      }
    }
    f++;
  }

  genf->Draw("same");
  mgrThry->Draw("c");
  legThry->AddEntry(genf,"F_{2}/F_{1} #propto ln^{2}(Q^{2}/#Lambda^{2})/Q^{2}","l");
  legThry->Draw();

  // draw a line at 1
  cn->Modified();

  cn->Update();
  cn->SaveAs(Form("%s.eps",psfile));
  cn->SaveAs(Form("%s.root",psfile));
  gSystem->Exec(Form("./replace_symbols.pl %s.eps",psfile));
  // now an overlay, hopefully matching dimensions

  // remove everything but the graph
  cn->Update();
  TList *clist = cn->GetListOfPrimitives();
  TFrame* frame = cn->GetFrame();
  for (int i=0; i<clist->GetSize(); ) {
    if (clist->At(i) != frame) {
      clist->RemoveAt(i);
    } else i++;
  }
  // draw markers in the corners
  TMarker *mkr = new TMarker(frame->GetX1(),frame->GetY1(),2);
  mkr->Draw();
  mkr = new TMarker(frame->GetX2(),frame->GetY1(),2);
  mkr->Draw();
  mkr = new TMarker(frame->GetX1(),frame->GetY2(),2);
  mkr->Draw();
  mkr = new TMarker(frame->GetX2(),frame->GetY2(),2);
  mkr->Draw();
  frame->SetLineColor(10);
  cn->Update();

  datafile_t miller = { "figure_input/Miller/lattice.GEn.rtf","Miller",
			"[0]","[1]","[1]-[3]","[2]-[1]",0,0,1,3,"F" };

  gr = OneGraph(&miller);
  TGraphAsymmErrors* egr = dynamic_cast<TGraphAsymmErrors*>(gr);
  if (egr && egr->GetEYhigh() && egr->GetEYhigh()[egr->GetN()/2]>0) {
    gr = toerror_band(egr);
    gr->SetLineStyle(1);
    gr->SetFillColor(gr->GetLineColor());
    gr->SetFillStyle(3000+miller.style);
  }
  
  gr->Draw("F");

  cn->Update();
  cn->SaveAs("gen_Miller_Overlay.eps");
  cn->SaveAs("gen_Miller_Overlay.root");
  
}

//.........这里部分代码省略.........
	      }
      }
    }
    f++;
  }

  mgrThry->Draw("c");
#ifdef PQCD
  ff->Draw("same");
  legThry->AddEntry(ff,"F2/F1 #propto ln^{2}(Q^{2}/#Lambda^{2})/Q^{2}","l");
#endif// PQCD
  legThry->Draw();



  cn->Modified();
  cn->Update();
  TFrame* frame = gPad->GetFrame();

  // draw a line at 1
  TLine *line1 = new TLine(frame->GetX1(),1.,frame->GetX2(),1.);
  line1->SetLineStyle(1);
  line1->Draw();

  return;
  cn->Update();
  cn->SaveAs(Form("%s.eps",psfile));
  cn->SaveAs(Form("%s.root",psfile));
  gSystem->Exec(Form("./replace_symbols.pl %s.eps",psfile));

  // now an overlay, hopefully matching dimensions
  // remove everything except the frame
  cn->Update();
  TList *clist = cn->GetListOfPrimitives();
  frame = cn->GetFrame();
  for (int i=0; i<clist->GetSize(); ) {
	  if (clist->At(i) != frame) {
		  clist->RemoveAt(i);
	  } else i++;
  }
  // draw markers in the corners
  TMarker *mkr = new TMarker(frame->GetX1(),frame->GetY1(),2);
  mkr->Draw();
  mkr = new TMarker(frame->GetX2(),frame->GetY1(),2);
  mkr->Draw();
  mkr = new TMarker(frame->GetX1(),frame->GetY2(),2);
  mkr->Draw();
  mkr = new TMarker(frame->GetX2(),frame->GetY2(),2);
  mkr->Draw();
  frame->SetLineColor(10);


  {
	  datafile_t goeckeler = { "figure_input/Goeckeler/goeckeler.dat","Goeckeler",
		  "[0]","[1]","[1]-[2]","[3]-[1]",0,0,1,4,"F" };

	  gr = OneGraph(&goeckeler);
	  TGraphAsymmErrors *egr = dynamic_cast<TGraphAsymmErrors*>(gr);
	  if (egr && egr->GetEYhigh() && egr->GetEYhigh()[egr->GetN()/2]>0) {
		  gr = toerror_band(egr);
		  gr->SetLineStyle(1);
		  gr->SetFillColor(gr->GetLineColor());
		  gr->SetFillStyle(3000+goeckeler.style);
	  }

	  gr->Draw("F");

void slopetest(){

   string dist = "maos210";
   double nevts = 2*50000.0;
   string dim = "mt";

   double p1=0, p2=0, p3=0;
   double dx = 0;
   double eps = 0.2;
   if( dim == "mt" ){
      p1 = 172.0; p2 = 172.5, p3 = 173.0;
      //p1 = 172.5; p2 = 173, p3 = 173.5;
      //p1 = 166.0; p2 = 166.5, p3 = 167.0;
      //p1 = 178.0; p2 = 178.5, p3 = 179.0;
   }
   if( dim == "jes" ){
      p1 = 0.999; p2 = 1.000, p3 = 1.001;
   }
   if( dist == "mbl" ){
      dx = 2.8;
   }
   if( dist == "mt2_221" ){
      dx = 1.9;
   }

   TFile *f = new TFile("../results/plotsTemplates.root");
   TDirectory *d = (TDirectory*)f->Get((dim+"shape_"+dist+"_gp").c_str());

   /*
   TCanvas *ctemp = (TCanvas*)f->Get(("c_"+dist+"_gp_signal_JSF1000").c_str());
   TH1 *h = (TH1*)(ctemp->GetListOfPrimitives()->At(1));
   dx = h->GetBinWidth(1);

   delete ctemp;
   */
   dx = 10;

   TGraph *gslope = new TGraph();
   TGraph *gint = new TGraph();

   double sig = 0;
   double d2sig = 0;

   TIter nextkey(d->GetListOfKeys());
   TKey *key;
   while( (key = (TKey*)nextkey()) ){

      string name = key->GetName();
      string classname = key->GetClassName();

      TCanvas *c = (TCanvas*)d->Get(name.c_str());
      TGraph *g = (TGraph*)(c->GetListOfPrimitives()->At(1));

      double y1=0, y2=0;
      double ycent = 0;
      for(int i=0; i < g->GetN(); i++){
         double x=0, y=0;
         g->GetPoint(i, x, y);
         if( x == p1 ) y1 = y;
         if( x == p2 ) ycent = y;
         if( x == p3 ) y2 = y;
      }
      double slope = y2-y1;

      string sbin = name;
      string ntemp = "c"+dist+"_gp";
      sbin.erase(sbin.begin(), sbin.begin()+ntemp.length());
      if( dist == "mt2_221" ) sbin.erase(sbin.end(), sbin.end()+1);
      double dbin = atof(sbin.c_str());

      double integrand = slope*slope/ycent;

      gint->SetPoint(gint->GetN(), dbin, integrand);
      gslope->SetPoint(gslope->GetN(), dbin, slope);

      sig += slope*slope*eps*eps*-0.5*dx/ycent;
      d2sig += integrand*dx;

   }

   cout << "Projected S(" << eps << ") = " << sig-2.0*nevts << endl;
   cout << "Projected sigma = " << sqrt(1.0/(d2sig*nevts)) << endl;

   TCanvas *c1 = new TCanvas();
   gslope->SetLineWidth(2);
   gslope->GetXaxis()->SetTitle("(GeV)");
   gslope->GetXaxis()->SetTitleSize(0.05);
   gslope->GetXaxis()->SetTitleOffset(0.8);
   gslope->Draw("AC");

   TCanvas *c2 = new TCanvas();
   gint->SetLineWidth(2);
   gint->GetXaxis()->SetTitle("(GeV)");
   gint->GetXaxis()->SetTitleSize(0.05);
   gint->GetXaxis()->SetTitleOffset(0.8);
   gint->Draw("AC");

   return;
}