C++ TH1F::Print方法代码示例

//testing how the TH1::Integral() values change when a hist is scaled
void scale()
{
	gSystem->Load("../RootTools/CommonTools.C.so");
	gStyle->SetOptStat("nemouri");
	
	TH1F *h = new TH1F("H1","Test Integral",100,0,10);
	h->Fill (1);
	h->Fill (1);
	h->Fill (1);
	h->Fill (1);
	//TH1F *h2 = new TH1F("H2","Test Integral copy",100,0,1);
	//h->FillRandom("gaus",10000);
	//h2->FillRandom("gaus",10000);

//	std::cout << "nbins, Entries, Integral = " << h->GetNbinsX() << ", " 
//				<< h->GetEntries() << ", " << h->Integral("width") << std::endl;
	//h = (TH1F*) MakeVariableBinHist(h, 0, 0.2, 0.6, 0.8, 1, 0.1, 0.2, 0.2, 0.2);
	//std::cout << "nbins, Entries, Integral = " << h->GetNbinsX() << ", " 
	//			<< h->GetEntries() << ", " << h->Integral("width") << std::endl;

	new TCanvas();
	h->Draw();
	//new TCanvas();
	//h2->Draw();

	//h->Scale(1,"width"); //this replaces the loop below
	/*
	for (int bin = 1; bin <= h->GetNbinsX(); ++bin)
	{
		const double val = h->GetBinContent(bin);
		const double err = h->GetBinError(bin);
		const double width = h->GetBinWidth(bin);
		h->SetBinContent(bin,val/width);
		h->SetBinError(bin,err/width);
	}
	*/
	//h->DrawClone("E");
	h->Print();
	
	std::cout << "nbins, Entries, Integral = " << h->GetNbinsX() << ", " 
				<< h->GetEntries() << ", " << h->Integral("width") << std::endl;


	//h->Scale(10);
	h->Scale(10/(double)h->Integral());
	//h_copy->Scale(h2_copy->Integral()/(double) h2->Integral());
	//h_copy->Draw();

	std::cout << "after scaling to " << 10 << std::endl;
	h->Print();
	std::cout << "nbins, Entries, Integral = " << h->GetNbinsX() << ", " 
				<< h->GetEntries() << ", " << h->Integral("width") << std::endl;
	
};

void real_emission_jackknife(TString dir="real_aajj_splitted_for_jackknife", TString nameroot="outphoton_aajj_effunf_mu1", int N=20){

  TFile *outfile = new TFile(Form("%s/%s_jackknifed.root",dir.Data(),nameroot.Data()),"recreate");
  outfile->mkdir("effunf");

  vector<TFile*> files;
  for (int i=1; i<=N; i++){
    files.push_back(new TFile(Form("%s/%s_splitted_%d.root",dir.Data(),nameroot.Data(),i),"read"));
  }

  for (std::vector<TString>::const_iterator it = diffvariables_list.begin(); it!=diffvariables_list.end(); it++){
    for (int reg=0; reg<3; reg++)
      {
	TString diffvar = *it;
	vector<TH1F*> histos;
	TString thisname = Form("htruth_%s_%d",diffvar.Data(),reg);
	TString thisnamedir = Form("effunf/%s",thisname.Data());
	for (uint i=0; i<files.size(); i++){
	  histos.push_back((TH1F*)(files.at(i)->Get(thisnamedir.Data())));
	}
	TH1F *out = JackknifammiStiHistos(histos);
	out->SetName(thisname.Data());
	out->Print();
	outfile->cd("effunf");
	out->Write();
#ifdef DEBUG
	return;
#endif
      } 
  }
  
  outfile->Close();
};

void fitTF1Self() {
	TFile *f = TFile::Open("plainHisto.root","read");
	TH1F *h = (TH1F*)f->Get("mbbReg;1");
	h->Print();
	TCanvas *c = new TCanvas("c1","c1");
	double params[7];
	fitTF1(c,*h,110,140,2.5,params,kRed);
}

int execDrawString()
{
   cernbuild(1);
#ifdef ClingWorkAroundMissingDynamicScope
   TTree *T; gFile->GetObject("T",T);
#endif
   T->Draw("Division");
#ifdef ClingWorkAroundMissingDynamicScope
   TH1F *htemp = (TH1F*)gROOT->FindObject("htemp");
#endif
   htemp->Print("all");
#ifdef ClingWorkAroundErracticValuePrinter
   printf("(int)0\n");
#endif
   return 0;
}

void efficiency( TString denfile="60-120.root",
                 TString numfile="orig.root",
                 TString outfile="res.root" ) {

  gROOT ->Reset();
  
  TFile * nume = new TFile(numfile);
  TFile * deno = new TFile(denfile);
  
  nume->cd("/");
  deno->cd("/");

  TString theName = "d02-x01-y01";

  const TH1F *myNume;
  nume->GetObject(theName,myNume);
  myNume->Sumw2();
  myNume->Print("all");

  const TH1F *myDeno;
  deno->GetObject(theName,myDeno);
  myDeno->Sumw2();
  myDeno->Print("all");

  TH1F *Acce = (TH1F*)myNume->Clone("efficiency");
  Acce->Reset();
  Double_t f1(1.),f2(1.);
  Option_t* opt("b");
  Acce->Divide(myNume,myDeno,f1,f2,opt);
  Acce->Print("all");

  TFile * out = new TFile(outfile,"NEW");
  Acce->Write();
  out->Close();

}

void merge_pbpb_pp_HLT(){
  
  TH1::SetDefaultSumw2();
  
  //TFile *fpbpb1 = TFile::Open("/mnt/hadoop/cms/store/user/rkunnawa/rootfiles/PbPb/2011/data/ntuple_2011_pbpbJet80_v2.root");
  //TFile *fpbpb2 = TFile::Open("/mnt/hadoop/cms/store/user/rkunnawa/rootfiles/PbPb/2011/data/ntuple_2011_pbpbJet65_v2.root");
  //TFile *fpbpb3 = TFile::Open("/mnt/hadoop/cms/store/user/rkunnawa/rootfiles/PbPb/2011/data/ntuple_2011_pbpbJet55_v2.root");
  
  //TFile *fpp1 = TFile::Open("/mnt/hadoop/cms/store/user/rkunnawa/rootfiles/PP/2013/data/ntuple_2013_ppJet80.root");
  //TFile *fpp2 = TFile::Open("/mnt/hadoop/cms/store/user/rkunnawa/rootfiles/PP/2013/data/ntuple_2013_ppJet40.root");
  
  TFile *fpp1_v2 = TFile::Open("/mnt/hadoop/cms/store/user/rkunnawa/rootfiles/PP/2013/data/ntuple_2013_JEC_applied_ppJet80_v2.root");
  TFile *fpp2_v2 = TFile::Open("/mnt/hadoop/cms/store/user/rkunnawa/rootfiles/PP/2013/data/ntuple_2013_JEC_applied_ppJet40_v2.root");

  /*
  TTree *jetpbpb1_v2 = (TTree*)fpbpb1->Get("jetR3");
  TTree *jetpbpb2_v2 = (TTree*)fpbpb2->Get("jetR3");
  TTree *jetpbpb3_v2 = (TTree*)fpbpb3->Get("jetR3");

  TTree *evtpbpb1_v2 = (TTree*)fpbpb1->Get("evt");
  TTree *evtpbpb2_v2 = (TTree*)fpbpb2->Get("evt");
  TTree *evtpbpb3_v2 = (TTree*)fpbpb3->Get("evt");  

  jetpbpb1_v2->AddFriend(evtpbpb1_v2);
  jetpbpb2_v2->AddFriend(evtpbpb2_v2);
  jetpbpb3_v2->AddFriend(evtpbpb3_v2);
  */
  //TTree *jetpp1 = (TTree*)fpp1->Get("ntjet");
  // TTree *jetpp2 = (TTree*)fpp2->Get("ntjet");

  TTree *jetpp1_v2 = (TTree*)fpp1_v2->Get("jetR3");
  TTree *jetpp2_v2 = (TTree*)fpp2_v2->Get("jetR3");

  TTree *evtpp1_v2 = (TTree*)fpp1_v2->Get("evt");
  TTree *evtpp2_v2 = (TTree*)fpp2_v2->Get("evt");

  jetpp1_v2->AddFriend(evtpp1_v2);
  jetpp2_v2->AddFriend(evtpp2_v2);

  //TCut pbpb3 = "abs(eta)<2&&jet55&&!jet65&&!jet80&&chMax/pt>0.01";
  TCut pp3 = "abs(eta)<2&&jet40&&!jet60&&!jet80&&chMax/pt>0.01";
  
  //TH1F *hpbpb1 = new TH1F("hpbpb1","",30,0,300);
  //TH1F *hpbpb2 = new TH1F("hpbpb2","",30,0,300);
  //TH1F *hpbpb3 = new TH1F("hpbpb3","",30,0,300);
  //TH1F *hpbpbComb = new TH1F("hpbpbComb","",30,0,300);
  
  TH1F *hpp1 = new TH1F("hpp1","",nbins_jetPtBin,boundaries_jetPtBin);
  TH1F *hpp2 = new TH1F("hpp2","",nbins_jetPtBin,boundaries_jetPtBin);
  TH1F *hpp3 = new TH1F("hpp3","",nbins_jetPtBin,boundaries_jetPtBin);
  TH1F *hppComb = new TH1F("hppComb","",nbins_jetPtBin,boundaries_jetPtBin);
  
  //get the prescl factor information. 
  //Float_t presclpbpb3 = (Float_t)jetpbpb1_v2->GetEntries("jet80")/jetpbpb1_v2->GetEntries("jet55&&jet80");
  //cout<<"pbpb prescl3 = "<<presclpbpb3<<endl;//1.99871
  Float_t presclpp3 = (Float_t)jetpp1_v2->GetEntries("jet80")/jetpp1_v2->GetEntries("jet40&&jet80");
  cout<<"pp prescl3 = "<<presclpp3<<endl; //9.24968
  /*
  jetpbpb1_v2->Project("hpbpb1","pt","abs(eta)<2&&jet80&&chMax/pt>0.01");
  hpbpb1->Print("base");
  divideBinWidth(hpbpb1);

  jetpbpb2_v2->Project("hpbpb2","pt","abs(eta)<2&&jet65&&!jet80&&chMax/pt>0.01");
  hpbpb2->Print("base");
  divideBinWidth(hpbpb2);

  jetpbpb3_v2->Project("hpbpb3","pt","1.9987"*pbpb3);
  hpbpb3->Print("base");
  divideBinWidth(hpbpb3);
  */
  jetpp1_v2->Project("hpp1","pt","abs(eta)<2&&jet80&&chMax/pt>0.01");
  hpp1->Print("base");
  //divideBinWidth(hpp1);
  //hpp1->Scale(1./3.083e11);
  //hpp1->Scale(1./4);

  jetpp2_v2->Project("hpp2","pt","abs(eta)<2&&jet60&&!jet80&&chMax/pt>0.01");
  hpp2->Print("base");
  //divideBinWidth(hpp2);
  //hpp2->Scale(1./3.083e11);
  //hpp2->Scale(1./4);

  jetpp2_v2->Project("hpp3","pt","9.25038"*pp3);
  //9.25038 - ak5
  //9.24955 - ak4
  //9.24968 - ak3
  hpp3->Print("base");
  //divideBinWidth(hpp3);
  //hpp3->Scale(1./3.083e11);
  //hpp3->Scale(1./4);

  //scale the PbPb histograms before adding them
  //we have to scale them according to the lumi of the Jet80 file. 
  // HLT file  |   Lumi
  // HLT_80    |   150 mub-1
  // HLT_65    |   12.1 mub-1
  // HLT_55    |   0.38 mub-1
  // 
  // therefore scale for HLT_55 = 150/0.38 = 394.73684
  // scale for HLT_65 = 150/12.1 = 12.3967
//.........这里部分代码省略.........

void test(int numbersigmas = 0, Bool_t debugtest = true)
{

    using namespace RooFit;
    using namespace std;

    TCanvas *canvas = new TCanvas("canvas","canvas",900,100,500,500);

    gSystem->Load("libRooFit");
    gSystem->AddIncludePath("-I$ROOFITSYS/include");


    float ptbinsarray[] = {20.,40.,60.,80.,100.,120.,200.,600.};

    std::vector<float> ptbins(ptbinsarray,ptbinsarray+sizeof(ptbinsarray)/sizeof(ptbinsarray[0]));

    std::vector<std::vector<float> > allbins;
    allbins.push_back(ptbins);

    std::vector<TString> VarString;
    VarString.push_back("VsPt");

    std::vector<TString> HistoNameString;
    HistoNameString.push_back("ptbin");

    std::vector<TString> GraphXTitleString;
    GraphXTitleString.push_back("p_{t} (GeV)");

    std::vector<TString> SideBandDefinitions;
    SideBandDefinitions.push_back("SideBand5_10");



    //   ------------FOR TESTING----------------
    unsigned int sidebandloopmax = 1;//5_10, 5_20, ...
    unsigned int templatevarsloopmax = 1;//sinin with conv safe veto, sinin, ch isol
    unsigned int binsloopmax = 1;//pt, eta, phi, pu

    //----------Open .root Templates

    //All the Jet Templates
    TFile *histojetfile = TFile::Open("/afs/cern.ch/user/c/ciperez/CMSSW_7_4_5/src/TemplateHistosJetCheckBinsEndcapsLoose.root");
    //TFile *histojetdenfile = TFile::Open("/afs/cern.ch/user/c/ciperez/CMSSW_7_4_5/src/Denominator_FREndCaps.root");
// TFile *histojetdatafile = TFile::Open("/afs/cern.ch/user/c/ciperez/CMSSW_7_4_5/src/Num_Templates.root");

    //Real Photon Templates
    TFile *historealmcfile = TFile::Open("/afs/cern.ch/user/c/ciperez/CMSSW_7_4_5/src/RealPhotonTemplatesEndCaps.root");

    //--- Write NEW .root Historams for Fake Rate
    TFile *FRhistosfile = new TFile("FakeRatePlotsCheckBins.root","recreate");

    //loop on error systematics

    for(unsigned int m = 0; m<templatevarsloopmax; m++) {
        for(unsigned int l = 0; l<binsloopmax; l++) {


            TMultiGraph *mg = new TMultiGraph();
            TLegend *legendAllGraphs = new TLegend(0.37,0.59,0.57,0.79);
            legendAllGraphs->SetTextSize(0.02);
            legendAllGraphs->SetFillColor(kWhite);
            legendAllGraphs->SetLineColor(kWhite);


            std::vector<float> fakeratevalues;
            std::vector<float> fakerateptvalues;
            std::vector<float> fakerateerrorvalues;

            for(unsigned int k = 0; k<allbins[0].size()-1; k++) {
// for(unsigned int k = 0;k<1;k++){
                float binlow = allbins[0][k];
                float binmax = allbins[0][k+1];

                TString binstring = TString::Format("%4.2f_%4.2f",binlow,binmax);
                binstring.ReplaceAll(".00","");
                binstring.ReplaceAll("-","m");
                binstring.ReplaceAll(".","p");
                binstring.ReplaceAll("10000","Inf");
                cout<<binstring.Data()<<endl;

                //Histograms for templates
                //Get histograms from each of the histojetfiles declared earlier

                //Numerator Fakes - FakePhotonNumEndCaps.root - FakePhoton_num2040...
                TH1F *h1 = (TH1F*)histojetfile->Get(("histoSininWithPixelSeedFakeJetptbin"+binstring).Data());// .Data() changes to char*
                h1->Print();

                //Numerator Real Photons - RealPhotonTemplatesEndCaps.root - EndCapsMCReal_20_40
                TH1F *h2 = (TH1F*)historealmcfile->Get(("EndCapsMCReal_"+binstring).Data());
                h2->Print();

                //Numerator Templates - Num_Templates.root - num2040...
                TH1F *hData = (TH1F*)histojetfile->Get(("histoSininWithPixelSeedDataJetptbin"+binstring).Data());
                hData->Print();

                //Denominator Templates - Denominator_FREndCaps.root -FakePhoton_den2040...
                TH1F *hnum = (TH1F*)histojetfile->Get(("histoSininWithPixelSeedTightAndFakeJetptbin"+binstring).Data());
                hnum->Print();

                //avoiding 0 entries in the histograms
//.........这里部分代码省略.........

void fitSignal(){

std::string shapes_file = "mlfit.root";
std::string data_file = "param_ws.root";
std::string channel = "ch1";
     
     TFile *dfile = TFile::Open(data_file.c_str());
     TFile *sfile = TFile::Open(shapes_file.c_str());

     TH1F *bkg   = (TH1F*)sfile->Get(Form("shapes_fit_b/%s/total",channel.c_str()));
     TH1F *data  = (TH1F*)dfile->Get("SR_data");		// TH1 for data 
     TH1F *signal= (TH1F*)dfile->Get("SR_signal");		// TH1 for signal 
     TH2F *covar = (TH2F*)sfile->Get(Form("shapes_fit_b/%s/total_covar",channel.c_str()));

     // bkg and covariance defined as pdf / GeV, so scale by bin widhts 
     //

     int nbins = data->GetNbinsX();

     if (!isTH1Input){
      for (int b=1;b<=nbins;b++){
       double bw = bkg->GetBinWidth(b);
       bkg->SetBinContent(b,bkg->GetBinContent(b)*bw);
       for (int j=1;j<=nbins;j++){
       	covar->SetBinContent(b,j,covar->GetBinContent(b,j)*bw*bw);
       }
      }
     }

     RooArgList xlist_;
     RooArgList olist_;
     RooArgList mu_;

     bkg->Print() ;
     covar->Print() ; 
     signal->Print() ;
     data->Print() ;  

     // Make a dataset (simultaneous)
     RooCategory sampleType("bin_number","Bin Number");
     RooRealVar  observation("observed","Observed Events bin",1);

     // You have to define the samples types first!, because RooFit suuuuucks!
     for (int b=1;b<=nbins;b++){
        sampleType.defineType(Form("%d",b-1),b-1);
        sampleType.setIndex(b-1);
     }

     RooArgSet   obsargset(observation,sampleType);
     RooDataSet obsdata("combinedData","Data in all Bins",obsargset);
     //obsdata.add(RooArgSet(observation,sampleType));

     for (int b=1;b<=nbins;b++){
        observation.setVal(data->GetBinContent(b));
        sampleType.setIndex(b-1);
        std::cout << sampleType.getLabel() << ", " << sampleType.getIndex() << std::endl;
        //RooArgSet localset(observation,sampleType);
   	//obsdata.add(localset);
        obsdata.add(RooArgSet(observation,sampleType));
	std::cout << " Observed at " << b << ", " << observation.getVal() << std::endl;
     }

     // make a constraint term for the background, and a RooRealVar for bkg 
     for (int b=1;b<=nbins;b++){
	double bkgy = (double)bkg->GetBinContent(b);
	RooRealVar *mean_ = new RooRealVar(Form("exp_bin_%d_In",b),Form("expected bin %d",b),bkgy); 
	mean_->setConstant(true);
	RooRealVar *x_ = new RooRealVar(Form("exp_bin_%d",b),Form("bkg bin %d",b),bkgy,0.2*bkgy,bkgy*4);
	std::cout << " Exp background At " << b << ", " << x_->getVal() << std::endl;
	xlist_.add(*x_);
	mu_.add(*mean_);
     }      

     // constraint PDF for background
     // Convert TH2 -> TMatrix 
     TMatrixDSym Tcovar(nbins);
     for (int i=0;i<nbins;i++){
      for (int j=0;j<nbins;j++){
	//if (i==j)Tcovar[i][j] = covar->GetBinContent(i+1,j+1);
	//else Tcovar[i][j] = 0;
	Tcovar[i][j] = covar->GetBinContent(i+1,j+1);
      }
     }
     std::cout<< "Made Covariance" << std::endl;
     RooMultiVarGaussian constraint_pdf("constraint_pdf","Constraint for background pdf",xlist_,mu_,Tcovar);
     std::cout<< "Made Covariance Gauss" << std::endl;
     
     // Make the signal component 
     RooRealVar r("r","r",1,-5,5);
     RooArgList signals_;
     for (int b=1;b<=nbins;b++) {
	RooProduct *sigF = new RooProduct(Form("signal_%d",b),"signal nominal",RooArgSet(r,RooFit::RooConst(signal->GetBinContent(b))));
	std::cout << " Signal At " << b << ", " << sigF->getVal() << std::endl;
	signals_.add(*sigF);
     }

     RooArgList plist_;
     RooArgList slist_;
 
     sampleType.setIndex(1); 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	hMeas_80->Scale(1./N_mb);
	hMeas_60->Scale(1./N_mb);
	
	TCanvas *cMerged = new TCanvas("cMerged","Merged PPb spectra",800,600);
	cMerged->SetLogy();

	hCombined->SetXTitle("Jet p_{T} [GeV/c]");
	hCombined->SetYTitle("1/N_mb d^2N/dp_t d eta");
	hCombined->SetMarkerColor(kBlack);
	hCombined->SetMarkerStyle(20);
	hCombined->SetAxisRange(30,300,"X");
	hCombined->Draw();
	hMeas_100->SetMarkerColor(kRed);
	hMeas_100->SetMarkerStyle(21);
	hMeas_100->Draw("same");
	hMeas_80->SetMarkerColor(kBlue);
	hMeas_80->SetMarkerStyle(22);
	hMeas_80->Draw("same");
	hMeas_60->SetMarkerColor(kGreen);
	hMeas_60->SetMarkerStyle(23);
	hMeas_60->Draw("same");
	
	TLegend *leg_PPb = myLegend(0.6,0.65,0.95,0.9);
	leg_PPb->SetTextSize(0.05);
	leg_PPb->AddEntry(hCombined,"Merged PPb spectra ","pl");
	leg_PPb->AddEntry(hMeas_100,"HLT_100 spectra","pl");
	leg_PPb->AddEntry(hMeas_80,"HLT_80 spectra","pl");
	leg_PPb->AddEntry(hMeas_60,"HLT_60 spectra","pl");
	leg_PPb->Draw();
	
	putCMSPrel(0.2,0.83,0.06);
	drawText("PPb AKPu3PF |eta|<1 |vz|<15",0.2,0.23,20);
	
	cMerged->SaveAs("pPb_merged.pdf","RECREATE");
	
	hMeas_100->Write();
	hCombined->Write();
	hMeas_80->Write();
	hMeas_60->Write();
	

	TFile *fYaxian = TFile::Open("AkPu3PFJetRpA.root");
	TH1F *Yaxian = (TH1F*)fYaxian->Get("DataJetWideBin;3");
	TH1F *test = (TH1F*)Yaxian->Clone("test");

	cout<<"hi"<<endl;
	//outfile->cd();
	Yaxian->Print("base");
	test->Print("base");
	Yaxian->Divide(hCombined);
	TCanvas *yaxian = new TCanvas("yaxian","",800,600);
	yaxian->Divide(2,1);
	yaxian->cd(1);
	Yaxian->SetTitle("ratio of Yaxian's measured pPb spectra to Mine");
	Yaxian->SetXTitle("Jet p_{T} [GeV/c]");
	Yaxian->SetMarkerColor(kBlack);
	Yaxian->SetMarkerStyle(23);
	Yaxian->Draw();
	
	yaxian->cd(2);
	yaxian->cd(2)->SetLogy();
	test->SetMarkerStyle(22);
	test->SetMarkerColor(kBlack);
	test->SetXTitle("Jet p_{T} [GeV/c]");
	test->SetYTitle("1/N_mb d^2N/dp_t d eta");
	test->Draw();
	hCombined->Draw("same");
	
	TLegend *leg = myLegend(0.6,0.65,0.95,0.9);
	leg->SetTextSize(0.05);
	leg->AddEntry(hCombined,"Merged PPb spectra ","pl");
	leg->AddEntry(test,"Yaxian's spectra","pl");
	leg->Draw();
	
	putCMSPrel(0.2,0.83,0.06);
	drawText("PPb AKPu3PF |eta|<1 |vz|<15",0.2,0.23,20);
	
	yaxian->SaveAs("Yaxian_Comparison_pPb_pt_spectra.root","RECREATE");
	
	outfile->cd();
	Yaxian->Write();
	test->Write();
	
	outfile->Write();
	outfile->Close();
	
	timer.Stop();
    float rtime  = timer.RealTime();
    float ctime  = timer.CpuTime();
    
    std::cout<<"\t"<<std::endl;
    std::cout<<Form("RealTime=%f seconds, CpuTime=%f seconds",rtime,ctime)<<std::endl;
    std::cout<<"\t"<<std::endl;
    std::cout<<"Good bye : " <<"\t"<<std::endl;
    
	//define the required histograms
	//static const Int_t nbins = 22;
	//static const Double_t bound[nbins+1] = {30.,40.,50.,60.,70.,80.,90.,100.,110.,120.,130.,140.,150.,160.,180.,200.,220.,260.,300.,350.,400.,450.,500.};
	
}

int main() {
  TH1::AddDirectory(0);
  ModTDRStyle();
  // string cfg;                                   // The configuration file

  // bool do_ratio = false;
  // bool do_logy = true;

  // vector<string> datacards;
  // string fitresult_file = "";
  // string parse_rule = "";

  // string output           = "";
  // string text1            = "";
  // string text2            = "";
  // bool postfit            = true;

  gSystem->Load("libHiggsAnalysisCombinedLimit.dylib");

  ch::CombineHarvester cmb;
  // cmb.SetVerbosity(2);
  // cmb.ParseDatacard(
  //     "output/cmshcg/summer2013/searches/hzz4l/125/hzz4l_4muS_8TeV_0.txt",
  //     "$MASS/$ANALYSIS_$CHANNEL_$ERA_$BINID.txt");
  cmb.ParseDatacard(
      "output/cmshcg/summer2013/searches/hzz4l/125/hzz4l_4muS_8TeV_1.txt",
      "$MASS/$ANALYSIS_$CHANNEL_$ERA_$BINID.txt");
  // cmb.ParseDatacard(
  //     "output/cmshcg/summer2013/searches/hzz4l/125/hzz4l_4muS_7TeV_0.txt",
  //     "$MASS/$ANALYSIS_$CHANNEL_$ERA_$BINID.txt");
  // cmb.ParseDatacard(
  //     "output/cmshcg/summer2013/searches/hzz4l/125/hzz4l_4muS_7TeV_1.txt",
  //     "$MASS/$ANALYSIS_$CHANNEL_$ERA_$BINID.txt");

  ch::SetStandardBinNames(cmb);

  RooFitResult fitresult = ch::OpenFromTFile<RooFitResult>(
      "output/cmshcg/summer2013/searches/hzz4l/125/mlfit.root:fit_s");
  auto fitparams = ch::ExtractFitParameters(fitresult);
  cmb.UpdateParameters(fitparams);

  cmb.process_rgx({".*SM"}, false);
  // cmb.PrintAll();

  TH1F sig_shape = cmb.cp().signals().GetShape();
  TH1F bkg_shape = cmb.cp().backgrounds().GetShape();
  bkg_shape.Print("range");
  std::cout << bkg_shape.Integral() << "\n";
  TH1F data = cmb.GetObservedShape();
  data.Print("range");
  int rebin = 1;
  sig_shape.Rebin(rebin);
  bkg_shape.Rebin(rebin);
  data.Rebin(rebin);

  sig_shape.Scale(data.GetBinWidth(1) / sig_shape.GetBinWidth(1));
  bkg_shape.Scale(data.GetBinWidth(1) / bkg_shape.GetBinWidth(1));

  TCanvas* canv = new TCanvas("hzz4l", "hzz4l");
  canv->cd();

  std::vector<TPad*> pads = OnePad();
  std::vector<TH1*> h = CreateAxisHists(1, &data);
  StandardAxes(h[0]->GetXaxis(), h[0]->GetYaxis(), "m_{4l}", "GeV");
  h[0]->Draw("axis");

  sig_shape.Add(&bkg_shape);
  sig_shape.SetLineColor(4);
  sig_shape.SetLineWidth(3);
  sig_shape.Draw("SAME HIST C");

  bkg_shape.SetLineColor(2);
  bkg_shape.SetLineWidth(3);
  bkg_shape.Draw("SAME HIST C");

  data.Draw("esamex0");

  FixTopRange(pads[0], GetPadYMax(pads[0]), 0.15);
  DrawCMSLogo(pads[0], "CMS", "Preliminary", 11, 0.045, 0.035, 1.2);
  DrawTitle(pads[0], "19.7 fb^{-1} (8 TeV)", 3);
  DrawTitle(pads[0], "H#rightarrowZZ", 1);
  canv->Update();
  pads[0]->RedrawAxis();
  pads[0]->GetFrame()->Draw();
  canv->SaveAs("hzz4l.pdf");

  return 0;
}

void plot_Validatin_DatavsMC(int radius = 4,
			     std::string coll = "PP",
			     std::string algo = "",
			     std::string jetType = "PF",
			     std::string DataFile = "PromptForestPP_DataJet80_ak4PF.root",
			     std::string MCFile = "PromptForestPP_MC_ak4PF.root")
{

  // get the data and MC histograms
  // these are simple histograms
  TFile * fData = TFile::Open(DataFile.c_str());
  TH1F * pt2overpt1_Data = (TH1F*)fData->Get("pt2overpt1");
  TH1F * hJetEta_Data = (TH1F*)fData->Get("hJetEta");
  TH1F * hJetPhi_Data = (TH1F*)fData->Get("hJetPhi");
  TH1F * hJetpT_Data = (TH1F*)fData->Get("hJetpT");
  hJetpT_Data->Print("base");
  TH1F * hAj_Data = (TH1F*)fData->Get("hAj");
  TH1F * hDeltaPhi_Data = (TH1F*)fData->Get("hDeltaPhi");

  TH1F * hevent = (TH1F*)fData->Get("hRunN_Vs_NJets");
  hevent->Print("base");
  TH1F * hJet80 = (TH1F*)fData->Get("hJet80");
  hJet80->Print("base");
  
  TFile * fMC = TFile::Open(MCFile.c_str());
  TH1F * pt2overpt1_MC = (TH1F*)fMC->Get("pt2overpt1");
  TH1F * hJetEta_MC = (TH1F*)fMC->Get("hJetEta");
  TH1F * hJetPhi_MC = (TH1F*)fMC->Get("hJetPhi");
  TH1F * hJetpT_MC = (TH1F*)fMC->Get("hJetpT");
  TH1F * hAj_MC = (TH1F*)fMC->Get("hAj");
  TH1F * hDeltaPhi_MC = (TH1F*)fMC->Get("hDeltaPhi");

  // get marta's histograms
  TFile * fMarta = TFile::Open("AnaResultsPFvsCaloJets4M.root");
  TList *lst = (TList*)fMarta->Get("anaPFvsCaloJet");
  TH3F * hpT3 = (TH3F*)lst->FindObject("fh3PtTrueEtaDeltaPt");
  int etabinlow = hpT3->GetYaxis()->FindBin(-2.00001);
  int etabinhigh = hpT3->GetYaxis()->FindBin(2.00001);
  TH1F * hPFCaloMatchpT = (TH1F*)hpT3->ProjectionX("hPFCaloMatchpT",etabinlow, etabinhigh);
  
  hPFCaloMatchpT->Print("base");
  
  // plot them on top of each other. 
  TCanvas * cAj = new TCanvas("Aj","",800,600);
  hAj_Data->SetMarkerStyle(20);
  hAj_Data->SetMarkerColor(kBlack);
  hAj_Data->SetAxisRange(0.0, 1.0, "X");
  hAj_Data->SetYTitle("Event Fraction");
  hAj_Data->SetXTitle("A_{j}");
  hAj_Data->SetTitle(" ");
  hAj_Data->DrawNormalized();
  
  hAj_MC->SetMarkerStyle(25);
  hAj_MC->SetMarkerColor(kRed);
  hAj_MC->DrawNormalized("same");

  putCMSPrel();
  //putPPLumi();
  
  TLegend * lAj = myLegend(0.4,0.6,0.7,0.9);
  lAj->AddEntry("","p_{T}^{lead}   > 100 GeV/c","");
  lAj->AddEntry("","p_{T}^{sublead}> 40  GeV/c","");
  lAj->AddEntry(hAj_Data,"Data (Jet 80 trigger)","pl");
  lAj->AddEntry(hAj_MC,"MC (pthat 80)","pl");
  lAj->Draw();

  cAj->SaveAs(Form("Aj_datavsMC_%s_ak%s%d%s.pdf",coll.c_str(), algo.c_str(), radius, jetType.c_str()),"RECREATE");

  // plot them on top of each other. 
  TCanvas * cDeltaPhi = new TCanvas("DeltaPhi","",800,600);
  hDeltaPhi_Data->SetMarkerStyle(20);
  hDeltaPhi_Data->SetMarkerColor(kBlack);
  hDeltaPhi_Data->SetAxisRange(0.0, 2.5, "X");
  hDeltaPhi_Data->SetYTitle("Event Fraction");
  hDeltaPhi_Data->SetXTitle("DeltaPhi");
  hDeltaPhi_Data->SetTitle(" ");
  hDeltaPhi_Data->DrawNormalized();
  
  hDeltaPhi_MC->SetMarkerStyle(25);
  hDeltaPhi_MC->SetMarkerColor(kRed);
  hDeltaPhi_MC->DrawNormalized("same");

  putCMSPrel();
  //putPPLumi();
  
  TLegend * lDeltaPhi = myLegend(0.4,0.6,0.7,0.9);
  lDeltaPhi->AddEntry("","p_{T}^{lead}   > 100 GeV/c","");
  lDeltaPhi->AddEntry("","p_{T}^{sublead}> 40  GeV/c","");
  lDeltaPhi->AddEntry(hDeltaPhi_Data,"Data (Jet 80 trigger)","pl");
  lDeltaPhi->AddEntry(hDeltaPhi_MC,"MC (pthat 80)","pl");
  lDeltaPhi->Draw();

  cDeltaPhi->SaveAs(Form("DeltaPhi_datavsMC_%s_ak%s%d%s.pdf",coll.c_str(), algo.c_str(), radius, jetType.c_str()),"RECREATE");

    // plot them on top of each other. 
  TCanvas * cpt2overpt1 = new TCanvas("pt2overpt1","",800,600);
  pt2overpt1_Data->SetMarkerStyle(20);
  pt2overpt1_Data->SetMarkerColor(kBlack);
  pt2overpt1_Data->SetAxisRange(0.0, 1, "X");
  pt2overpt1_Data->SetYTitle("Event Fraction");
//.........这里部分代码省略.........

void apply_unfolding(){

  TFile *unf_file = new TFile("unfolding.root","read");

  TFile *outfile = new TFile("Unfolding_Results.root","recreate");  

  for (int i=0; i<3; i++){
  for (std::vector<TString>::const_iterator diffvariable = diffvariables_list.begin(); diffvariable!=diffvariables_list.end(); diffvariable++){

    if (*diffvariable=="dR") continue;

    TString reg;
    if (i==0) reg="EBEB"; else if (i==1) reg="EBEE"; else if (i==2) reg="EEEE";

    RooUnfoldResponse *resp = NULL;
    unf_file->GetObject(Form("response_%s_%s",reg.Data(),diffvariable->Data()),resp);
    resp->Print();

    TFile *histo_file = new TFile(Form("plots/histo_xsec_%s_%s.root",diffvariable->Data(),reg.Data()),"read");
    TH1F *histo_reco = NULL;
    histo_file->GetObject("xsec_ngammagammayield",histo_reco);
    histo_reco->Print();

    RooUnfoldBayes *algo = new RooUnfoldBayes(resp,histo_reco,4);
    TH1D *_unfolded = (TH1D*)(algo->Hreco());
    _unfolded->Print();

    TH1F unfolded;
    _unfolded->Copy(unfolded);
    
    std::map<TString,TString> translation2;
    translation2.insert(std::pair<TString,TString>(TString("invmass"),TString("mgg")));
    translation2.insert(std::pair<TString,TString>(TString("diphotonpt"),TString("pt")));
    translation2.insert(std::pair<TString,TString>(TString("costhetastar"),TString("costt")));
    translation2.insert(std::pair<TString,TString>(TString("dphi"),TString("phi")));
    translation2.insert(std::pair<TString,TString>(TString("dR"),TString("dR")));
    unfolded.SetName(Form("Unfolding_Nevt_%s_%s",translation2[TString(diffvariable->Data())].Data(),reg.Data()));
    unfolded.SetTitle(Form("Unfolding_Nevt_%s_%s",translation2[TString(diffvariable->Data())].Data(),reg.Data()));
    
    unfolded.Print();

    TH1F unfolded_relsyserr;
    unfolded.Copy(unfolded_relsyserr);
    unfolded_relsyserr.SetName(Form("Unfolding_RelativeSysErr_%s_%s",translation2[TString(diffvariable->Data())].Data(),reg.Data()));
    unfolded_relsyserr.SetTitle(Form("Unfolding_RelativeSysErr_%s_%s",translation2[TString(diffvariable->Data())].Data(),reg.Data()));
    unfolded_relsyserr.Reset();
    unfolded_relsyserr.Print();

    outfile->cd();
    unfolded.Write();
    unfolded_relsyserr.Write();

  }
  }

  outfile->Close();





}