C++ TH1::Rebin方法代码示例

// -----------------------------------------------------------------------------
//
TH1* getHisto( std::string nameFile,
	       std::string nameHist,
	       std::string Dirname, 
	       int rebin ) {
  std::string name = nameFile;
  TFile* file =  new TFile(name.c_str());
  if (file) { std::cout << "Opened file: " << file->GetName() << std::endl; }
  else { 
    std::cout << "Could not find file: " << name << std::endl; 
    return 0; 
  }
  
  TDirectory* dir = (TDirectory*)file->Get(Dirname.c_str());
  if (dir) { std::cout << "Opened dir: " << dir->GetName() << std::endl; }
  else { 
    std::cout << "Could not find dir: " << Dirname << std::endl; 
    return 0; 
  }
  
  int low = 375;
  TH1* hist = 0;
  if ( false || nameHist.find("HtMultiplicity_HT375") == std::string::npos ) { 
    hist = (TH1*)dir->Get(nameHist.c_str());
  } else {
    
    for ( uint ii = low; ii <= 975; ii+=100 ) {
      std::stringstream tmp; tmp << "HtMultiplicity_HT" << ii << nameHist.substr(20);
      if ( !hist ) { 
	dir->cd();
	TH1D* temp = (TH1D*)dir->Get( "HtMultiplicity_HT375_aT0" );
	//TH1D* temp = (TH1D*)file->Get( tmp.str().c_str() );
	if (temp) { hist = (TH1D*)temp->Clone(); } 
	else { std::cout << "1 Unable to retrieve histo with name " << tmp.str() << std::endl; }
      } else { 
	dir->cd();
	TH1D* temp = (TH1D*)dir->Get( tmp.str().c_str() );
	if (temp) { hist->Add( (TH1D*)temp ); } 
	else { std::cout << "2 Unable to retrieve histo with name " << tmp.str() << std::endl; }
      }
    }

  }
  if (hist) { std::cout << "Opened histo: " << hist->GetName() << std::endl; }
  else { 
    std::cout << "Could not find histo: " << nameHist << std::endl; 
    return 0; 
  }

  hist->SetLineWidth(3);
  if ( rebin > 0 ) { hist->Rebin(rebin); }
  hist->GetXaxis()->SetTitleSize(0.055);
  hist->GetYaxis()->SetTitleSize(0.055);
  hist->GetXaxis()->SetLabelSize(0.05);
  hist->GetYaxis()->SetLabelSize(0.05);
  hist->SetStats(kFALSE);
  return hist;
}

TH1* loadHistogram(TFile* inputFile, const std::string& directory, const std::string& histogramName)
{
  std::string histogramName_full = Form("%s/%s", directory.data(), histogramName.data());
  TH1* histogram = dynamic_cast<TH1*>(inputFile->Get(histogramName_full.data()));
  if ( !histogram ) {
    std::cerr << "Failed to load histogram = " << histogramName_full << " from file = " << inputFile->GetName() << " !!" << std::endl;
    assert(0);
  }
  if ( !histogram->GetSumw2N() ) histogram->Sumw2();
  histogram->Rebin(4);
  return histogram;
}

TH1* getHistogram(TFile* inputFile, const TString& dqmDirectory, const TString& meName)
{  
  TString histogramName = dqmDirectory;
  if ( histogramName.Length() > 0 && !histogramName.EndsWith("/") ) histogramName.Append("/");
  histogramName.Append(meName);

  TH1* histogram = (TH1*)inputFile->Get(histogramName.Data());
  std::cout << "histogramName = " << histogramName.Data() << ": histogram = " << histogram;
  if ( histogram ) std::cout << ", integral = " << histogram->Integral();
  std::cout << std::endl; 

  if ( !histogram->GetSumw2N() ) histogram->Sumw2();

  if ( histogram->GetDimension() == 1 ) histogram->Rebin(5);

  return histogram;
}

// A function that get histogram and sets contents to 0 
// if entries are too small
TH1 * getHisto(TFile * file, const char * name, unsigned int rebin) {
  TObject * h = file->Get(name);
  if(h == 0)
    throw edm::Exception(edm::errors::Configuration) 
      << "Can't find object " << name << "\n";
  TH1 * histo = dynamic_cast<TH1*>(h);
  if(histo == 0)
    throw edm::Exception(edm::errors::Configuration) 
      << "Object " << name << " is of type " << h->ClassName() << ", not TH1\n";
  histo->Rebin(rebin);  
  for(int i = 1; i <= histo->GetNbinsX(); ++i) {
    if(histo->GetBinContent(i) < 0.1) {
      histo->SetBinContent(i, 0.0);
      histo->SetBinError(i, 0.0);
    }
  }
  return histo;
}

/** 
 * Draw the Poisson estimate of the occupancy in a given ring.
 * 
 * @param p            List 
 * @param d            Detector
 * @param r            Ring
 * 
 * @return The occupancy (in percent)
 *
 * @deprecated Use QATrender instead
 * @ingroup pwglf_forward_scripts_qa
 */
Double_t
DrawRingOccupancy(TList* p, UShort_t d, Char_t r)
{
  if (!p) return 0;

  TList* ring = static_cast<TList*>(p->FindObject(Form("FMD%d%c",d,r)));
  if (!ring) { 
    Error("DrawOccupancy", "List FMD%d%c not found in %s",d,r,p->GetName());
    return 0;
  }
  
  TH1* corr = static_cast<TH1*>(ring->FindObject("occupancy"));
  if (!corr) { 
    Error("DrawRingOccupancy", "Histogram occupancy not found in FMD%d%c",
	  d, r);
    return 0;
  }
  corr->Rebin(4);

  TPad* pad = static_cast<TPad*>(gPad);
  pad->SetGridy();
  pad->SetGridx();
  pad->SetLogy();
  pad->SetFillColor(0);
    pad->SetRightMargin(0.01);
#if 0
  if (d == 3) { 
    pad->SetPad(pad->GetXlowNDC(), pad->GetYlowNDC(), .99, 
		 pad->GetYlowNDC()+pad->GetHNDC());
    pad->SetRightMargin(0.15);
  }
#endif

  corr->Draw("hist");

  TLatex* ltx = new TLatex(.95, .95, Form("FMD%d%c", d, r));
  ltx->SetNDC();
  ltx->SetTextAlign(33);
  ltx->SetTextSize(.08);
  ltx->Draw();

  return corr->GetMean();
}

TH1* getMonitorElement(TFile* inputFile, const TString& dqmDirectory, const char* dqmSubDirectory, const TString& meName)
{
  TString meName_full = TString("DQMData").Append("/");
  if ( dqmDirectory != "") meName_full.Append(dqmDirectory).Append("/");
  meName_full.Append(dqmSubDirectory).Append("/").Append(meName);
  std::cout << "meName_full = " << meName_full << std::endl;
  
  TH1* me = (TH1*)inputFile->Get(meName_full);
  std::cout << "me = " << me <<  std::endl;
  
  //if ( !me->GetSumw2() ) me->Sumw2();
  me->Sumw2();

  me->Rebin(2);

  me->Scale(1./me->Integral());

  me->SetMaximum(1.);
  me->SetStats(false);

  return me;
}

TCanvas* overlay_Merged_RecoSmeared(const vector<string>& folders, 
								const hist_t& h
								)
{

	TLegend *leg  = new TLegend(0.6,0.65,0.9,0.9);
	leg->SetTextFont(42);
	vector<TH1*> hists;
	
	vector<string> jetcoll;
	jetcoll.push_back("reco");
	jetcoll.push_back("gen");
	jetcoll.push_back("smeared");

	stringstream title;
	const string njets("3-5");
	//const string eta("2.5");
	const string eta("5.0");
	title << njets
			<< " Jets, MHT from Jets with P_{T}>30 GeV, |#eta |<" 
			<< eta << ", L = 10 fb^{-1}" << ";" << h.title ;

	for (unsigned j=0; j< jetcoll.size(); ++j)
	{
		TH1* Hist = 0;
		for (unsigned i = 0; i < folders.size(); ++i)
		{
			stringstream histname;
			histname << folders.at(i) << "/" << jetcoll.at(j) << h.name;
			cout << __LINE__ << ": Looking for hist: " << histname.str().c_str() << endl;
			TH1* htemp = GetHist(histname.str());
			
			if (Hist == 0) Hist = htemp;
			else Hist->Add(htemp);
		}
		
		Hist->Rebin(h.rebin);
		Hist->SetTitle(title.str().c_str());
		Hist->SetMarkerStyle(20+j);
		Hist->SetLineWidth(2);
		Hist->SetStats(0);

		stringstream legname;
		if (j==0) 
		{
			legname << "Reco"; 
		} else if (j==1) 
		{
			legname << "Gen"; 
			Hist->SetLineColor(kBlue);
			Hist->SetMarkerColor(kBlue);
		} else if (j==2)
		{
			legname << "Smeared"; 
			Hist->SetLineColor(kRed);
			Hist->SetMarkerColor(kRed);
		}

		const double sum = Hist->Integral(); 
		legname << " (" << sum << ")";
		if (j!=1) leg->AddEntry(Hist, legname.str().c_str());

		hists.push_back(Hist);

	} //end jetcoll
	
	TH1* ratio = dynamic_cast<TH1*> (hists.at(2)->Clone("ratio"));
	ratio->GetYaxis()->SetTitle("Smeared/Reco");
	ratio->SetTitle("");
	ratio->Divide(hists.at(0));
	ratio->GetYaxis()->SetRangeUser(-0.01,2.01);
	//ratio->SetTickLength (+0.01,"Y");

   //TCanvas *c1 = new TCanvas("c1", "c1",15,60,550,600);
   TCanvas *c1 = new TCanvas("c1");
   c1->Range(0,0,1,1);
   c1->SetBorderSize(2);
   c1->SetFrameFillColor(0);
  
// ------------>Primitives in pad: c1_1
   TPad *c1_1 = new TPad("c1_1", "c1_1",0.01,0.30,0.99,0.99);
   c1_1->Draw();
   c1_1->cd();
   c1_1->SetBorderSize(2);
   c1_1->SetTickx(1);
   c1_1->SetTicky(1);
   c1_1->SetTopMargin(0.1);
   c1_1->SetBottomMargin(0.0);
   //c1_1->SetFrameFillColor(3);
	c1_1->SetLogy();
   
	hists.at(0)->GetYaxis()->CenterTitle(1);
	hists.at(0)->SetLabelFont(42,"XYZ");
	hists.at(0)->SetTitleFont(42,"XYZ");
	hists.at(0)->GetYaxis()->SetTitleOffset(0.8);
	hists.at(0)->SetLabelSize(0.05,"XYZ");
	hists.at(0)->SetTitleSize(0.06,"XYZ");
   hists.at(0)->Draw("P");
   hists.at(2)->Draw("same P");
	leg->Draw();
//.........这里部分代码省略.........

//______________________________________________________________________________
void CheckTime(const Char_t* loc)
{
    // Check time.

    Char_t t[256];

    // number of runs
    Int_t nRuns = gFiles->GetSize();

    // create arrays
    const Int_t nCh = 352;
    Double_t** pedPos = new Double_t*[nCh];
    Double_t* runNumbersD = new Double_t[nRuns];
    for (Int_t i = 0; i < nCh; i++) pedPos[i] = new Double_t[nRuns];

    // open the output files
    TFile* fROOTout = new TFile("/tmp/tagger_time.root", "RECREATE");

    // create directories
    for (Int_t i = 0; i < nCh; i++)
    {
        sprintf(t, "%03d", i);
        fROOTout->mkdir(t);
    }
    
    TF1* func = new TF1("func", "gaus(0)+pol1(3)", -1000 , 1000);

    // loop over runs
    for (Int_t i = 0; i < nRuns; i++)
    {   
        // get the file
        TFile* f = (TFile*) gFiles->At(i);

        // extract run number
        Int_t runNumber;
        sprintf(t, "%s/ARHistograms_CB_%%d.root", loc);
        sscanf(f->GetName(), t, &runNumber);
        runNumbersD[i] = (Double_t)runNumber;

        printf("Processing run %d (%d/%d)\n", runNumber, i+1, nRuns);

        fROOTout->cd();
        
        TH2* h2 = (TH2*) f->Get("CaLib_Tagger_Time");

        // loop over channels
        for (Int_t j = 0; j < nCh; j++)
        {
            // load histogram
            sprintf(t, "%03d", j);
            fROOTout->cd(t);
            sprintf(t, "%d_%d", i, j);
            TH1* h = h2->ProjectionX(t, j+1, j+1);
            h->Rebin(2);
            
            sprintf(t, "Run_%d", runNumber);
            TCanvas* c = new TCanvas(t, t);
            TLine* tline = 0;
            
            // check entries
            if (h->GetEntries())
            {
                // fit gaussian to peak
                Double_t maxPos = h->GetXaxis()->GetBinCenter(h->GetMaximumBin());
                func->SetParameters(1, maxPos, 0.5, 1, 0.1);
                func->SetRange(maxPos - 2, maxPos + 2);
                func->SetParLimits(0, 0, 1000);
                for (Int_t k = 0; k < 10; k++)
                    if (!h->Fit(func, "RBQ")) break;

                // save position in file and memory
                maxPos = func->GetParameter(1);
                pedPos[j][i] = maxPos;

                h->GetXaxis()->SetRangeUser(maxPos - 10, maxPos + 10);
                h->Draw();
                
                tline = new TLine(maxPos, 0, maxPos, 10000);
                tline->SetLineColor(kRed);
                tline->SetLineWidth(2);
                tline->Draw();
            }
            else
            {
                h->Draw();
            }

            c->Write(c->GetName(), TObject::kOverwrite);
    
            delete h;
            delete c;
            if (tline) delete tline;
        }
        
        delete h2;
    }
    
    // create pedestal evolution graphs
    fROOTout->cd();
//.........这里部分代码省略.........

void comparisonJetMCData(string plot,int rebin){
  string tmp;

  string dir="/gpfs/cms/data/2011/Observables/Approval/";
	
  if (isAngularAnalysis){
    mcfile=dir+"MC_zjets"+version;
    back_w=dir+"MC_wjets"+version;
    back_ttbar=dir+"MC_ttbar"+version;
    WW=dir+"MC_diW"+version;
    ZZ=dir+"MC_siZ"+version;
    WZ=dir+"MC_diWZ"+version;
    datafile=dir+"DATA"+version;
    mcfiletau=dir+"MC_zjetstau"+version;
  }
  // List of files

  TFile *dataf = TFile::Open(datafile.c_str()); //data file
  TFile *mcf = TFile::Open(mcfile.c_str()); //MC file
  TFile *mcftau = TFile::Open(mcfiletau.c_str()); //MC file
  TFile *ttbarf = TFile::Open(back_ttbar.c_str()); //MC background file
  TFile *wf = TFile::Open(back_w.c_str());


  TFile *qcd23emf = TFile::Open(qcd23em.c_str());
  TFile *qcd38emf = TFile::Open(qcd38em.c_str());
  TFile *qcd817emf = TFile::Open(qcd817em.c_str());

  TFile *qcd23bcf = TFile::Open(qcd23bc.c_str());
  TFile *qcd38bcf = TFile::Open(qcd38bc.c_str());
  TFile *qcd817bcf = TFile::Open(qcd817bc.c_str());

  TFile *WZf = TFile::Open(WZ.c_str());
  TFile *ZZf = TFile::Open(ZZ.c_str());
  TFile *WWf = TFile::Open(WW.c_str());


  // Canvas
  if (CanvPlot) delete CanvPlot;
  CanvPlot = new TCanvas("CanvPlot","CanvPlot",0,0,800,600);

  // Getting, defining ...
  dataf->cd("validationJEC");
  if (isMu && isAngularAnalysis) dataf->cd("validationJECmu");
	
  TObject * obj;
  gDirectory->GetObject(plot.c_str(),obj);

  TH1 *data;
  TH2F *data2; 
  TH1D *data3; 

  THStack *hs = new THStack("hs","Total MC");


  int flag=-1;
  if ((data = dynamic_cast<TH1F *>(obj)) ){
    flag=1;
    gROOT->Reset();
    gROOT->ForceStyle();
    gStyle->SetPadRightMargin(0.03);
    gPad->SetLogy(1);
    gPad->Modified();
    gPad->Update();
  }
  if ((data2 = dynamic_cast<TH2F *>(obj)) ){
    flag=2;
    gStyle->SetPalette(1);
    gStyle->SetPadRightMargin(0.15);
    gPad->Modified();
  }



  //===================
  // Dirty jobs :)
  if (flag==1){
    CanvPlot->cd();
    TPad *pad1 = new TPad("pad1","pad1",0.01,0.33,0.99,0.99);
    pad1->Draw();
    pad1->cd();
    pad1->SetTopMargin(0.1);
    pad1->SetBottomMargin(0.01);
    pad1->SetRightMargin(0.1);
    pad1->SetFillStyle(0);
    pad1->SetLogy(1);
    TString str=data->GetTitle();
    if (str.Contains("jet") && !str.Contains("zMass") && !str.Contains("Num") && !str.Contains("Eta") && !str.Contains("Phi") && !str.Contains("eld") && !str.Contains("meanPtZVsNjet")) {
      if (!isAngularAnalysis) rebin=1;
    }

    //======================
    // DATA
    Double_t dataint = data->Integral();
    data->SetLineColor(kBlack);
    data->Rebin(rebin);
    if(str.Contains("nJetVtx")) data->GetXaxis()->SetRangeUser(0,10);	
    if(str.Contains("zMass")) data->GetXaxis()->SetRangeUser(70,110);	
    data->SetMinimum(1.);
    data->Sumw2();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	exp_stats->SetY1NDC(ymin);
	exp_stats->SetY2NDC(ymax);
	exp_stats->Draw("same");

	//fit quality
	stringstream gaus_fit_res, exp_fit_res;
	const float gausFit_goodness =  gausFit->GetChisquare()/gausFit->GetNDF();
	gaus_fit_res << "#chi^{2}/ndof = " << gausFit->GetChisquare()
			 << "/"<< gausFit->GetNDF() << " = " << gausFit_goodness;
	const float expFit_goodness =  expFit->GetChisquare()/expFit->GetNDF();
	exp_fit_res << "#chi^{2}/ndof = " << expFit->GetChisquare()
			 << "/"<< expFit->GetNDF() << " = " << expFit_goodness;

	cout << gaus_fit_res.str() << endl;
	cout << exp_fit_res.str() << endl;

	TPaveText *pt1 = new TPaveText(0.5,0.8,0.7,0.9);
	pt1->AddText(gaus_fit_res.str().c_str());
	pt1->AddText(exp_fit_res.str().c_str());
//	pt1->Draw("same");

	//stringstream epsname;
	//epsname << "factnomht/HTbin" << HTbin << "_fitrange_" << fitrange_xmin << "to" << fitrange_xmax << ".eps";
	//epsname << "factnomht_" << HTbinlabel << ".eps";
	//gPad->Print(epsname.str().c_str());
	gPad->Print("ratios.eps");
	//return;

	 //make a predicion
	TH1 *sigHist      = GetHist(signalHistName, scaleTo);
	TH1 *controlHist  = GetHist(controlHistName, scaleTo);

	//temp hack to get # of bins the same for easy debugging
	sigHist->Rebin(2);

	cout << red << "sigHist/controlHist bins = " << sigHist->GetNbinsX() << "/" << controlHist->GetNbinsX() << clearatt << endl; 
//	new TCanvas();
//	sigHist->SetLineColor(kRed);
//	sigHist->Rebin(50);
//	controlHist->Rebin(50);
//	sigHist->Draw();
//	controlHist->Draw("same");


	//collect results
/*	Predictions_t pred_gaus       = GetPredictions(controlHist, gausFit2, sigHist); 
	Predictions_t pred_gaus_sigma = GetPredictions(controlHist, gausFit_sigma, sigHist); 
	Predictions_t pred_exp        = GetPredictions(controlHist, expFit2, sigHist); 
	Predictions_t pred_exp_sigma  = GetPredictions(controlHist, expFit_sigma, sigHist); 

//	PrintExclResults(pred_gaus);
//	PrintExclResults(pred_exp);
	PrintExclResults(pred_gaus, pred_exp, pred_gaus_sigma, pred_exp_sigma);
*/
	/*****************************************/
	// TEMP HACK TO GET RESULTS for all HT bins 
	// using these inclusive fits
	/*****************************************/
	vector<pair<float, float> > htBins_temp;
	pair<float, float> htbin1(500,750);	
	pair<float, float> htbin2(750,1000);	
	pair<float, float> htbin3(1000,1250);	
	pair<float, float> htbin4(1250,1500);	
	pair<float, float> htbin5(1500,8000);	

	htBins_temp.push_back(htbin1);

 void MejMassFits(TH1F* HistName) {

   // set up RooFit
   gSystem->Load("libRooFit");
   using namespace RooFit;
   
   gROOT->SetStyle("Plain");
   
   // get histogram and convert it to a RooFit understandble format
   TH1* hh = HistName->Clone("hh");
   hh->Rebin(1);
   RooRealVar Mej("Mej","Mej LQ",40.,1000,"GeV");
   RooDataHist RDH_data("RDH_data","dataset with Mej mass",Mej,hh);
   RooHistPdf Pdf_data("dataMej","dataMej" ,RooArgList(Mej),RDH_data); 

   //gen data accordingly to histo
   RooDataSet *data = Pdf_data.generate(Mej,5000) ;
   
   // make signal function
   //    RooRealVar mean("mean","mean",250,200,300);
   //    RooRealVar sigma("sigma","sigma",100,0.,200.);
   //    RooGaussian sig("sig","signal pdf",Mej,mean,sigma);
   
   RooRealVar mean("mean","mean",650,600.,700.);
   RooRealVar sigma("sigma","sigma",30,0.,100);
   RooRealVar alpha("alpha","alpha",0.1,0.,1.);
   RooRealVar n("n","n",0.1,0.,+100);
   RooCBShape sig("sig","signal pdf",Mej,mean,sigma,alpha,n);


   // make background function
   RooRealVar c0("c0","coeff #0",1.0,-5.,5.);
   RooRealVar c1("c1","coeff #1",0.1,-5.,5.);
   RooRealVar c2("c2","coeff #2",0.1,-5.,5.);
   RooRealVar c3("c3","coeff #3",0.1,-5.,5.); 
   RooChebychev bkg("bkg","bkgd pdf",Mej,RooArgList(c0,c1,c2,c3));

   // --- Build Argus background PDF ---
   //    RooRealVar argpar("argpar","argus shape parameter",-20.0,-100.,-1.) ;
   //    RooRealVar cutoff("cutoff","argus cutoff",1001,1000,2000) ;
   //    RooRealVar p("p","p",0.5,0,10);
   //    RooArgusBG bkg("argus","Argus PDF",Mej,cutoff,argpar,p);
   
   // we are interested in getting the number of signal/background events
   RooRealVar nsig("nsig","weel reco signal",500,0.,10000.);
   RooRealVar nbkg("nbkg","bad reco signal",500,0.,50000.);
   
   //define the model
   RooAddPdf model("model","signal+bkgd",RooArgList(sig,bkg),RooArgList(nsig,nbkg));
   
   //fit (notice use of Extended(true) - this is so that we can get nsig/nbkg
   RooFitResult *fitres = model.fitTo(*data,Extended(kTRUE),Minos(kFALSE),PrintLevel(0),Save(kTRUE));
   
   // make new canvas  
   TCanvas *TopMass = new TCanvas("TopMass","Top Mass");
   
   // make new frame and plot histogram and fit results on it 
   RooPlot *Massframe = Mej.frame();
   
   data.plotOn(Massframe,MarkerColor(kBlack));
   model.plotOn(Massframe,Components(sig),LineStyle(kDashed),LineColor(kRed));
   model.plotOn(Massframe,Components(bkg),LineStyle(kDashed),LineColor(kGreen));
   model.plotOn(Massframe,LineColor(kBlack));
   
   Massframe->Draw();
   
   // print the results of the fit   
   fitres->Print(); 
   
   // You have to get the chi-square of the fit from MassFrame - the 7 tells RooFit to take into account the 7 fit parameters when calculating the number of degrees of freedom
   cout<<" Fit chi square/dof = "<<Massframe->chiSquare(10)<<endl; 
   
 }

int ratio5() {

  // Constants and arrays

  Int_t multi = 2;

  const Int_t n_at = 3;
  Int_t at[n_at] = { 500, 510, 550 };
  //for ( int ii = 0; ii < n_at; ++ii ) { at[ii] = 500 + ii * 10; }

  TString eq = "Gt";

  const Int_t n = 4;
  float pt[n]     = { 50., 40., 30., 20. };
  Int_t colour[n] = { 1, 2, 3, 4 };

  const Int_t m = 2;
  Int_t style[m]  = { kOpenSquare, kFullSquare };
  
  const Int_t ngr = 1000;
  double x3[ngr];
  double r[ngr];

  int count = 0;

  // General style

  gStyle->SetOptStat(0);
  
//   // Canvas for RECO curves 
//   TCanvas* reco_canvas = new TCanvas("Reco");
//   reco_canvas->SetFillColor(0);
//   reco_canvas->SetLineColor(0); 
//   reco_canvas->SetLogy();
//   TLegend* reco_legend = new TLegend( 0.5, 0.7, 0.88, 0.88, NULL, "brNDC" );
//   reco_legend->SetFillColor(0);
//   reco_legend->SetLineColor(0); 
//   bool empty = true;
//   double reco_max = 1.e-15.;
//   double reco_min = 1.e15;
  
  // Loop through pt bins

  for ( Int_t i = 0; i < 1; ++i ) {

    std::stringstream pt_can;
    pt_can << "PtBin" << pt[i];
    
    // Canvas for Pt bin
    TCanvas* pt_canvas = new TCanvas(TString(pt_can.str()),"");
    pt_canvas->SetFillColor(0);
    pt_canvas->SetLineColor(0); 
    pt_canvas->SetLogy();
    TLegend* pt_legend = new TLegend( 0.82, 0.5, 0.98, 0.9, NULL, "brNDC" );

    pt_legend->SetFillColor(0);
    pt_legend->SetLineColor(0); 
    bool empty = true;
    double pt_max = 1.e-15.;
    double pt_min = 1.e15;
    std::vector<TH1*> pt_ratio;

    pt_canvas->SetRightMargin(0.2);

    // Open files
    std::stringstream ss;
    ss << "results/4/Reco" << pt[i] << "_QCDPythia6.root";
    TString name(ss.str());
    TFile* file = new TFile(name);
    if ( file->IsZombie() || !(file->IsOpen()) ) { continue; }
    file->cd();
    
    // Loop through AlphaT thresolds
    for ( Int_t iat = 0; iat < n_at; ++iat ) {
    
      // Loop through RECO and GEN
      for ( Int_t j = 0; j < m; ++j ) {
	
	// Define names of histos to open 
	std::stringstream pre;
	std::stringstream post;
	if ( j == 0 ) {
	  pre << "Ratio" << at[iat] << "/GenHt" << eq << "PreAlphaT" << at[iat] << "_" << multi;
	  post << "Ratio" << at[iat] << "/GenHt" << eq << "PostAlphaT" << at[iat] << "_" << multi;
	  std::cout << pre.str() << std::endl;
	  std::cout << post.str() << std::endl;
	} else if ( j == 1 ) {
	  pre << "Ratio" << at[iat] << "/Ht" << eq << "PreAlphaT" << at[iat] << "_" << multi;
	  post << "Ratio" << at[iat] << "/Ht" << eq << "PostAlphaT" << at[iat] << "_" << multi;
	  std::cout << pre.str() << std::endl;
	  std::cout << post.str() << std::endl;
	}
	
	// Create ratio histo
	TH1* denominator = his( (TH1*)file->Get(TString(pre.str())), 45, 200., 650. );
	TH1* numerator = his( (TH1*)file->Get(TString(post.str())), 45, 200., 650. );
	int rebin = 5;
	numerator->Rebin(rebin);
	denominator->Rebin(rebin);
	TH1* ratio = (TH1*)numerator->Clone();
//.........这里部分代码省略.........

void MakePi0Analysis(){

  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadRightMargin(0.01);
  gStyle->SetPadTopMargin(0.09);
  gStyle->SetPadBottomMargin(0.11);
  //gStyle->SetOptStat(0);
  //gStyle->SetOptTitle(1);
  //gStyle->SetPadTickX(1);
  //gStyle->SetPadTickY(1);
  TGaxis::SetMaxDigits(3);
  
  Double_t bins[] = {1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0, 3.2, 3.4, 3.6,3.8,4.0,4.5, 5.0, 5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,10.0,11,12}; //Nicolas LHC12d PHI
  const Int_t binnum = sizeof(bins)/sizeof(Double_t) - 1;
  
  TString cwd = gSystem->WorkingDirectory();

  TH1F* fHistPeakMean  = new TH1F("fHistPeakMean","",binnum,bins);
  TH1F* fHistPeakWidth = new TH1F("fHistPeakWidth","",binnum,bins);
  TH1F* fHistRawYield  = new TH1F("fHistRawYield","",binnum,bins);
  TH1F* fHistRawYieldPerEvent  = new TH1F("fHistRawYieldPerEvent","",binnum,bins);
  
  TCanvas *c1[binnum];
  TCanvas *c2[binnum];
  TCanvas *c3;

  TFile* ef     = TFile::Open("../AnalysisResults.root");
  TList* list   = (TList*)ef->Get("list_kINT7_Pi0");
  TList *listEv = (TList*)ef->Get("list_kINT7_Event");
  
  TH1F* fHistEvents = (TH1F*)listEv->FindObject("fHistAnalyzedEvents");
  Double_t nEvt = fHistEvents->GetEntries();
  
  TString tofName[] = {"","_TOFcut1","_TOFcut2"};
  TString fitName[] = {"CrystalBall","AsymmGauss"};
  TString modName[] = {"","_M1","_M3"};
  
  Int_t ntof = 3;
  
  Double_t signal_range_min=0.11;
  Double_t signal_range_max=0.15;
  Double_t bg_range_min=0.050;
  Double_t bg_range_max=0.250;
  
  for(Int_t itof=0; itof<ntof; ++itof){
    for(Int_t iMod=0; iMod<1; iMod++){
      for(Int_t iFit=0; iFit<2; iFit++){
	for(Int_t iPol=0; iPol<3; iPol++){
	  
	  Int_t pol = iPol;

	  TString fNameSame = "fHistMassTwoGammas"+tofName[itof]+modName[iMod];
	  TString fNameMix  = "fHistMixMassTwoGammas"+tofName[itof]+modName[iMod];
	  
	  THnSparse *fTHnSparseRecPi0;
	  THnSparse *fTHnSparseRecMixPi0;
	  THnSparse *fTHnSparseGeneratePi0;
	  
	  TH2* hPi0A08    = 0x0;
	  TH2* hMixPi0A08 = 0x0;
	  hPi0A08    = (TH2*)list->FindObject(fNameSame)->Clone();
	  hMixPi0A08 = (TH2*)list->FindObject(fNameMix)->Clone();

	  Int_t sbin = 0;
	  Int_t lbin = 0;
      
	  TH1F* fHistFinalRatio[binnum];
	  TH1F* fHistFinalSame[binnum];
	  TH1F* fHistFinalBG[binnum];
	  TH1F* fHistFinalSignal[binnum];
	  TH1F* fHistOnlyFitSignal[binnum];
	  
	  TF1* fFitFinalRatio[binnum];
	  TF1* fFitFinalSignal[binnum];
	  TF1* fFitOnlyFitSignal[binnum];

	  for(Int_t i=0; i<binnum; i++){    
	    
	    Double_t bin_width = (bins[i+1]-bins[i]);
	    Double_t bin_mean  = (bins[i+1]+bins[i])/2;
	    
	    Int_t rebin = 1;
	    
	    if(bin_mean>5.0){
	      rebin=10;
	    }
	    else{
	      rebin=5;
	    }
	    
	    Int_t sproj_bin = hPi0A08->GetYaxis()->FindBin(bins[i]);
	    Int_t lproj_bin = hPi0A08->GetYaxis()->FindBin(bins[i+1])-1;
	    
	    TH1* fHistBasicSame = (TH1*)hPi0A08   ->ProjectionX(Form("fHistBasicSame_No%d",i+1),sproj_bin,lproj_bin,""); 
	    TH1* fHistBasicMix  = (TH1*)hMixPi0A08->ProjectionX(Form("fHistBasicMix_No%d",i+1),sproj_bin,lproj_bin,""); 
	    fHistBasicSame->Rebin(rebin);
	    fHistBasicMix->Rebin(rebin);
	    fHistBasicSame->Sumw2();
	    fHistBasicMix->Sumw2();
	    fHistBasicSame->GetYaxis()->SetTitle(Form("dN/dM per %.0f MeV/c^{2}",fHistBasicSame->GetBinWidth(1)*1000));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
			/*	if ((hist2print.at(ihist).Name()).find("Jet"))
				{
				reco_hist_name << folder.str() << "reco" << hist2print.at(ihist).Name() << "_copy";
				reco_hist << folder.str() << "reco" << hist2print.at(ihist).Name();
				smear_hist_name << folder.str() << "smeared" << hist2print.at(ihist).Name() << "_copy";
				smear_hist << folder.str() << "smeared" << hist2print.at(ihist).Name();
				gen_hist_name << folder.str() << "gen" << hist2print.at(ihist).Name() << "_copy";
				gen_hist << folder.str() << "gen" << hist2print.at(ihist).Name();
				} else
				*/	{
					reco_hist_name << folder.str() << "reco_" << hist2print.at(ihist).Name() << "_copy";
					reco_hist << folder.str() << "reco_" << hist2print.at(ihist).Name();
					smear_hist_name << folder.str() << "smeared_" << hist2print.at(ihist).Name() << "_copy";
					smear_hist << folder.str() << "smeared_" << hist2print.at(ihist).Name();
					gen_hist_name << folder.str() << "gen_" << hist2print.at(ihist).Name() << "_copy";
					gen_hist << folder.str() << "gen_" << hist2print.at(ihist).Name();
				}

				TH1* hreco = (TH1*) (outRootFile->Get(reco_hist.str().c_str()));
				if (hreco == NULL) { cout << "hreco = " << reco_hist.str() << " was not found!" << endl; assert(false); } 
				hreco->SetDirectory(0);
				TH1* hsmear = (TH1*) (outRootFile->Get(smear_hist.str().c_str()));
				if (hsmear == NULL) { cout << "hsmear = " << smear_hist.str() << " was not found!" << endl; assert(false); } 
				hsmear->SetDirectory(0);
				TH1* hgen = (TH1*) (outRootFile->Get(gen_hist.str().c_str()));
				//->Clone(gen_hist_name.str().c_str()));
				if (hgen == NULL) { cout << "hgen = " << gen_hist.str() << " was not found!" << endl; assert(false); } 
				hgen->SetDirectory(0);

				hreco->Sumw2();
				hsmear->Sumw2();
				hgen->Sumw2();

				const int rebin = hist2print.at(ihist).Rebin();
				const string title = hist2print.at(ihist).Title();
				const double xmin = hist2print.at(ihist).Xmin();
				const double xmax = hist2print.at(ihist).Xmax();

				if (rebin>1)
				{
					hreco->Rebin(rebin);
					hsmear->Rebin(rebin);
					hgen->Rebin(rebin);
				}
				if (title.length()>0)
				{
					hreco->SetTitle(title.c_str());
					hsmear->SetTitle(title.c_str());
					hgen->SetTitle(title.c_str());
				}
				if (xmin != LargeNegNum || xmax != LargeNegNum)
				{
					hreco->GetXaxis()->SetRangeUser(xmin,xmax);
					hsmear->GetXaxis()->SetRangeUser(xmin,xmax);
					hgen->GetXaxis()->SetRangeUser(xmin,xmax);
				}

				const double reco_max_y  = hreco->GetBinContent(hreco->GetMaximumBin());
				const double smear_max_y = hsmear->GetBinContent(hsmear->GetMaximumBin());
				const double y_max = max(reco_max_y, smear_max_y);
				double y_min = 9999.0;
				for (unsigned bin=1; bin<hreco->GetNbinsX(); ++bin)
				{
					const double v1 = hreco->GetBinContent(bin);
					const double v2 = hsmear->GetBinContent(bin);
					const double minv = min(v1,v2);

void makeSplitQCDhist_PythiaBinned(vector<string> folders, const string histname, 
											const string htrange, const string htbinlabel,
											const hist_t histinfo)
{
	const float scaleTo = fDATA_LUMI; // pb

	TLegend *leg  = new TLegend(0.6,0.65,0.9,0.9);
	leg->SetTextFont(42);
	vector<TH1*> hists;
	new TCanvas();
	gPad->SetLogy();
	gPad->SetTickx();
	gPad->SetTicky();

		stringstream title;
		title << htrange << " [" << histinfo.name << "];" << histinfo.title;
		cout << title.str() << endl;
		TH1* Hist = GetHist(histname);

		Hist->SetTitle(title.str().c_str());
		Hist->SetMarkerStyle(20);
		Hist->SetStats(0);
		Hist->Rebin(histinfo.rebin);
		Hist->Draw("P");


/*
	for (unsigned i = 0; i < folders.size(); ++i)
	{	
		string njet("");
		if (i==0) njet += "[2-3]";
		else if (i==1) njet += "[4-5]";
		else if (i==2) njet += "[6-7]";
		else if (i==3) njet += "#geq 8";
*/		/*if (i==0) njet += "3";
		else if (i==1) njet += "4";
		else if (i==2) njet += "5";
		else if (i==3) njet += "6";
		else if (i==4) njet += "7";
		else if (i==5) njet += "8";*/

/*		stringstream title, histName;
		title << htrange << ";" << histinfo.title;
		cout << title.str() << endl;
		histName << folders.at(i) << "/" << histname;;
		hists.push_back(GetHist(histName.str()));
		hists.at(i)->SetTitle(title.str().c_str());
		hists.at(i)->SetMarkerStyle(20+i);
		hists.at(i)->SetMarkerColor(1+i*2);
		hists.at(i)->SetStats(0);
		hists.at(i)->Rebin(histinfo.rebin);
		if (histinfo.normalizeByBinWidth) NormByBinWidth(hists.at(i));

		stringstream legname;
		legname << "Njets " << njet;
		leg->AddEntry(hists.at(i), legname.str().c_str());
		if (i==0) hists.at(i)->Draw("P");
		else hists.at(i)->Draw("same P");
	}
	leg->Draw();
*/
}