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

/**
 * Helper function to extract the histograms from the canvas file. <br>
 * Takes references to pointers in order to fill them.
 */
void getHistograms(const TString canvasName, TH1F * & histo1, TH1D * & histo2, const TString & resonance)
{
  std::cout << "canvasName = " << canvasName << std::endl;
  TFile * inputFile = new TFile("plotMassOutput.root");
  TCanvas * canvas = (TCanvas*)inputFile->Get(canvasName);
  TString resonanceNum("_1");
  if( resonance == "Upsilon3S" ) resonanceNum = "_2";
  if( resonance == "Upsilon2S" ) resonanceNum = "_3";
  if( resonance == "Upsilon" ) resonanceNum = "_4";
  if( resonance == "Psi2S" ) resonanceNum = "_5";
  if( resonance == "JPsi" ) resonanceNum = "_6";

  if( resonance == "Psis" ) resonanceNum = "_1";
  if( resonance == "Upsilons" ) resonanceNum = "_2";
  if( resonance == "LowPtResonances" ) resonanceNum = "_3";
  if( resonance == "AllResonances" ) resonanceNum = "_4";

  TPad * pad = (TPad*)canvas->GetPrimitive(canvasName+resonanceNum);
  histo1 = (TH1F*)pad->GetPrimitive("hRecBestResAllEvents_Mass");
  if( resonance == "Z" || resonance == "AllResonances" ) histo2 = (TH1D*)pad->GetPrimitive("Mass_PProf");
  else histo2 = (TH1D*)pad->GetPrimitive("Mass_fine_PProf");
  // if( resonance == "Z" || resonance == "AllResonances" ) histo2 = (TH1D*)pad->GetPrimitive("Mass_Probability");
  // else histo2 = (TH1D*)pad->GetPrimitive("Mass_fine_Probability");

  // cout << "histo1 = " << histo1 << ", histo2 = " << histo2 << endl;
  // cout << "histo1 = " << histo1->GetEntries() << ", histo2 = " << histo2->GetEntries() << endl;
}

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);
}

void analyze3MuonEvents(const string& version)
{
  gROOT->Reset();

  vector<string> MTBins;
  MTBins.push_back("_lowMT");
  MTBins.push_back("_highMT");

  vector<string> tauIso;
  tauIso.push_back("Iso");
  tauIso.push_back("NonIso");

  vector<string> muIso;
  muIso.push_back("data");
  muIso.push_back("nonIsoWData");

  vector<string> methods;
  methods.push_back("ShareTrack");
  methods.push_back("SoftMu");
  methods.push_back("SoftMu5GeV");
  methods.push_back("SoftMu15GeV");
  methods.push_back("SoftMu20GeV");

  for (vector<string>::const_iterator iMTBin = MTBins.begin(); iMTBin != MTBins.end(); ++iMTBin) {

    for (vector<string>::const_iterator iTauIso = tauIso.begin(); iTauIso != tauIso.end(); 
	 ++iTauIso) {

      for (vector<string>::const_iterator iMuIso = muIso.begin(); iMuIso != muIso.end(); 
	   ++iMuIso) {

	string prefix("");
	if (*iMuIso == "nonIsoWData") prefix = "nonIsoW_";
	TFile file(("/data1/yohay/" + *iMuIso + "/analysis/" + prefix + "muHad" + *iTauIso + 
		    "Analysis" + *iMTBin + "_SingleMu_" + version + ".root").c_str());

	for (vector<string>::const_iterator iMethod = methods.begin(); iMethod != methods.end(); 
	     ++iMethod) {

	  TCanvas* canvas = NULL;
	  file.GetObject(("muHadMass3Mu" + *iMethod + "Canvas").c_str(), canvas);
	  if (canvas != NULL) {
	    TH1F* hist = (TH1F*)canvas->GetPrimitive(("muHadMass3Mu" + *iMethod).c_str());
	    if (hist != NULL) {
	      cout << *iMTBin << " " << *iTauIso << " " << *iMuIso << " " << *iMethod << " ";
	      cout << hist->Integral(0, 4);
	      if ((*iTauIso == "Iso") && (*iMuIso == "data")) {
		cout << endl;
	      }
	      else cout << " (" << hist->Integral(5, -1) << ")\n";
	    }
	  }
	}

	file.Close();
      }
    }
  }
}

void Shape(string var = "topness",  string region = "baseline", string channel = "allChannels", string filename = "plots/plotsProducer/1DStack.root"){
	TFile* fplots = new TFile(filename.c_str(),"READ");
	//retrieve canvas
	string cname = channel+"/"+region+"/"+var;
	TCanvas* c = fplots->Get(cname.c_str());
	//c->Draw();
	
	//signal
	string splot_name = "v:"+var+"|p:ttbar_2l|r:"+region+"|c:"+channel+"|t:1DEntries";
	TH1F* histo= (TH1F*) c->GetPrimitive(splot_name.c_str());
	cout<<"mean = "<<histo->GetMean()<<" rms = "<<histo->GetRMS()<<" skewness = "<<histo->GetSkewness()<<" kurtosis = "<<histo->GetKurtosis()<<endl;
	
}

/*compare the mu+had mass bin contents one to one for two versions of the analysis and print 
  discrepant bins*/
void compareVersions(const vector<vector<string> >& versions)
{
  const string histogramName("muHadMass");
  const string canvasName(histogramName + "Canvas");
  vector<vector<vector<Float_t> > > muHadMass;
  for (vector<vector<string> >::const_iterator iVersion = versions.begin(); 
       iVersion != versions.end(); ++iVersion) {
    vector<vector<Float_t> > muHadMassThisVersion;
    for (vector<string>::const_iterator iFile = iVersion->begin(); iFile != iVersion->end(); 
	 ++iFile) {
      vector<Float_t> muHadMassThisFile;
      TFile file(iFile->c_str());
      if (file.IsOpen()) {
	TCanvas* canvas = NULL;
	file.GetObject(canvasName.c_str(), canvas);
	if (canvas != NULL) {
	  TH1F* histogram = NULL;
	  histogram = (TH1F*)canvas->GetPrimitive(histogramName.c_str());
	  if (histogram != NULL) {
	    for (Int_t iBin = 0; iBin <= (histogram->GetNbinsX() + 1); ++iBin) {
	      muHadMassThisFile.push_back(histogram->GetBinContent(iBin));
	    }
	  }
	  else cerr << "Null histogram pointer\n";
	}
	else cerr << "Null canvas pointer\n";
	file.Close();
      }
      else cerr << "File not open\n";
      muHadMassThisVersion.push_back(muHadMassThisFile);
    }
    muHadMass.push_back(muHadMassThisVersion);
  }
  cout << "old new % difference\n";
  for (unsigned int iFile = 0; iFile < versions[0].size(); ++iFile) {
    cout << versions[0][iFile] << endl;
    for (unsigned int iBin = 0; iBin < muHadMass[0][0].size(); ++iBin) {
      if (muHadMass[0][iFile][iBin] != muHadMass[1][iFile][iBin]) {
	cout << "Bin " << iBin << ": " << muHadMass[0][iFile][iBin] << " ";
	cout << muHadMass[1][iFile][iBin] << " ";
	cout << ((muHadMass[0][iFile][iBin] - muHadMass[1][iFile][iBin])/
		 muHadMass[0][iFile][iBin])*100 << endl;
      }
    }
  }
}

TCanvas* getExtrapolFak(TString plotName, TString label, int verbose, TString outputFileFull, TString outputFileParton, TString outputFileHadron){

  // ============================ 
  //  Set Root Style
  // ============================

  TStyle myStyle("HHStyle","HHStyle");
  setHHStyle(myStyle);
  TGaxis::SetMaxDigits(2);
  myStyle.cd();
  gROOT->SetStyle("HHStyle");
  
  // open files
  TFile* fileFull   = TFile::Open(outputFileFull  , "READ");
  TFile* fileParton = TFile::Open(outputFileParton, "READ");
  TFile* fileHadron = TFile::Open(outputFileHadron, "READ");
  // dont associate new objects with file to be able to close it in the end
  gROOT->cd();
  // get canvas^3 for chosen cross section
  TCanvas* canvasFull   = (TCanvas*)(fileFull  ->Get("xSec/sysNo/"+plotName+"Norm")->Clone());
  TCanvas* canvasParton = (TCanvas*)(fileParton->Get("xSec/sysNo/"+plotName+"Norm")->Clone());
  TCanvas* canvasHadron = (TCanvas*)(fileHadron->Get("xSec/sysNo/"+plotName+"Norm")->Clone());
  
  // get data histos 
  TH1F* dataFull      = killEmptyBins((TH1F*)((canvasFull  ->GetPrimitive(plotName+"kData"))->Clone()));
  TH1F* dataRawParton = killEmptyBins((TH1F*)((canvasParton->GetPrimitive(plotName+"kData"))->Clone()));
  TH1F* dataRawHadron = killEmptyBins((TH1F*)((canvasHadron->GetPrimitive(plotName+"kData"))->Clone()));
  
  // use always the correct PS definition:
  //     hadron level for b-quarks and lepton
  //     parton level for all others
  TH1F* dataPS = ( (plotName.Contains("bq")||plotName.Contains("lep")) ? (TH1F*)dataRawHadron->Clone() : (TH1F*)dataRawParton->Clone() );

  // adjust style and labels
  TString PSlabel = ( (plotName.Contains("bq")||plotName.Contains("lep")) ? "hadron" : "parton" );
  int color = kBlue;
  if(PSlabel=="hadron") color-=4;
  dataPS->SetLineColor(color);
  dataPS->SetMarkerColor(color);
  histogramStyle(*dataFull, kData);
  dataFull->SetLineWidth(3);
  dataPS->SetLineWidth(3);
  if     (plotName=="lepPt") dataFull->GetXaxis()->SetRangeUser(0.,199.);
  else if(plotName=="bqPt" ) dataFull->GetXaxis()->SetRangeUser(0.,399.);
  else                       setXAxisRange(dataFull, plotName);
  double max=dataFull->GetMaximum();
  if(max<dataPS->GetMaximum()) max=dataPS->GetMaximum();
  dataFull->SetMaximum(1.3*max);
  dataFull->GetXaxis()->SetTitle(xSecLabelName(plotName));
  TString label2=label;
  TString label3="";
  if(label.Contains("/[GeV]")){
    label2.ReplaceAll("/[GeV]","");
    label3=" / [GeV]";
  }
  label2.ReplaceAll("/ ","");
  dataFull->GetYaxis()->SetTitle("#frac{1}{#sigma} #frac{d#sigma}{d"+label2+"}"+label3);
  dataFull->GetYaxis()->SetNoExponent(false);
  dataFull->GetXaxis()->SetNoExponent(true);
  dataFull->SetTitle("");

  // create legend
  TLegend *leg0 = new TLegend(0.65, 0.762, 0.95, 0.89);
  leg0->SetFillStyle(0);
  leg0->SetBorderSize(0);
  leg0->SetHeader("phase spaces");
  leg0->AddEntry(dataFull, "extrapol. parton lv","L");
  leg0->AddEntry(dataPS  , "restricted "+PSlabel+" lv","L");

  // create label
  TPaveText *headerlabel = new TPaveText();
  headerlabel -> SetX1NDC(gStyle->GetPadLeftMargin());
  headerlabel -> SetY1NDC(1.0-gStyle->GetPadTopMargin());
  headerlabel -> SetX2NDC(1.0-gStyle->GetPadRightMargin());
  headerlabel -> SetY2NDC(1.0);
  headerlabel -> SetTextFont(42);
  headerlabel -> AddText("comparing 2011 data results");
  headerlabel->SetFillStyle(0);
  headerlabel->SetBorderSize(0);
  headerlabel->SetTextSize(0.04);
  headerlabel->SetTextAlign(32);

  // create extrapolation factor / ratio canvas
  std::vector<TCanvas*> plotCanvas_;
  addCanvas(plotCanvas_);
  plotCanvas_[0]->cd();
  plotCanvas_[0]->Draw();
  dataFull->Draw("hist");
  dataPS->Draw("hist same");
  leg0->Draw("same");
  headerlabel->Draw("same");
  DrawDecayChLabel("e/#mu + Jets Combined");
  drawRatio(dataPS, dataFull, 0., 2.4, myStyle, verbose, std::vector<double>(0), PSlabel+" PS", "extrapolated", "hist", kBlack);

  // close files
  fileFull  ->Close();
  fileParton->Close();
  fileHadron->Close();
  
  // return
//.........这里部分代码省略.........

TCanvas* getRatio(TString plotName, int verbose, TString outputFile){
  // ============================
  //  Set Root Style
  // ============================
  TStyle myStyle("HHStyle","HHStyle");
  setHHStyle(myStyle);
  TGaxis::SetMaxDigits(2);
  myStyle.cd();
  gROOT->SetStyle("HHStyle");
  
  // draw data uncertainties as bands
  bool bands=true;

  // compare current with old result
  bool compare=false;
  TString oldResult="/afs/naf.desy.de/group/cms/scratch/tophh/tmp/OldCombination/";

  // open file
  TFile* file = TFile::Open(outputFile, "READ");
  // dont associate new objects with file to be able to close it in the end
  gROOT->cd();
  // get canvas for chosen cross section
  TCanvas* canvas = (TCanvas*)(file->Get("finalXSec/"+plotName+"Norm")->Clone());
  // GET DATA: with final errors from canvas
  TGraphAsymmErrors* dataRaw  = (TGraphAsymmErrors*)canvas->GetPrimitive("dataTotalError");
  TGraphAsymmErrors* dataStat = (TGraphAsymmErrors*)canvas->GetPrimitive("dataStatError");
  // GET DATA: create rebinned histo
  std::map< TString, std::vector<double> > binning_ = makeVariableBinning();
  int Nbins = std::abs(binning_[plotName][binning_[plotName].size()-1]-binning_[plotName][0])*100;
  if(plotName.Contains("topY")) Nbins/=10;
  if(verbose>1) std::cout << Nbins << std::endl;
  TH1F* datatemp= new TH1F("data"+plotName, "data"+plotName, Nbins, binning_[plotName][0], binning_[plotName][binning_[plotName].size()-1]);
  reBinTH1F(*datatemp, binning_[plotName], 0);
  // GET DATA: refill TGraphAsymmErrors to rebinned histo
  for(int bin=1; bin<=datatemp->GetNbinsX(); ++bin){
    if(verbose>1){
      std::cout << "bin: " << bin << std::endl;
      std::cout << dataRaw->GetY()[bin];
    }
    datatemp->SetBinContent(bin, dataRaw->GetY()[bin]);
    double err=dataRaw->GetErrorYhigh(bin);
    if(err<dataRaw->GetErrorYlow(bin)) err=dataRaw->GetErrorYlow(bin);
    if(verbose>1) std::cout << " +- " << err << std::endl;
    datatemp->SetBinError(bin, err);
  }
  // GET DATA: delete empty bins
  TH1F* data=killEmptyBins((TH1F*)datatemp->Clone(), verbose);
  data->GetXaxis()->SetTitle(xSecLabelName(plotName));
  if(verbose>1){
    for(int bin=1; bin<=data->GetNbinsX(); ++bin){
      std::cout << "bin: " << bin << std::endl;
      std::cout << data->GetBinContent(bin) << " +- " << data->GetBinError(bin) << std::endl;
    }
  }

  // GET THEORY: binned curves from canvas
  TH1F* plotNNLO     = (TH1F*)canvas->GetPrimitive(plotName+"nnlo"   );
  TH1F* plotMadGraph = (TH1F*)canvas->GetPrimitive(plotName          );
  TH1F* plotmcatnlo  = (TH1F*)canvas->GetPrimitive(plotName+"[email protected]");
  if(!plotmcatnlo) plotmcatnlo  = (TH1F*)canvas->GetPrimitive(plotName+"[email protected]");
  TGraphAsymmErrors* plotmcatnloerror = (TGraphAsymmErrors*)canvas->GetPrimitive(plotName+"[email protected]");
  TH1F* plotpowheg   = (TH1F*)canvas->GetPrimitive(plotName+"POWHEG");
  TH1F* plotpowhegherwig = (TH1F*)canvas->GetPrimitive(plotName+"POWHEGHERWIG");
  std::vector<TH1F*>hist_;
  // GET THEORY: delete empty bins
  // a) peturbative QCD
  TH1F* finalNNLO=0;
  if(plotNNLO        ){
    // delete empty bins
    TH1F* tempNNLO=killEmptyBins(plotNNLO, verbose);
    std::cout << tempNNLO->GetName() << std::endl;  
    // delete bins put of range
    int Nnnlobins = std::abs(binning_[plotName][binning_[plotName].size()-1]-binning_[plotName][0])*10;
    finalNNLO=new TH1F(tempNNLO->GetName(),tempNNLO->GetTitle(), Nnnlobins, binning_[plotName][0], binning_[plotName][binning_[plotName].size()-1]);
    reBinTH1F(*finalNNLO, binning_[plotName], 0);
    for(int bin=0; bin<=tempNNLO->GetNbinsX()+1; ++bin){
      double binlowedge=tempNNLO->GetBinLowEdge(bin); 
      if(plotName.Contains("topPt")&&binlowedge==1.) binlowedge=0.;
      //std::cout << "binlowedge: " << binlowedge << std::endl;      
      for(int binf=0; binf<=finalNNLO->GetNbinsX()+1; ++binf){
	//std::cout << "scanlowedge: " << finalNNLO->GetBinLowEdge(binf) << std::endl;
	if(binlowedge==finalNNLO->GetBinLowEdge(binf)){
	  //std::cout << "fits!" << std::endl;
	  finalNNLO->SetBinContent(binf, tempNNLO->GetBinContent(bin)) ;
	  finalNNLO->SetBinError(binf, tempNNLO->GetBinError(bin)); 
	}
      }
    }
  }
  // b) [email protected] errorbands
  if(plotmcatnloerror&&plotmcatnlo){
    TH1F*   plotmcatnloerror1 =(TH1F*)((killEmptyBins(plotmcatnlo, verbose))->Clone((TString)plotmcatnloerror->GetName()+"Up"));
    TH1F*   plotmcatnloerror2 =(TH1F*)((killEmptyBins(plotmcatnlo, verbose))->Clone((TString)plotmcatnloerror->GetName()+"Dn"));
    for(int p=0; p<plotmcatnloerror->GetN(); ++p){
      plotmcatnloerror1->SetBinContent(p, plotmcatnloerror->GetErrorYhigh(p)+plotmcatnloerror->GetY()[p]);
      plotmcatnloerror2->SetBinContent(p, plotmcatnloerror->GetY()[p]-plotmcatnloerror->GetErrorYlow(p));
    }
    plotmcatnloerror1->SetLineStyle(1);
    plotmcatnloerror2->SetLineStyle(1);
    hist_.push_back( killEmptyBins(plotmcatnloerror1, verbose) );
//.........这里部分代码省略.........

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);
//.........这里部分代码省略.........

void rootMacro_CleanJets_NTauDM()
{
  gStyle->SetOptStat(kFALSE);

  TFile infileCJ("/afs/cern.ch/user/k/ktos/GroupDir/CMSSW_7_6_3/src/AnalyzerGeneratorRecoVariousFunctions/Analyzer/BSUB/ggH125a9_GenTauDecayID_IndivCJ_NewDMFind_TauDMPlots_Scale_APR6/ggH125a9_GenTauDecayID_IndivCJ_NewDMFind_TauDMPlots_Scale_APR6_Plots.root");
  TFile infileRECO("/afs/cern.ch/user/k/ktos/GroupDir/CMSSW_7_6_3/src/AnalyzerGeneratorRecoVariousFunctions/Analyzer/BSUB/ggH125a9_GenTauDecayID_IndivRECO_NewDMFind_TauDMPlots_Scale_APR6/ggH125a9_GenTauDecayID_IndivRECO_NewDMFind_TauDMPlots_Scale_APR6_Plots.root");

  TFile *outFile = new TFile("combHist_CleanJets_h125a9_NTauDecayMode.root", "RECREATE");

cout << "Files Created" << endl;

  TCanvas *NTauDecayModeCJRecoCanvas = (TCanvas*)infileCJ.Get("NTauDecayModeRECO");
  TCanvas *NTauDecayModeRECORecoCanvas = (TCanvas*)infileRECO.Get("NTauDecayModeRECO");


  TCanvas *NTauDecayModeCJGenCanvas = (TCanvas*)infileCJ.Get("NTauDecayModeGEN");
  TCanvas *NTauDecayModeRECOGenCanvas = (TCanvas*)infileRECO.Get("NTauDecayModeGEN");

cout << "Got Canvases" << endl;

  TH1F* NTauDecayModeCJReco_ = (TH1F*)NTauDecayModeCJRecoCanvas->GetPrimitive("NTauDecayModeRECO");
  TH1F* NTauDecayModeRECOReco_ = (TH1F*)NTauDecayModeRECORecoCanvas->GetPrimitive("NTauDecayModeRECO");

  TH1F* NTauDecayModeCJGen_ = (TH1F*)NTauDecayModeCJGenCanvas->GetPrimitive("NTauDecayModeGEN");
  TH1F* NTauDecayModeRECOGen_ = (TH1F*)NTauDecayModeRECOGenCanvas->GetPrimitive("NTauDecayModeGEN");

cout << "Histograms assigned." << endl; 

  TCanvas NTauDecayModeRECO("NTauDecayModeRECO","",600,600);
  TCanvas NTauDecayModeGEN("NTauDecayModeGEN","",600,600);

cout << "Canvases created" << endl;

  NTauDecayModeCJReco_->SetLineColor(kBlack);
  NTauDecayModeRECOReco_->SetLineColor(kRed);
  
  NTauDecayModeCJGen_->SetLineColor(kBlack);
  NTauDecayModeRECOGen_->SetLineColor(kRed);
  
  NTauDecayModeCJReco_->SetMarkerColor(kBlack);
  NTauDecayModeRECOReco_->SetMarkerColor(kRed);
  
  NTauDecayModeCJGen_->SetMarkerColor(kBlack);
  NTauDecayModeRECOGen_->SetMarkerColor(kRed);
 
  NTauDecayModeCJReco_->GetXaxis()->SetTitle("Tau Decay Mode RECO");
  NTauDecayModeRECOReco_->GetXaxis()->SetTitle("Tau Decay Mode RECO");
 
  NTauDecayModeCJGen_->GetXaxis()->SetTitle("Tau Decay Mode Gen");
  NTauDecayModeRECOGen_->GetXaxis()->SetTitle("Tau Decay Mode GEN");
 
cout << "Attributes set." << endl;  

  leg = new TLegend(0.1,0.7,0.25,0.9);
  leg->AddEntry(NTauDecayModeCJReco_, "No CJ","L");
  leg->AddEntry(NTauDecayModeRECOReco_, "CJ","L");

  NTauDecayModeRECO.cd();
  NTauDecayModeCJReco_->Draw();
  NTauDecayModeRECOReco_->Draw("SAME");
  leg->Draw();

  NTauDecayModeGEN.cd();
  NTauDecayModeCJGen_->Draw();
  NTauDecayModeRECOGen_->Draw("SAME");
  leg->Draw();

cout << "Histograms Drawn" << endl;

  outFile->cd();
  NTauDecayModeRECO.Write();
  NTauDecayModeGEN.Write();
  outFile->Write();
  outFile->Close();

cout << "end" << endl;

}//rootMacro_BBA_combine

void makeResultTables(std::string decayChannel = "combined", bool extrapolate=true, bool hadron=false, bool addCrossCheckVariables=false, bool useBCC=false, int verbose=1){
  
  // ============================
  //  Set Root Style
  // ============================

  TStyle myStyle("HHStyle","HHStyle");
  setHHStyle(myStyle);
  TGaxis::SetMaxDigits(2);
  myStyle.cd();
  gROOT->SetStyle("HHStyle");

  // ============================
  //  Open file
  // ============================
  if(extrapolate==true) hadron=false; 
  TString filename="diffXSecTopSemi";
  if(decayChannel=="combined") filename+="Lep";
  else if(decayChannel=="electron") filename+="Elec";
  else if(decayChannel=="muon"    ) filename+="Mu";
  if(extrapolate) filename+="Parton";
  else{
    if(hadron) filename+="Hadron";
    else filename+="Parton";
    filename+="PhaseSpace";
  }
  filename+=".root";
  if(verbose>0) std::cout << "opening file " << filename << std::endl;
  TFile* file = TFile::Open(filename, "READ");
  if(!file){
    std::cout << "ERROR: can not open file " << filename << std::endl;
    exit(0);
  }
  // dont associate new objects with file to be able to close it in the end
  gROOT->cd();
  // ============================
  //  Get plots
  // ============================
  // variables to be processed
  std::vector<TString> xSecVariables_;
  // a) top and ttbar quantities
  if(!hadron){
    xSecVariables_.insert(xSecVariables_.end(), xSecVariablesKinFit, xSecVariablesKinFit + sizeof(xSecVariablesKinFit)/sizeof(TString));
  }
  // b) lepton and b-jet quantities
  if(hadron||!extrapolate){
    xSecVariables_.insert(xSecVariables_.end(), xSecVariablesFinalState    , xSecVariablesFinalState     + sizeof(xSecVariablesFinalState    )/sizeof(TString));
  }
  // c) cross check variables presently only available for parton level cross-sections
  if (addCrossCheckVariables && !hadron){
    xSecVariables_.insert( xSecVariables_.end(),   xSecVariablesCCVar,     xSecVariablesCCVar     + sizeof(xSecVariablesCCVar    )/sizeof(TString)    );
    xSecVariables_.insert( xSecVariables_.end(),   xSecVariablesCCVarNorm, xSecVariablesCCVarNorm + sizeof(xSecVariablesCCVarNorm)/sizeof(TString));
  }

  for(unsigned int i=0; i<xSecVariables_.size(); ++i){
    TString plotName=xSecVariables_[i];
    if(verbose>0) std::cout << std::endl << "variable: " << plotName << std::endl;
    // get canvas for chosen cross section
    TCanvas* canvas = (TCanvas*)(file->Get("finalXSec/"+plotName+"Norm")->Clone());
    if(!canvas){
      std::cout << "ERROR: can not load canvas finalXSec/"+plotName+"Norm" << std::endl;
      exit(0);
    }
    // GET DATA: with final errors from canvas
    TGraphAsymmErrors* dataTot  = (TGraphAsymmErrors*)canvas->GetPrimitive("dataTotalError");
    TGraphAsymmErrors* dataStat = (TGraphAsymmErrors*)canvas->GetPrimitive("dataStatError" );
    TH1F* binned = (TH1F*)canvas->GetPrimitive(plotName);
    TH1F* binnedMCatNLO = (TH1F*)canvas->GetPrimitive(plotName+"[email protected]");
    TH1F* binnedPowheg  = (TH1F*)canvas->GetPrimitive(plotName+"POWHEG");
    TH1F* binnedPowhegHerwig  = (TH1F*)canvas->GetPrimitive(plotName+"POWHEGHERWIG");
    TH1F* binnedNNLO    = (TH1F*)canvas->GetPrimitive(plotName+"nnlo");
    if(!dataTot){
      std::cout << "ERROR: can not load TGraphAsymmErrors dataTotalError in canvas finalXSec/"+plotName+"Norm" << std::endl;
      exit(0);
    }
    if(!dataStat){
      std::cout << "ERROR: can not load TGraphAsymmErrors dataStatError in canvas finalXSec/"+plotName+"Norm" << std::endl;
      exit(0);
    }
    if(!binned){
      std::cout << "ERROR: can not load TH1F topPt in canvas finalXSec/"+plotName << std::endl;
      exit(0);
    }
    // define range of relevant plots
    // INFO: keep this consistent with the range as defined in setXAxisRange
    //       and makevariableBinning in basicFunctions.h
    double xMin=-999999999;
    double xMax= 999999999;
    if(plotName.Contains     ("topPt"    )){ xMin=0.   ; xMax=401. ;} 
    else if(plotName.Contains("topY"     )){ xMin=-2.51; xMax=2.51 ;}
    else if(plotName.Contains("ttbarY"   )){ xMin=-2.51; xMax=2.51 ;}
    else if(plotName.Contains("ttbarMass")){ xMin=344. ; xMax=1601.;}
    else if(plotName.Contains("ttbarPt"  )){ xMin=0.   ; xMax=301. ;}
    else if(plotName.Contains("lepPt"    )){ xMin=29   ; xMax=201. ;}
    else if(plotName.Contains("lepEta"   )){ xMin=-2.11; xMax=2.11 ;}
    else if(plotName.Contains("bqPt"     )){ xMin=29.  ; xMax=401. ;}
    else if(plotName.Contains("bqEta"    )){ xMin=-2.41; xMax=2.41 ;}
    else if(plotName.Contains("bbbarPt"  )){ xMin=0.   ; xMax=800. ;}
    else if(plotName.Contains("bbbarMass")){ xMin=0.   ; xMax=1200.;}
    else if(plotName.Contains("ttbarDelPhi" )){ xMin=0. ; xMax=3.16;}
//.........这里部分代码省略.........

void estimateSignalFitPerformance()
{
  //open the ROOT efficiency file
  TFile ROOTFile(efficiencyFile.c_str());
  if (!ROOTFile.IsOpen()) {
    cerr << "Error opening file " << efficiencyFile << ".\n";
    return;
  }

  //get numBackgroundFail, numBackgroundPass, numSignalAll, and efficiency
  RooFitResult* fitResult = NULL;
  ROOTFile.GetObject("PhotonToIDEB/unbinned/probe_eta_bin0__probe_nJets05_bin0__gaussPlusLinear/fitresults", fitResult);
  Double_t efficiency = 0.0;
  Double_t efficiencyError = 0.0;
  Double_t numBackgroundFail = 0.0;
  Double_t numBackgroundFailError = 0.0;
  Double_t numBackgroundPass = 0.0;
  Double_t numBackgroundPassError = 0.0;
  Double_t numSignalAll = 0.0;
  Double_t numSignalAllError = 0.0;
  if (fitResult != NULL) {
    RooRealVar* theEfficiency = (RooRealVar*)fitResult->floatParsFinal().find("efficiency");
    efficiency = theEfficiency->getVal();
    efficiencyError = theEfficiency->getError();
    RooRealVar* theNumBackgroundFail = (RooRealVar*)fitResult->floatParsFinal().find("numBackgroundFail");
    numBackgroundFail = theNumBackgroundFail->getVal();
    numBackgroundFailError = theNumBackgroundFail->getError();
    RooRealVar* theNumBackgroundPass = (RooRealVar*)fitResult->floatParsFinal().find("numBackgroundPass");
    numBackgroundPass = theNumBackgroundPass->getVal();
    numBackgroundPassError = theNumBackgroundPass->getError();
    RooRealVar* theNumSignalAll = (RooRealVar*)fitResult->floatParsFinal().find("numSignalAll");
    numSignalAll = theNumSignalAll->getVal();
    numSignalAllError = theNumSignalAll->getError();
  }
  else {
    cerr << "Error getting RooFitResult PhotonToIDEB/unbinned/probe_eta_bin0__probe_nJets05_bin0__gaussPlusLinear/fitresults from file ";
    cerr << efficiencyFile << ".\n";
  }

  //get integrals of tag-pass and tag-fail distributions
  TCanvas* fitCanvas = NULL;
  ROOTFile.GetObject("PhotonToIDEB/unbinned/probe_eta_bin0__probe_nJets05_bin0__gaussPlusLinear/fit_canvas", fitCanvas);
  Double_t tagPassIntegral = 0;
  Double_t tagFailIntegral = 0;
  if (fitCanvas != NULL) {
    fitCanvas->cd(1);
    RooHist* tagPassDistribution = NULL;
    tagPassDistribution = (RooHist*)((TCanvas*)fitCanvas->GetPrimitive("fit_canvas_1"))->GetPrimitive("h_data");
    fitCanvas->cd(2);
    RooHist* tagFailDistribution = NULL;
    tagFailDistribution = (RooHist*)((TCanvas*)fitCanvas->GetPrimitive("fit_canvas_2"))->GetPrimitive("h_data");
    RooHist* blah = NULL;
    blah = (RooHist*)((TCanvas*)fitCanvas->GetPrimitive("fit_canvas_3"))->GetPrimitive("h_data");
    cerr << blah->Integral() << endl;
    if ((tagPassDistribution != NULL) && (tagFailDistribution != NULL)) {
      tagPassIntegral = tagPassDistribution->Integral()*/*1.796*/1.844;
      tagFailIntegral = tagFailDistribution->Integral()*/*1.796*/1.844;
    }
    else cerr << "Error: could not get RooPlots.\n";
  }
  else {
    cerr << "Error getting TCanvas PhotonToIDEB/unbinned/probe_eta_bin0__probe_nJets05_bin0__gaussPlusLinear/fit_canvas from file ";
    cerr << efficiencyFile << ".\n";
  }

  //close file
  ROOTFile.Close();

  //subtract fitted background from integral
  Double_t tagPassNumBkgSubtractedEvts = tagPassIntegral - numBackgroundPass;
  Double_t tagPassNumBkgSubtractedEvtsError = numBackgroundPassError;
  Double_t tagFailNumBkgSubtractedEvts = tagFailIntegral - numBackgroundFail;
  Double_t tagFailNumBkgSubtractedEvtsError = numBackgroundFailError;

  //calculate fitted signal
  Double_t tagPassNumFittedSignal = efficiency*numSignalAll;
  Double_t tagPassNumFittedSignalError = tagPassNumFittedSignal*sqrt(((efficiencyError*efficiencyError)/(efficiency*efficiency)) + 
								     ((numSignalAllError*numSignalAllError)/(numSignalAll*numSignalAll)));
  Double_t tagFailNumFittedSignal = (1.0 - efficiency)*numSignalAll;
  Double_t tagFailNumFittedSignalError = tagFailNumFittedSignal*
    sqrt(((efficiencyError*efficiencyError)/((1.0 - efficiency)*(1.0 - efficiency))) + 
	 ((numSignalAllError*numSignalAllError)/(numSignalAll*numSignalAll)));

  //calculate difference between signal fit result and background subtracted integral
  Double_t tagPassDifference = tagPassNumBkgSubtractedEvts - tagPassNumFittedSignal;
  Double_t tagPassDifferenceError = sqrt(tagPassNumBkgSubtractedEvtsError*tagPassNumBkgSubtractedEvtsError + 
					 tagPassNumFittedSignalError*tagPassNumFittedSignalError);
  Double_t tagFailDifference = tagFailNumBkgSubtractedEvts - tagFailNumFittedSignal;
  Double_t tagFailDifferenceError = sqrt(tagFailNumBkgSubtractedEvtsError*tagFailNumBkgSubtractedEvtsError + 
					 tagFailNumFittedSignalError*tagFailNumFittedSignalError);

  //compare signal fit result to background subtracted integral
  cout << "Tag pass signal fit: " << tagPassNumFittedSignal << " +/- " << tagPassNumFittedSignalError << endl;
  cout << "Tag pass background subtracted integral: " <<  tagPassNumBkgSubtractedEvts << " +/- " << tagPassNumBkgSubtractedEvtsError;
  cout << endl;
  cout << "Difference: " << tagPassDifference << " +/- " << tagPassDifferenceError << endl;
  cout << "Tag fail signal fit: " << tagFailNumFittedSignal << " +/- " << tagFailNumFittedSignalError << endl;
  cout << "Tag fail background subtracted integral: " <<  tagFailNumBkgSubtractedEvts << " +/- " << tagFailNumBkgSubtractedEvtsError;
  cout << endl;
  cout << "Difference: " << tagFailDifference << " +/- " << tagFailDifferenceError << endl;
//.........这里部分代码省略.........