C++ THStack::Write方法代码示例

void createInputs(int n = 2) 
{
   for(UInt_t i = 0; i < (UInt_t)n; ++i ) {
      TFile *file = TFile::Open(TString::Format("input%d.root",i),"RECREATE");
      TH1F * h = new TH1F("h1","",10,0,100);
      h->Fill(10.5); h->Fill(20.5);
 
      Int_t nbins[5];
      Double_t xmin[5];
      Double_t xmax[5];
      for(UInt_t j = 0; j < 5; ++j) {
         nbins[j] = 10; xmin[j] = 0; xmax[j] = 10;
      }
      THnSparseF *sparse = new THnSparseF("sparse", "sparse", 5, nbins, xmin, xmax);
      Double_t coord[5] = {0.5, 1.5, 2.5, 3.5, 4.5};
      sparse->Fill(coord);
      sparse->Write();
      
      THStack *stack = new THStack("stack","");
      h = new TH1F("hs_1","",10,0,100);
      h->Fill(10.5); h->Fill(20.5);
      h->SetDirectory(0);
      stack->Add(h);
      h = new TH1F("hs_2","",10,0,100);
      h->Fill(30.5); h->Fill(40.5);
      h->SetDirectory(0);
      stack->Add(h);
      stack->Write();

      TGraph *gr = new TGraph(3);
      gr->SetName("exgraph");
      gr->SetPoint(0,1,1);
      gr->SetPoint(1,2,2);
      gr->SetPoint(2,3,3);
      
      gr->Write();
      
      TTree *tree = new TTree("tree","simplistic tree");
      Int_t data = 0;
      tree->Branch("data",&data);
      for(Int_t l = 0; l < 2; ++l) {
         data = l;
         tree->Fill();
      }
      
      file->Write();
      delete file;
   }
}

//.........这里部分代码省略.........
    h_data->SetMarkerSize(1);
    h_data->SetMarkerStyle(24);
    h_data->SetTitle("");
    h_data->Draw("e");
    if(xmin>-9999 && xmax>-9999) {
        h_data->GetXaxis()->SetRangeUser(xmin,xmax);
    }
    cout << "Data integral = " << h_data->Integral() << endl;
    hs->Draw("histsame");
    h_data->Draw("esame");
    leg->Draw("same");


    std::string dirName = "compareDataMC_" + leptonName;
    gSystem->mkdir(dirName.data());

    std::string filename;
    std::string psname ;
    psname = dirName+ "/overlay_" + histoName;
    filename = psname + ".eps";
    c1->Print(filename.data());
    filename = psname + ".gif";
    c1->Print(filename.data());
    filename = psname + ".pdf";
    c1->Print(filename.data());


    // study the ratios

    TCanvas* c2 = new TCanvas("c2","",700,0,700,1000);
    c2->Divide(1,2,0.01,0);
    c2->cd(1);
    gPad->SetTopMargin(0.01);
    gPad->SetBottomMargin(0);
    gPad->SetRightMargin(0.04);

    h_data->Draw("e");
    if(xmin>-9999 && xmax>-9999) {
        h_data->GetXaxis()->SetRangeUser(xmin,xmax);
    }
    cout << h_data->GetName() << " integral = " << h_data->Integral() << endl;
    hs->Draw("histsame");
    h_data->Draw("esame");
    leg->Draw("same");


    c2->cd(2);
    gStyle->SetStatW       (0.3);
    gStyle->SetStatH       (0.3);
    gStyle->SetStatX       (0.879447);
    gStyle->SetStatY       (0.939033);
    gStyle->SetStatFontSize(0.05);
    gStyle->SetStatBorderSize(0);
    gPad->SetRightMargin(0.04);
    gPad->SetTopMargin(0);
    gPad->SetBottomMargin(0.2);
    gPad->SetTickx();
    gStyle->SetOptFit(1);

    TH1D* hratio = (TH1D*)h_data->Clone("hratio");
    hratio->Reset();
    hratio->Divide(h_data,h_all,1.0,1.0);
    hratio->SetTitle("");
    hratio->SetMaximum(1.5);
    hratio->SetMinimum(0.5);
    hratio->SetTitleOffset(1.2,"Y");
    hratio->GetXaxis()->SetTitle(h_data->GetXaxis()->GetTitle());
    hratio->Draw("e1");

    cout << "( " << h_data->GetBinContent(maxBin) << "+-" << h_data->GetBinError(maxBin) << " )/("
         << h_all->GetBinContent(maxBin) << "+-" << h_all->GetBinError(maxBin) <<  ")= "
         << hratio->GetBinContent(maxBin) << "+-" << hratio->GetBinError(maxBin) << endl;

    psname = dirName+ "/ratio_" + histoName;
    filename = psname + ".eps";
    c2->Print(filename.data());
    filename = psname + ".gif";
    c2->Print(filename.data());
    filename = psname + ".pdf";
    c2->Print(filename.data());


    std::string remword  =".txt";
    size_t pos  = inputFile.find(remword);

    if(pos!= std::string::npos)
        inputFile.replace(pos,remword.length(),"");

    std::string command = "recreate";

    if(update)command ="update";

    TFile* outFile = new TFile(Form("combined_%s.root",inputFile.data()),command.data());
    h_all->Write();
    hs->Write();
    h_data->Write();
    outFile->Close();


}

//.........这里部分代码省略.........
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("Unrolled%s.pdf",Label.c_str()));

  // Unrolled plots in percentages
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.85,0.20,0.95,0.80);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);

  for (Int_t i = histMRRsq.size()-1 ; i >= 0; --i) {
    if (hasSignal && i==0) {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "L");
    } else {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "F");
    }
  }
  stackUnrolledPercentage->Draw();
  stackUnrolledPercentage->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage->GetHists()->At(0)))->GetXaxis()->GetTitle());
  // stackUnrolledPercentage->GetHistogram()->GetXaxis()->SetRangeUser(0, 35);
  stackUnrolledPercentage->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage->GetHists()->At(0)))->GetYaxis()->GetTitle());
  if(hasSignal) histUnrolledPercentage[0]->Draw("same hist"); 
  legend->Draw();
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("UnrolledPercentage%s.pdf",Label.c_str()));

  // Unrolled plots in sideband vs signal box
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.85,0.20,0.95,0.80);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);

  for (Int_t i = histMRRsq.size()-1 ; i >= 0; --i) {
    if (hasSignal && i==0) {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "L");
    } else {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "F");
    }
  }
  stackUnrolled2bins->Draw();
  stackUnrolled2bins->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackUnrolled2bins->GetHists()->At(0)))->GetXaxis()->GetTitle());
  stackUnrolled2bins->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackUnrolled2bins->GetHists()->At(0)))->GetYaxis()->GetTitle());
  if(hasSignal) histUnrolled2bins[0]->Draw("same hist");
  legend->Draw();
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("Unrolled2bins%s.pdf",Label.c_str()));

  // Unrolled plots in sideband vs signal box in fractions
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.7,0.23,0.90,0.88);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);

  for (Int_t i = histMRRsq.size()-1 ; i >= 0; --i) {
    if (hasSignal && i==0) {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "L");
    } else {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "F");
    }
  }

  stackUnrolledPercentage2bins->Draw();
  stackUnrolledPercentage2bins->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage2bins->GetHists()->At(0)))->GetXaxis()->GetTitle());
  stackUnrolledPercentage2bins->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage2bins->GetHists()->At(0)))->GetYaxis()->GetTitle());
  stackUnrolledPercentage2bins->GetHistogram()->GetXaxis()->SetBinLabel(1, "Sideband");
  stackUnrolledPercentage2bins->GetHistogram()->GetXaxis()->SetBinLabel(2, "Signal Sensitive Region");
  if(hasSignal) histUnrolledPercentage2bins[0]->Draw("same hist");
  legend->Draw();
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("UnrolledPercentage2bins%s.pdf",Label.c_str()));


   //--------------------------------------------------------------------------------------------------------------
  // Output
  //==============================================================================================================
  TFile *file = TFile::Open(("RazorPlots"+Label+".root").c_str(), "RECREATE");
  file->cd();

  for(int i=0; i<int(inputfiles.size()); i++) {
    file->WriteTObject(histMR[i], Form("histMR_%s",processLabels[i].c_str()), "WriteDelete");
    file->WriteTObject(histRsq[i], Form("histRsq_%s",processLabels[i].c_str()), "WriteDelete");
    file->WriteTObject(histMRRsq[i], Form("histMRRsq_%s",processLabels[i].c_str()), "WriteDelete");
    histUnrolled[i]->Write();  
    histUnrolled2bins[i]->Write();  
    histUnrolledPercentage[i]->Write();  
    histUnrolledPercentage2bins[i]->Write();  
  }
  
  stackUnrolled->Write();
  stackUnrolled2bins->Write();
  stackUnrolledPercentage->Write();
  stackUnrolledPercentage2bins->Write();
 }

void Stacker_onePoint(TString path,TString cut,double ME_Br,double EM_Br)
{
	TString MCSamples[11]={"wt", "ZZ","H2WWleptonic", "WlepZmue","WincZtautau",
			"H2tt", "ttbar", "WWleptonic",  "Z2tt","H2tm","H2te"};
	Int_t   MCcolors[11]={ kMagenta+3, kMagenta,kBlue+1,kOrange+1,kOrange,
			kGreen, kRed, kYellow, kCyan, kBlack, kBlack};


	TCanvas* c1 = new TCanvas("canvasLog"+cut,"canvasLog"+cut,600,600);
	c1->SetLogy();
//	TCanvas* c2 = new TCanvas("c_diff","c_diff",600,600);
	TCanvas* c3 = new TCanvas("c_ratio"+cut,"c_ratio"+cut,600,600);
	TCanvas* c5 = new TCanvas("canvasStack2"+cut,"canvasStack2"+cut,600,600); c5->SetLogy();
	TCanvas* c4 = new TCanvas("canvasStack"+cut,"canvasStack"+cut,600,600);	c4->SetLogy();
	TCanvas* c6 = new TCanvas("canvasEM_ME","canvasEM_ME",600,600); c6->SetLogy();

	TLegend* leg = new TLegend(0.5,0.7,0.7,0.9);
	leg->SetFillColor(kWhite);
	leg->SetBorderSize(1);

	THStack* hs = new THStack("ME_EM","ME_EM");
	THStack* hsStackedME = new THStack("stacked_ME",";"
			"M_{Collinear} (GeV);Events / 4 GeV");
	THStack* hsStackedEM = new THStack("stacked_EM","e#mu Channel;"
			"M_{Collinear} (GeV);Events / 4 GeV");
	THStack* hsEM_ME = new THStack("stacked_EM_ME",";"
				"M_{Collinear} (GeV);Events / 4 GeV");

	TH1D* EM_sum = new TH1D("EM_sum","bla",250,0,500);
	TH1D* ME_sum = new TH1D("ME_sum","bla",250,0,500);

	TFile* fs = new TFile(path+"H2tm"+cut+".root");
	TH1D* signal_ME = (TH1D*)fs->Get("ME_Mcoll"); signal_ME->SetLineColor(kBlack);
	signal_ME->SetLineStyle(2);
	signal_ME->SetName("signal_ME_Mcoll");
	TH1D* signal_EM = (TH1D*)fs->Get("EM_Mcoll"); signal_EM->SetLineColor(kBlack);
	signal_EM->SetLineStyle(2);
	signal_EM->SetName("signal_EM_Mcoll");

	double signalME_c = ME_Br*10;
	double signalEM_c = EM_Br*10;

	double c[11]={1, 1,1,1,1,1,1,1,1,signalME_c,signalEM_c};

	//group BG
	//diboson
	TFile* f_ZZ = new TFile(path+"ZZ"+cut+".root");
	TFile* f_WZ1 = new TFile(path+"WlepZmue"+cut+".root");
	TFile* f_WZ2 = new TFile(path+"WincZtautau"+cut+".root");
	TFile* f_WW = new TFile(path+"WWleptonic"+cut+".root");

	TH1D* h_ME_ZZ = (TH1D*)f_ZZ->Get("ME_Mcoll");h_ME_ZZ->Scale(c[1]);
	TH1D* h_EM_ZZ = (TH1D*)f_ZZ->Get("EM_Mcoll");h_EM_ZZ->Scale(c[1]);
	TH1D* h_ME_WZ1 = (TH1D*)f_WZ1->Get("ME_Mcoll");h_ME_WZ1->Scale(c[3]);
	TH1D* h_EM_WZ1 = (TH1D*)f_WZ1->Get("EM_Mcoll");h_EM_WZ1->Scale(c[3]);
	TH1D* h_ME_WZ2 = (TH1D*)f_WZ2->Get("ME_Mcoll");h_ME_WZ2->Scale(c[4]);
	TH1D* h_EM_WZ2 = (TH1D*)f_WZ2->Get("EM_Mcoll");h_EM_WZ2->Scale(c[4]);
	TH1D* h_ME_WW = (TH1D*)f_WW->Get("ME_Mcoll");h_ME_WW->Scale(c[7]);
	TH1D* h_EM_WW = (TH1D*)f_WW->Get("EM_Mcoll");h_EM_WW->Scale(c[7]);

	TH1D* Diboson_ME = new TH1D("Diboson_ME","Diboson_ME",250,0,500);
	Diboson_ME = (TH1D*)h_ME_ZZ->Clone("ME_Mcoll");
	Diboson_ME->Add(h_ME_WZ1);Diboson_ME->Add(h_ME_WZ2);Diboson_ME->Add(h_ME_WW);
	TH1D* Diboson_EM = new TH1D("Diboson_EM","Diboson_EM",250,0,500);
	Diboson_EM = (TH1D*)h_ME_ZZ->Clone("EM_Mcoll");
	Diboson_EM->Add(h_EM_WZ1);Diboson_EM->Add(h_EM_WZ2);Diboson_EM->Add(h_EM_WW);

	TFile *diboson_out = new TFile(path+"Diboson"+cut+".root","RECREATE");
	Diboson_ME->Write();
	Diboson_EM->Write();
	diboson_out->Close();

	//Higgs SM
	TFile* f_H2WW = new TFile(path+"H2WWleptonic"+cut+".root");
	TFile* f_H2tt = new TFile(path+"H2tt"+cut+".root");
	TH1D* h_ME_H2WW = (TH1D*)f_H2WW->Get("ME_Mcoll");h_ME_H2WW->Scale(c[2]);
	TH1D* h_EM_H2WW = (TH1D*)f_H2WW->Get("EM_Mcoll");h_EM_H2WW->Scale(c[2]);
	TH1D* h_ME_H2tt = (TH1D*)f_H2tt->Get("ME_Mcoll");h_ME_H2tt->Scale(c[5]);
	TH1D* h_EM_H2tt = (TH1D*)f_H2tt->Get("EM_Mcoll");h_EM_H2tt->Scale(c[5]);

	TH1D* SMHiggs_ME = new TH1D("SMHiggs_ME","SMHiggs_ME",250,0,500);
	SMHiggs_ME = (TH1D*)h_ME_H2WW->Clone("ME_Mcoll");
	SMHiggs_ME->Add(h_ME_H2tt);
	TH1D* SMHiggs_EM = new TH1D("SMHiggs_EM","SMHiggs_EM",250,0,500);
	SMHiggs_EM = (TH1D*)h_EM_H2WW->Clone("EM_Mcoll");
	SMHiggs_EM->Add(h_EM_H2tt);

	TFile *SMHiggs_out = new TFile(path+"SMHiggs"+cut+".root","RECREATE");
	SMHiggs_ME->Write();
	SMHiggs_EM->Write();
	SMHiggs_out->Close();

	TString MCSamplesGrouped[5]={"wt",	"SMHiggs", "ttbar", "Diboson",  "Z2tt"};
	Int_t   MCcolorsGrouped[5]={ kMagenta+3,kGreen, kRed, kYellow, kCyan};
	double cGrouped[5]={1, 1,1,1,1};


	//BG
	for(int i=0; i<5; i++)
		{
//.........这里部分代码省略.........

int main(){

  char tempname[200];
  vector<TFile *> T_inputfilevec,WJet_inputfilevec, TTbar_inputfilevec, ZJet_inputfilevec, QCD_inputfilevec;
  map<int, string> cutname, histname, Hname;
  TFile *file, *file2, *file30, *file3;
  TH1D *temphist, *temphist2, * temphist30, *temphistI, *temphistII, *temphistIII, *temphistI_lowDphi, *temphistII_lowDphi, *temphistIII_lowDphi;
  THStack * tempstack;
  TDirectory *cdtoitt, *cdtoit;
  Selection2 * sel = new Selection2();
  cutname = sel->cutName();
  map<int, string> Ttype, WJettype, TTbartype, ZJettype, QCDtype;
  WJettype[0]="allEvents";
  TTbartype[0]="allEvents";
  Ttype[0]="allEvents";
  ZJettype[0]="allEvents";
  QCDtype[0]="allEvents";

  int qcdHT=7, zjnHT=4, ttbarnHT=1, wjnHT=7, tnHT=4;



// ..............................................................................................................//
// QCD Section
// .............................................................................................................//

  for(int i=1; i<=qcdHT ; i++){
    if(i==1)sprintf(tempname,"../results_filelist_Spring15_QCD_HT_200_300_.root");
    else if(i==2)sprintf(tempname,"../results_filelist_Spring15_QCD_HT_300_500_.root");
    else if(i==3)sprintf(tempname,"../results_filelist_Spring15_QCD_HT_500_700_.root");
    else if(i==4)sprintf(tempname,"../results_filelist_Spring15_QCD_HT_700_1000_.root");
    else if(i==5)sprintf(tempname,"../results_filelist_Spring15_QCD_HT_1000_1500_.root");
    else if(i==6)sprintf(tempname,"../results_filelist_Spring15_QCD_HT_1500_2000_.root");
    else if(i==7)sprintf(tempname,"../results_filelist_Spring15_QCD_HT_2000_Inf_.root");
    else{cout << " Error!! There are only 4 QCD ht binned sample " << endl;}
    QCD_inputfilevec.push_back(TFile::Open(tempname,"R"));
  }//end of loop over HTbins


  // Stack
  tempstack = new THStack("stack","Binned Sample Stack");
  sprintf(tempname,"results_filelist_Spring15_QCD_stacked.root");
  file = new TFile(tempname,"RECREATE");
  histname.clear();
  histname[0]="weight";
  histname[1]="HT";
  histname[2]="MHT";
  histname[3]="NJet";
  histname[4]="NBtag";

  Hname.clear();
  Hname[0]="yield_tauId";
  Hname[1]="yield_tauId_trk";
  Hname[2]="cutflow_preselection";


  for(int j=0; j< Hname.size(); j++){

    for(int i=0; i<qcdHT ; i++){                                                  // loop over different HT bins

      sprintf(tempname,"%s",(Hname[j]).c_str());
      temphist = (TH1D *) QCD_inputfilevec.at(i)->Get(tempname)->Clone();
      //if (luminosity>0&&!doScale) temphist->Scale(scalefactor);
      temphist->SetFillColor(i+2);
      tempstack->Add(temphist);

    }//end of loop over HTbins 1..7
    sprintf(tempname,"%s",(Hname[j]).c_str());
    tempstack->Write(tempname);

    delete tempstack;
    tempstack = new THStack("stack","Binned Sample Stack");

  }

    for(map<int , string >::iterator itt=QCDtype.begin(); itt!=QCDtype.end();itt++){        // loop over different event types

    cdtoitt = file->mkdir((itt->second).c_str());
    cdtoitt->cd();

    for(map<int , string >::iterator it=cutname.begin(); it!=cutname.end();it++){   // loop over different cutnames

      cdtoit =  cdtoitt->mkdir((it->second).c_str());
      cdtoit->cd();

      for(int j=0; j<histname.size(); j++){                                        // loop over different histograms

        for(int i=0; i<qcdHT ; i++){                                                  // loop over different HT bins

          //cout << "================================" << endl;
          //cout << "HT#: " <<i << ", WJtype: " << itt->second << ", cutname: " << it->second << ", hist#: " << j << endl;
          sprintf(tempname,"%s/%s/%s_%s_%s",(itt->second).c_str(),(it->second).c_str(),(histname[j]).c_str(),(it->second).c_str(),(itt->second).c_str());
          temphist = (TH1D *) QCD_inputfilevec.at(i)->Get(tempname)->Clone();
          //if (luminosity>0&&!doScale) temphist->Scale(scalefactor);
          temphist->SetFillColor(i+2);
          tempstack->Add(temphist);

        }//end of loop over HTbins 1..7

        sprintf(tempname,"%s_%s_%s",histname[j].c_str(),(it->second).c_str(),(itt->second).c_str());
//.........这里部分代码省略.........

//.........这里部分代码省略.........
//sprintf(tempname,"PhaseII4_t2cc350340_14TEV_140PileUp_00.root");
//sprintf(tempname,"PhaseII4_t2cc350310_14TEV_140PileUp_00.root");
//sprintf(tempname,"PhaseII4_Stop_CharmLSPv4_14TEV_140PileUp_00.root");
sprintf(tempname,"PhaseII4_StauC_14TEV_140PileUp.root");
file = new TFile(tempname,"RECREATE");
 for(map<int , string >::iterator itt=sigtype.begin(); itt!=sigtype.end();itt++){        // loop over different event types
    cdtoitt = file->mkdir((itt->second).c_str());
    cdtoitt->cd();
int c=0;
    for(map<int , string >::iterator it=cutname.begin(); it!=cutname.end();it++){   // loop over different cutnames
      cdtoit =  cdtoitt->mkdir((it->second).c_str());
      cdtoit->cd();
      for(int j=0; j<histname.size(); j++){                                        // loop over different histograms
        for(int i=0; i<Sig_nHT ; i++){                                                  // loop over different HT bins
sprintf(tempname,"%s/%s/%s_%s_%s",(itt->second).c_str(),(it->second).c_str(),(histname[j]).c_str(),(it->second).c_str(),(itt->second).c_str());
temphist = (TH1D *) Sig_inputfilevec.at(i)->Get(tempname)->Clone();
temphist->Scale(Sig_scalevec[i]);
if(histname[j]=="MHT"){
Sig_numberofevents+=(double)temphist->GetSumOfWeights();
} 
temphist->SetFillColor(i+2);
tempstack->Add(temphist);


               }//end of loop over HTbins 1..7
if(histname[j]=="MHT"){
if(itt->second=="allEvents"){
yieldmap[c].push_back(Sig_numberofevents);
}

}
Sig_numberofevents=0;
        sprintf(tempname,"%s_%s_%s",histname[j].c_str(),(it->second).c_str(),(itt->second).c_str());
        tempstack->Write(tempname);
        delete tempstack;
        tempstack = new THStack("stack","Binned Sample Stack");
      }//end of loop over histograms
 c+=1;   }//end of loop over cutnames
  }//end of loop over event types
file->Close();

//BJ Section//BJ Section//BJ Section//BJ Section//BJ Section//BJ Section//BJ Section//BJ Section//BJ Section//BJ Section//BJ Section

 //build a vector of scale factors
  //first load the cross sections into a vector
  BJ_xs_vec.push_back(34409.92339);
  BJ_xs_vec.push_back(2642.85309);
  BJ_xs_vec.push_back(294.12311);
  BJ_xs_vec.push_back(25.95000);
  BJ_xs_vec.push_back(2.42111);
  BJ_xs_vec.push_back(0.22690);
  BJ_xs_vec.push_back(0.02767);

  double BJ_numberofevents =0;
  const int bjnHT = 7;   // Total number of HT bin samples


  for(int i=1; i<=bjnHT ; i++){
    sprintf(tempname,"../Results/results_PhaseII4_BJ_14TEV_HT%d_140PileUp.root",i);
    file = new TFile(tempname, "R");
    sprintf(tempname,"allEvents/RA2nocut/MHT_RA2nocut_allEvents");
    tempvalue = (luminosity*BJ_xs_vec[i-1])/((* (TH1D* ) file->Get(tempname)).GetEntries());
    BJ_scalevec.push_back(tempvalue);
  }//end of loop over HTbins 
  std::cout << "normalization scale factor determination done" << std::endl;
for(int i=1; i<=bjnHT; i++){

//.........这里部分代码省略.........
  //////////////////
  stackMR->Add(histMRAllBkg);
  stackMR->Add(histMRQCD);
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.50,0.54,0.90,0.84);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);
  legend->SetFillStyle(0);
  stackMR->Draw();
  stackMR->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackMR->GetHists()->At(0)))->GetXaxis()->GetTitle());
  stackMR->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackMR->GetHists()->At(0)))->GetYaxis()->GetTitle());
  stackMR->Draw("");
  histMRData->Draw("same PE");
  legend->Draw();
  cv->SetLogy();
  cv->SaveAs(Form("MRStack_QCD_%s.pdf",Label.c_str()));

  ///////////////////////
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.50,0.54,0.90,0.84);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);
  legend->SetFillStyle(0);

  histMRAllBkg->SetLineColor(kRed);
  histMRAllBkg->GetYaxis()->SetTitle("Number of Events");
  histMRAllBkg->GetYaxis()->SetTitleOffset(1.2);
  histMRData->SetMarkerStyle(8);

  legend->AddEntry(histMRAllBkg, "All Backgrounds", "L");
  legend->AddEntry(histMRData, "Data", "L");

  histMRAllBkg->Add(histMRQCD, 1.0);

  histMRAllBkg->Draw("hist");
  histMRData->Draw("PE same");

  legend->Draw();
  cv->SetLogy();
  cv->SaveAs(Form("MR_QCD_%s.pdf",Label.c_str()));


  //////

  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.50,0.54,0.90,0.84);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);
  legend->SetFillStyle(0);

  histRsqAllBkg->SetLineColor(kRed);
  histRsqAllBkg->GetYaxis()->SetTitle("Number of Events");
  histRsqAllBkg->GetYaxis()->SetTitleOffset(1.2);
  histRsqData->SetMarkerStyle(8);

  legend->AddEntry(histRsqAllBkg, "All Backgrounds", "L");
  legend->AddEntry(histRsqData, "Data", "L");

  histRsqAllBkg->Add(histRsqQCD, 1.0);

  histRsqAllBkg->Draw("hist");
  histRsqData->Draw("PE same");

  legend->Draw();
  cv->SetLogy();
  cv->SaveAs(Form("Rsq_QCD_%s.pdf",Label.c_str()));

  //////////////////
  histRsqQCD->SetFillColor(kAzure+4);
  histRsqAllBkg->SetFillColor(kMagenta);

  stackRsq->Add(histRsqAllBkg);
  stackRsq->Add(histRsqQCD);
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.50,0.54,0.90,0.84);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);
  legend->SetFillStyle(0);
  stackRsq->Draw();
  stackRsq->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackRsq->GetHists()->At(0)))->GetXaxis()->GetTitle());
  stackRsq->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackRsq->GetHists()->At(0)))->GetYaxis()->GetTitle());
  stackRsq->Draw();
  histRsqData->Draw("same PE");
  legend->Draw();
  cv->SetLogy();
  cv->SaveAs(Form("RsqStack_QCD_%s.pdf",Label.c_str()));
 
   //--------------------------------------------------------------------------------------------------------------
  // Output
  //==============================================================================================================
  TFile *file = TFile::Open(("RazorPlots"+Label+".root").c_str(), "RECREATE");
  file->cd();

  for(int i=0; i<int(inputfiles.size()); i++) {
    file->WriteTObject(histMRRsq[i], Form("histMRRsq_%s",processLabels[i].c_str()), "WriteDelete");
    histUnrolled[i]->Write();  
  }
  
  stackUnrolled->Write();
 }

//.........这里部分代码省略.........
  histname[0]="weight";
  histname[1]="HT";
  histname[2]="MHT";
  histname[3]="NJet";
  histname[4]="NBtag";
  histname[5]="MuonPt";
  histname[6]="MtW";

  for(map<int , string >::iterator itt=WJtype.begin(); itt!=WJtype.end();itt++){        // loop over different event types

    cdtoitt = file->mkdir((itt->second).c_str());
    cdtoitt->cd();

    for(map<int , string >::iterator it=cutname.begin(); it!=cutname.end();it++){   // loop over different cutnames

      cdtoit =  cdtoitt->mkdir((it->second).c_str());
      cdtoit->cd();

      for(int j=0; j<histname.size(); j++){                                        // loop over different histograms

        for(int i=0; i<wjnHT ; i++){                                                  // loop over different HT bins

          //cout << "================================" << endl;
          //cout << "HT#: " <<i << ", WJtype: " << itt->second << ", cutname: " << it->second << ", hist#: " << j << endl;  
          sprintf(tempname,"%s/%s/%s_%s_%s",(itt->second).c_str(),(it->second).c_str(),(histname[j]).c_str(),(it->second).c_str(),(itt->second).c_str());
          temphist = (TH1D *) WJ_inputfilevec.at(i)->Get(tempname)->Clone();
          temphist->Scale(WJ_scalevec[i]);
          temphist->SetFillColor(i+2);
          tempstack->Add(temphist);

        }//end of loop over HTbins 1..7

        sprintf(tempname,"%s_%s_%s",histname[j].c_str(),(it->second).c_str(),(itt->second).c_str());
        tempstack->Write(tempname);

        delete tempstack;
        tempstack = new THStack("stack","Binned Sample Stack");

      }//end of loop over histograms

    }//end of loop over cutnames 

  }//end of loop over event types

  file->Close();
  printf("WJ main histograms stacked \n ");






// .....................................................................................................................................................//
// TTbar Section
// .....................................................................................................................................................//

  //build a vector of scale factors
  //first load the cross sections into a vector
  vector<double> TTbar_xs_vec;

  TTbar_xs_vec.push_back(806.1); // 

  const int ttbarnHT = (int) TTbar_xs_vec.size();   // Total number of HT bin samples

  for(int i=1; i<=ttbarnHT ; i++){
    if(i==1)sprintf(tempname,"../Results/results_TTbar_.root");

//.........这里部分代码省略.........
  //  type[8]="photon";
//    type[9]="H";
   type[10]="TTbar";
   type[11]="TTSingLep";
   type[12]="TTdiLep";
   type[13]="TThadronic";



  //KH
  histname[0]="weight";
  histname[1]="METAsys";
  histname[2]="MET";
  histname[3]="NJet";
  histname[4]="j1Pt";
  histname[5]="Jet1Eta";
  histname[6]="Jet1Phi";
histname[7]="j2Pt";
histname[8]="Jet2Eta";
histname[9]="Jet2Phi";
histname[10]="j3Pt";
histname[11]="Jet3Eta";
histname[12]="Jet3Phi";
histname[13]="DelPhij1j2";
histname[14]="NLep";
histname[15]="NElec";
histname[16]="NMuon";
histname[17]="NTau";
  ///end of initialization of the maps

  //build a vector of scale factors
  //first load the cross sections into a vector
  xs_vec.push_back(530.89358);
  xs_vec.push_back(42.55351);
  xs_vec.push_back(4.48209);
  xs_vec.push_back(0.52795);
  xs_vec.push_back(0.05449);

  double numberofevents =0;
  const int ttnHT = 5;   // Total number of HT bin samples
  const int nHist = 18; // Number of histograms in each TDirectory


  for(int i=1; i<=ttnHT ; i++){
    sprintf(tempname,"../Results/results_PhaseI_TT_14TEV_HT%d_NoPileUp_00.root",i);
    file = new TFile(tempname, "R");
    sprintf(tempname,"allEvents/nocut/MET_nocut_allEvents");
    tempvalue = (luminosity*xs_vec[i-1])/((* (TH1D* ) file->Get(tempname)).GetEntries());
    scalevec.push_back(tempvalue);
  }//end of loop over HTbins 
  std::cout << "normalization scale factor determination done" << std::endl;
for(int i=1; i<=ttnHT; i++){
sprintf(tempname,"../Results/results_PhaseI_TT_14TEV_HT%d_NoPileUp_00.root",i);
inputfilevec.push_back(TFile::Open(tempname,"R"));
}

tempstack = new THStack("stack","Binned Sample Stack");
file = new TFile("stack.root","RECREATE");
 for(map<int , string >::iterator itt=type.begin(); itt!=type.end();itt++){        // loop over different event types
    cdtoitt = file->mkdir((itt->second).c_str());
    cdtoitt->cd();
    for(map<int , string >::iterator it=cutname.begin(); it!=cutname.end();it++){   // loop over different cutnames
      cdtoit =  cdtoitt->mkdir((it->second).c_str());
      cdtoit->cd();
      for(int j=0; j<histname.size(); j++){                                        // loop over different histograms
        for(int i=0; i<ttnHT ; i++){                                                  // loop over different HT bins
sprintf(tempname,"%s/%s/%s_%s_%s",(itt->second).c_str(),(it->second).c_str(),(histname[j]).c_str(),(it->second).c_str(),(itt->second).c_str());
temphist = (TH1D *) inputfilevec.at(i)->Get(tempname)->Clone();
temphist->Scale(scalevec[i]);
if(histname[j]=="MET"){numberofevents+=(double)temphist->GetSumOfWeights();} //all the histograms in one directory have the same number of events
//if(histname[j]=="MET"){cout << " temphist->GetSumOfWeights() " << temphist->GetSumOfWeights() << endl;}
/*if(i==0){
cout << "" << endl;
cout << "type: " << (itt->second).c_str() << ",  cutname: " << (it->second).c_str()<< ", histname: " << histname[j].c_str() << ", bin#: " << i << endl;
cout << "temphist->GetEntries():  " << temphist->GetEntries() << endl;
cout << "temphist->GetSumOfWeights(): " << temphist->GetSumOfWeights() << endl;
cout << " ===============================================================" << endl;
}
*/
temphist->SetFillColor(i+2);
          tempstack->Add(temphist);
               }//end of loop over HTbins 1..7
if(histname[j]=="MET"){
cout << " +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
cout << "type: " << (itt->second).c_str() << ",  cutname: " << (it->second).c_str() << endl;
cout << "Number of events:  " << numberofevents << endl;
cout << " +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
}
numberofevents=0;
        sprintf(tempname,"%s_%s_%s",histname[j].c_str(),(it->second).c_str(),(itt->second).c_str());
        tempstack->Write(tempname);
        delete tempstack;
        tempstack = new THStack("stack","Binned Sample Stack");
      }//end of loop over histograms
    }//end of loop over cutnames
  }//end of loop over event types
file->Close();


}//end of the constructor

//.........这里部分代码省略.........
  // local histograms
  TH1D* h_all;
  THStack *hs = new THStack("hs","Stacked 1D histograms");

  // for local debugging
  TH1D* h_deno[nfiles];
  int COLORCODE[nfiles];

  COLORCODE[0] = kRed-7;
  COLORCODE[1] = kRed-10;
  COLORCODE[2] = kRed-6;
  COLORCODE[3] = kMagenta-2;
  COLORCODE[4] = kMagenta-6;
  COLORCODE[5] = kBlue-7;
  COLORCODE[6] = kBlue-9;
  for(int i=7; i<nfiles;i++)
    COLORCODE[i] = kGreen-6;

  cout << "opening " << myMCFiles[0].filename << endl;
  TFile *f1 = TFile::Open(myMCFiles[0].filename.data());

  TH1D* h_template = (TH1D*)(f1->Get(Form("%s",histoName.data())));
  h_template->Reset();

  h_all = (TH1D*)h_template->Clone(Form("%s_all",histoName.data()));

  h_all -> Reset();

  //-------------------------------------------------------------------------------
  // combine
  //-------------------------------------------------------------------------------

  

  for(int ifile=nfiles-1; ifile>=0; ifile--){

    cout << "File " << ifile << endl;
    TFile *f_temp = TFile::Open(myMCFiles[ifile].filename.data());
    cout << "opening " << myMCFiles[ifile].filename << endl;


    h_deno[ifile] = (TH1D*)(f_temp->Get(Form("%s",histoName.data())));
    h_deno[ifile] -> SetName(Form("h_deno_%d",ifile));
    h_deno[ifile] -> Rebin(rebin);
    h_deno[ifile] -> SetLineColor(COLORCODE[ifile]);
    h_deno[ifile] -> SetFillColor(COLORCODE[ifile]);
    h_deno[ifile] -> SetFillStyle(1001);
    h_deno[ifile] -> SetMarkerColor(COLORCODE[ifile]);
    h_deno[ifile] -> Sumw2();
    double weight = myMCFiles[ifile].scaleFactor;
    h_deno[ifile] -> Scale(weight);
    if(ifile==nfiles-1)
      {
	h_all    -> Rebin(rebin);
	h_all    -> Sumw2();
      }

    h_all    -> Add(h_deno[ifile]);
    hs       -> Add(h_deno[ifile]);

    // to be used with TEfficiency methods
    cout << h_deno[ifile]->GetEntries() << endl;
  } // end of loop over files

  
  TCanvas* c1 = new TCanvas("c1",inputFile.data(),0,0,500,500);
  h_all->SetMarkerSize(1);
  h_all->SetMarkerStyle(24);
  h_all->Draw("e");
  h_all->SetXTitle(xtitle.data());
  h_all->GetXaxis()->SetRangeUser(xmin,xmax);
  for(unsigned int ifile=0; ifile < nfiles; ifile++){
    h_deno[ifile]->Draw("hist,same");
  }
  
  std::string remword  =".txt";
  size_t pos  = inputFile.find(remword);

  if(pos!= std::string::npos)
    inputFile.swap(inputFile.erase(pos,remword.length()));

  std::string command = "recreate";

  if(update)command ="update";

  TCanvas* c2 = new TCanvas("c2",inputFile.data(),500,0,500,500);
  hs->Draw("hist");
  hs->GetXaxis()->SetRangeUser(xmin,xmax);
  hs->GetXaxis()->SetTitle(xtitle.data());
  hs->Draw("hist");

  

  TFile* outFile = new TFile(Form("combined_%s.root",inputFile.data()),command.data());               
  h_all->Write();
  hs->Write();
  outFile->Close();     


}

void Interpolate(const TString& trigger="INEL")
{
  if (gSystem->Getenv("FWD"))
    fwd = gSystem->Getenv("FWD");
  else 
    fwd = gSystem->ExpandPathName("$ALICE_PHYSICS/PWGLF/FORWARD/analysis2");
  gROOT->SetMacroPath(Form("%s/dndeta:%s", gROOT->GetMacroPath(),fwd));

  if (!gROOT->GetClass("Drawer"))  gROOT->LoadMacro("Drawer.C+");

  TH1* h0900 = GetOne( 900, trigger);
  TH1* h2760 = GetOne(2760, trigger);
  TH1* h7000 = GetOne(7000, trigger);
  TH1* h8000 = GetOne(8000, trigger);
  Info("","900: %p 2760: %p 7000: %p 8000: %p", 
       h0900, h2760, h7000, h8000);
  Double_t e8000 = (trigger.EqualTo("INEL") ? 0.852 : 0.93);
  h8000->Scale(e8000);

  TFile* out = TFile::Open("trends.root", "RECREATE");
  THStack* sOrig = new THStack("orig", Form("pp - %s", trigger.Data()));
  sOrig->Add(h8000);
  sOrig->Add(h7000);
  sOrig->Add(h2760);
  sOrig->Add(h0900);

  TCanvas* cOrig = new TCanvas("cOrig", "Original", 1200, 1200);
  cOrig->SetTopMargin(0.01);
  cOrig->SetRightMargin(0.01);
  sOrig->Draw("nostack");
  sOrig->GetHistogram()->SetYTitle("1/#it{N} d#it{N}_{ch}/d#it{#eta}");
  sOrig->GetHistogram()->SetXTitle("#it{#eta}");
  sOrig->DrawClone("nostack");
  sOrig->Write();
  
  TLegend* l = cOrig->BuildLegend(.35, .2, .55, .6, "#sqrt{s}");
  l->SetFillColor(0);
  l->SetFillStyle(0);
  l->SetBorderSize(0);

  cOrig->Modified();
  cOrig->Update();
  cOrig->cd();
  cOrig->Write();
  Info("", "Wrote original"); 

  TCanvas* cG = new TCanvas("cG", "one", 1200, 1200);
  cG->SetTopMargin(0.01);
  cG->SetRightMargin(0.01);
  
  Info("","Creating tuple");
  TNtuple* tuple = new TNtuple("tuple", "Tuple", 
			       "eta:deta:"
			       "v0900:e0900:v2760:e2760:"
			       "v7000:e7000:v8000:e8000");
  TMultiGraph* mg = new TMultiGraph;
  Int_t n = h0900->GetNbinsX();
  Info("","Loop over bins %d", n);

  for (Int_t i = 1; i <= n; i++) {
    Info("", "Getting one bin %d,%p,%p,%p,%p,%p,%p",
	 i, h0900,h2760,h7000,h8000,mg,tuple);
    OneBin(i, h0900, h2760, h7000, h8000, mg, tuple);
  }
  mg->Draw("alp");

  cG->Modified();
  cG->Update();
  cG->cd();

  TPrincipal* p =tuple->Principal("v0900:v2760:v7000:v8000","eta<0", "npdhc");
}