C++ TH1D::GetName方法代码示例

void mc2ibd_NE(TChain *tMC, TFile *fBgnd, TCanvas *cv)
{
	char str[1024];
	TLatex *txt = new TLatex();

	TH1D *hExp = (TH1D *) fBgnd->Get("hNEA-diff");
	if (!hExp) {
		printf("Histogram hNEA-diff not found in %s\n", fBgnd->GetName());
		return;
	}
	hExp->SetTitle("Delayed event energy;MeV;Events/200 keV");
	hExp->SetLineColor(kBlack);
	hExp->SetLineWidth(3);

	gROOT->cd();
	TH1D *hMC = new TH1D("hNEMC", "Delayed event energy (MC);MeV;Events/200 keV", 45, 3, 12);
	hMC->SetLineColor(kBlue);
	tMC->Project(hMC->GetName(), "NeutronEnergy", cX && cY && cZ && cR && c20 && cGamma && cGammaMax && cPe);
	hMC->Sumw2();
	hMC->Scale(hExp->Integral(15, 45) / hMC->Integral(15,45));
	
	cv->Clear();
	hMC->Draw("hist");
	hExp->DrawCopy("same");
	TLegend *lg = new TLegend(0.65, 0.8, 0.89, 0.89);
	lg->AddEntry(hExp, "IBD", "LE");
	lg->AddEntry(hMC, "MC", "L");
	lg->Draw();
}

void PrintHistogram(TH1D h) {
    gStyle->SetOptStat("eimr");
    TCanvas c("c","c",800,500);
    h.Draw();
    system("mkdir -p parms_timewalk");
    c.Print("parms_timewalk/"+TString(h.GetName())+".gif");
}

void mc2ibd_R2(TChain *tMC, TFile *fBgnd, TCanvas *cv)
{
	char str[1024];
	TLatex *txt = new TLatex();

	TH1D *hExp = (TH1D *) fBgnd->Get("hR2A-diff");
	if (!hExp) {
		printf("Histogram hR2A-diff not found in %s\n", fBgnd->GetName());
		return;
	}
	hExp->SetTitle("Distance between positron and neutron, 3D case;cm;Events/4cm");
	hExp->SetLineColor(kBlack);
	hExp->SetLineWidth(3);

	gROOT->cd();
	TH1D *hMC = new TH1D("hR2MC", "Distance between positron and neutron, 3D case (MC);cm;Events/4cm", 40, 0, 160);
	hMC->SetLineColor(kBlue);
	tMC->Project(hMC->GetName(), "Distance", cX && cY && cZ && cRXY && c20 && cGamma && cGammaMax && cPe && cN);
	hMC->Sumw2();
	hMC->Scale(hExp->Integral() / hMC->Integral());
	
	cv->Clear();
	hExp->DrawCopy();
	hMC->Draw("hist,same");
	TLegend *lg = new TLegend(0.65, 0.8, 0.89, 0.89);
	lg->AddEntry(hExp, "IBD", "LE");
	lg->AddEntry(hMC, "MC", "L");
	lg->Draw();
}

void Draw_KL_Test(){

  TChain* ch = new TChain("Tree");
  TChain* ch1 = new TChain("Tree");
  
  TH1D* his  = new TH1D("Klong6g","Klong6g",20,450,550);
  TH1D* his1 = new TH1D("Klong4g","Klong4g",20,450,550);
  TH1D* his2 = new TH1D("Klong4gAll","Klong4gAll",60,250,550);

  for( int i = 0; i< 68; i++){
    ch->Add(Form("klongRootFile/kl%d.root" ,4162+i));
    ch1->Add(Form("klongRootFile/ks%d.root",4162+i));
  }
  ch->Project(his->GetName()  ,"KlongMass[0]","CutCondition==0");
  ch1->Project(his1->GetName(),"KlongMass[0]","CutCondition==0");
  ch1->Project(his2->GetName(),"KlongMass[0]","CutCondition==0");
  
  TF1* func = new TF1("func","gaus(0)+expo(3)",0,550);
  func->SetParameter(1,498);
  func->SetParameter(2,5);
  TF1* func2 = new TF1("func2","gaus(0)",0,550);
  func2->SetParameter(1,498);
  func2->SetParameter(2,5);



  TCanvas* can = new TCanvas("can","",1200,600);
  can->Divide(2,1);
  can->cd(1);
  his2->Fit(func->GetName(),"","",450,550);
  his2->Draw();
  TF1* func1 = new TF1("Test","gaus",450,550);
  func1->SetParameter(0,func->GetParameter(0));
  func1->SetParameter(1,func->GetParameter(1));
  func1->SetParameter(2,func->GetParameter(2));

  can->cd(2);
  his1->SetLineColor(2);  
  his->Draw();
  his->Fit(func2->GetName(),"","",450,550);
  func->Draw("same");
  his1->Draw("same");

  std::cout<< func2->GetParameter(0) << " " 
	   << func->GetParameter(0)  << " " 
	   << func->GetParameter(0)/func2->GetParameter(0)<< std::endl;
  std::cout<< func2->Integral(450,550) << " " 
	   << func1->Integral(450,550)  << " " 
	   << func1->Integral(450,550)/func2->Integral(450,550)
	   << std::endl;
  //ch->Draw("KlongPt[0]:KlongMass[0]>>(400,200,600,50,0,20)","(CutCondition&(1|2|4|8))==0","colz");
  gPad->SetLogz();
  TText* text = new TText(0.5,0.5,"");
  TText* text1 = new TText(0.5,0.5,"");
  text->DrawTextNDC(0.5,0.5,Form("Integral:%2.3lf",func1->Integral(450,550)));
  text1->DrawTextNDC(0.5,0.6,Form("Integral:%2.3lf",func2->Integral(450,550)));
  
}

scalePair energyCorrectionDiff(TCut centCut1, float lowPt, float highPt, float lowEta, float highEta,TCut addCut) { 
   
   TString fname1 = "forest/barrelHiForestPhoton_MCphoton50_51k.root";
   if ( lowPt > 90 ) 
      fname1 = "forest/barrelHiForestPhoton_MCphoton80_28k.root";
   TFile *f1  =new TFile(fname1.Data());
   TTree *photon1 = (TTree*)f1->Get("yongsunPhotonTree");
   photon1->AddFriend("yEvt=yongsunHiEvt"    ,fname1.Data());
   photon1->AddFriend("ySkim=yongsunSkimTree"   ,fname1.Data());
   photon1->AddFriend("yHlt=yongsunHltTree"     ,fname1.Data());
   
   TCut collisionCut = "ySkim.pcollisionEventSelection==1";
   
   TCut hoeCut = "hadronicOverEm<0.2";
   TCut isoCut = "cc4 + cr4 + ct4j20 < 5 && sigmaIetaIeta<0.011";
   TCut ptCut  = Form("genMatchedPt>%.f && genMatchedPt <%.f",lowPt, highPt);
   TCut etaCut = Form("abs(eta)>%f && abs(eta)<%f",lowEta,highEta);
      
   TCut finalCut1 = genMatchCut1 && collisionCut && centCut1 && hoeCut && isoCut && ptCut && etaCut && addCut ; 
   TString variable1 = "pt/genMatchedPt";
   
   TH1D* hScale = new TH1D("hScale","",100,.5,1.5);
   TH1D* hdpt = new TH1D("hdpt","",100,-20,20);
   photon1->Draw(Form("%s>>%s",variable1.Data(),hScale->GetName()), finalCut1);
   photon1->Draw(Form("pt-genMatchedPt>>%s",hdpt->GetName()),       finalCut1);
   cout << "cut = " << finalCut1.GetTitle() <<endl;
   hScale->Draw();
   hdpt->Draw();
   TF1* ff =  cleverGaus(hScale);
   
   scalePair ret;
   double *ps = ff->GetParameters();
   ret.val = ps[1];
   ret.err = ff->GetParError(1);
   // resErr = ff->GetParError(2);
   ret.res  = ps[2];
   ret.resErr = ff->GetParError(2);  
   ret.absVal = hdpt->GetMean();
   cout <<"scale = " << ret.val << " +-" << ret.err << endl;
   return ret;
}

void makeplot(const char *name, TTree *tree, TCut weight, const char *drawstring, const char *xlabel, int nbins, double xlow, double xhigh) {
  
  //this is for NLO with FXFX merging
//    TCut mult0 = "LHEEvent.npNLO()==0";
//    TCut mult1 = "LHEEvent.npNLO()==1";
//    TCut mult2 = "LHEEvent.npNLO()==2";
//    TCut mult3 = "LHEEvent.npNLO()==3";
  
  //this is for LO with MLM
   TCut mult0 = "GenEvent.nMEPartons()==0";
   TCut mult1 = "GenEvent.nMEPartons()==1";
   TCut mult2 = "GenEvent.nMEPartons()==2";
   TCut mult3 = "GenEvent.nMEPartons()==3";

  //this is for LO with MLM (plotting partons after excluding non-matched partons in wbb/vbf type processes)
  //   TCut mult0 = "GenEvent.nMEPartonsFiltered()==0";
  //   TCut mult1 = "GenEvent.nMEPartonsFiltered()==1";
  //   TCut mult2 = "GenEvent.nMEPartonsFiltered()==2";
  
  TH1D *hall = new TH1D(TString::Format("hall_%s",name),"",nbins,xlow,xhigh);
  TH1D *hmult0 = new TH1D(TString::Format("hmult0_%s",name),"",nbins,xlow,xhigh);
  TH1D *hmult1 = new TH1D(TString::Format("hmult1_%s",name),"",nbins,xlow,xhigh);
  TH1D *hmult2 = new TH1D(TString::Format("hmult2_%s",name),"",nbins,xlow,xhigh);
  TH1D *hmult3 = new TH1D(TString::Format("hmult3_%s",name),"",nbins,xlow,xhigh);
  
  hmult0->SetLineColor(kBlue);
  hmult1->SetLineColor(kRed);
  hmult2->SetLineColor(kMagenta);
  hmult3->SetLineColor(kGreen+1);
  
  tree->Draw(TString::Format("%s>>%s",drawstring,hall->GetName()),weight,"goff");
  tree->Draw(TString::Format("%s>>%s",drawstring,hmult0->GetName()),weight*mult0,"goff");
  tree->Draw(TString::Format("%s>>%s",drawstring,hmult1->GetName()),weight*mult1,"goff");
  tree->Draw(TString::Format("%s>>%s",drawstring,hmult2->GetName()),weight*mult2,"goff");
  tree->Draw(TString::Format("%s>>%s",drawstring,hmult3->GetName()),weight*mult3,"goff");
  
  hall->GetXaxis()->SetTitle(xlabel);
  
  TCanvas* c = new TCanvas(name,name);
  c->cd();
  hall->SetLineWidth(2);
  hall->Draw("EHIST");
  hmult0->SetLineWidth(2);
  hmult0->Draw("EHISTSAME");
  hmult1->SetLineWidth(2);
  hmult1->Draw("EHISTSAME");
  hmult2->SetLineWidth(2);
  hmult2->Draw("EHISTSAME");
  hmult3->SetLineWidth(2);
  hmult3->Draw("EHISTSAME");
  c->SetLogy();
  c->SaveAs(TString::Format("%s.pdf",name));
}

void photonSpectra()
{
  TH1::SetDefaultSumw2();

  const Double_t PTBINS[] = {40, 50, 60, 80, 120, 300};
  const Int_t nPTBINS = 5;

  const TString DATA_FILE = "gammaJets_pA_Data.root";
  const TString MC_FILE = "gammaJets_pA_MC_allQCDPhoton.root";
  TFile *dataFile = TFile::Open(DATA_FILE);
  TTree *dataTree = (TTree*)dataFile->Get("photonTree");
  TFile *mcFile = TFile::Open(MC_FILE);
  TTree *mcTree = (TTree*)mcFile->Get("photonTree");

  TH1D *photonSpectra = new TH1D("photonSpectra","photonSpectra", nPTBINS, PTBINS);
  TH1D *correctedPhotonSpectra = (TH1D*)photonSpectra->Clone("correctedPhotonSpectra");

  const TCut photonEtaCut = "abs(photonTree.eta) < 1.44";
  const TCut isoCut = "ecalRecHitSumEtConeDR04 < 4.2  &&  hcalTowerSumEtConeDR04 < 2.2  &&  trkSumPtHollowConeDR04 < 2 && hadronicOverEm<0.1";
  const TCut genCut = "genCalIsoDR04<5 && abs(genMomId)<=22";
  const TCut sbCut = "(cc4+cr4+ct4PtCut20>10) && (cc4+cr4+ct4PtCut20<20) && hadronicOverEm<0.1";
  const TCut candidateCut = "sigmaIetaIeta<0.01";
  const TCut decayCut = "(sigmaIetaIeta>0.011) && (sigmaIetaIeta<0.017)";

  TCut photonPtCut = "photonTree.corrPt>40 && photonTree.corrPt<300";

  dataTree->Project(photonSpectra->GetName(),"corrPt",isoCut && photonEtaCut && photonPtCut && candidateCut);

  for(int i = 0; i<nPTBINS; ++i)
  {
    photonPtCut = Form("photonTree.corrPt>%f && photonTree.corrPt<%f", PTBINS[i], PTBINS[i+1]  );

    TCut dataCandidateCut = isoCut && photonEtaCut && photonPtCut;
    TCut sidebandCut =  sbCut && photonEtaCut && photonPtCut;
    TCut mcSignalCut = dataCandidateCut && genCut;

    fitResult fitr = getPurity(dataTree, mcTree,
  			       dataCandidateCut, sidebandCut,
  			       mcSignalCut);
    Double_t purity = fitr.purity;
    correctedPhotonSpectra->SetBinContent(i+1,photonSpectra->GetBinContent(i+1)*purity);
  }

  // TCanvas *c1 = new TCanvas();
  // photonSpectra->Draw();
  // correctedPhotonSpectra->Draw("same");

  TFile *outFile = new TFile("photonSpectra.root", "RECREATE");
  outFile->cd();
  photonSpectra->Write();
  correctedPhotonSpectra->Write();
  outFile->Close();
}

void MakeHQTHiggsPtDistribution( string inputFilename ,  const string Label) {

  //--------------------------------------------------------------------------------------------------------------
  // Settings 
  //==============================================================================================================
  
  Double_t lumi;              // luminosity (pb^-1)
  string label = "";
  if (Label != "") label = "_" + Label;
 


  //--------------------------------------------------------------------------------------------------------------
  // Main analysis code 
  //==============================================================================================================  
  
  //
  // Access samples and fill histograms
  //  
  ifstream infile(inputFilename.c_str());
  vector<double> lowEdges;
  vector<double> binContent;
  double x_min;
  double x_max;
  double y;
  Int_t binNumber = 1;
  while(infile >> x_min) {
    infile >> x_max >> y;
    lowEdges.push_back(x_min);
    binContent.push_back(y);
  }

  Double_t *xLowEdges = new Double_t[lowEdges.size()];
  for (UInt_t i=0; i < lowEdges.size();i++) {
    xLowEdges[i] = lowEdges[i];
  }


  //--------------------------------------------------------------------------------------------------------------
  // Histograms
  //==============================================================================================================  
  TH1D *bosonSystemPt = new TH1D((string("HiggsPt")+ label).c_str(), ";  Higgs p_{T} [GeV/c]; Number of Events", lowEdges.size()-1, xLowEdges);

  for (UInt_t i=1;i<lowEdges.size();i++) {
    bosonSystemPt->SetBinContent(i,binContent[i]);
  }

  TFile *file = new TFile("HwwAcceptanceSystematics.root", "UPDATE");
  file->WriteTObject(bosonSystemPt, bosonSystemPt->GetName(), "WriteDelete");
  file->Close();

}

Spectrum roothist_to_spectrum(const TH1D & hist, const TH2D * cov, bool use_hist_uncertainties){
    const int n = hist.GetNbinsX();
    if(cov){
        if(cov->GetDimension() != 2 || cov->GetNbinsX() != n || cov->GetNbinsY() != n){
            throw runtime_error(string("covariance histogram '") + cov->GetName() + "' has wrong dimension");
        }
    }
    Spectrum s(n);
    for(int i=0; i<n; ++i){
        double c = hist.GetBinContent(i+1);
        if(!isfinite(c)){
            throw runtime_error(string("non-finite entry in histogram '") + hist.GetName() + "'");
        }
        s[i] = c;
        if(use_hist_uncertainties){
            double e = hist.GetBinError(i+1);
            if(!isfinite(e)){
                throw runtime_error(string("non-finite error in histogram '") + hist.GetName() + "'");
            }
            s.cov()(i,i) += e*e;
        }
    }
    if(cov){
        for(int i=0; i<n; ++i){
            for(int j=0; j<n; ++j){
                double c_ij = cov->GetBinContent(i+1, j+1);
                double c_ji = cov->GetBinContent(j+1, i+1);
                if(!isfinite(c_ij)){
                    throw runtime_error(string("covariance histogram '") + cov->GetName() + "' does have non-finite entry");
                }
                if(fabs(c_ij - c_ji) > 1e-8 * max(fabs(c_ij), 1.0)){
                    cerr << "covariance histogram '" << cov->GetName() <<  "' is not symmetric: C(" << i << ","<< j << ") = " << c_ij << "; transposed: " << c_ji << "; diff = " << (c_ij - c_ji) << endl;
                }
                s.cov()(i,j) += c_ij;
            }
        }
    }
    return s;
}

void CompareCFLADiffCuts()
{
  TFile newFile("~/Analysis/lambda/AliAnalysisLambda/Results/2016-01/15-NoOppChargeCut/All/CFs.root");
  TDirectory *newDir = newFile.GetDirectory("Merged");
  vector<TH1D*> newCFs;
  newCFs.push_back((TH1D*)newDir->Get("CFLamALam010"));
  newCFs.push_back((TH1D*)newDir->Get("CFLamALam1030"));
  newCFs.push_back((TH1D*)newDir->Get("CFLamALam3050"));

  TFile oldFile("~/Analysis/lambda/AliAnalysisLambda/Results/2016-01/08-NewAvgSepCuts/All/CFs.root");
  TDirectory *oldDir = oldFile.GetDirectory("Merged");
  vector<TH1D*> oldCFs;
  oldCFs.push_back((TH1D*)oldDir->Get("CFLamALam010"));
  oldCFs.push_back((TH1D*)oldDir->Get("CFLamALam1030"));
  oldCFs.push_back((TH1D*)oldDir->Get("CFLamALam3050"));

  TFile *compareFile = new TFile("Compare.root","update");
  TDirectory *dir = compareFile->GetDirectory("Delta");
  if(!dir) dir = compareFile->mkdir("Delta");
  for(UInt_t i = 0; i < newCFs.size(); i++) {
    // TH1D *ratio = (TH1D*)newCFs[i]->Clone();
    // TString name = ratio->GetName();
    // ratio->SetName(name + "Ratio");
    // ratio->SetTitle(name + "Ratio");
    // ratio->Divide(oldCFs[i]);

    // TH1D *barlowRatio = ComputeRogerBarlowRatio(newCFs[i], oldCFs[i]);
    // barlowRatio->SetName(name + "BarlowRatio");
    // barlowRatio->SetTitle(name + "BarlowRatio");

    TString name = newCFs[i]->GetName();
    TH1D *barlowDifference = ComputeRogerBarlowDifference(newCFs[i], oldCFs[i]);
    barlowDifference->SetName(name + "BarlowDifference");
    barlowDifference->SetTitle(name + "BarlowDifference");
    
    dir->cd();
    // ratio->Write();
    // barlowRatio->Write();
    barlowDifference->Write(barlowDifference->GetName(), TObject::kOverwrite);
    Chi2TestWithZero(barlowDifference);
    FitWithConstant(barlowDifference, compareFile);
    // LookAtMean(barlowDifference);
    RebinHist(barlowDifference, compareFile);
    ManuallyRebin(barlowDifference, compareFile);
  }
    
  compareFile->Close();
}

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

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

void ManuallyRebin(TH1D* h1, TFile *out)
{
  cout<<"Taking weighted average"<<endl;
  Int_t nCombine = 40; // Merge 100 bins into 1

  Int_t totalBins = h1->GetNbinsX();
  Int_t finalBins = totalBins/nCombine;
  cout<<"Final N bins:\t"<<finalBins<<endl;
  TString newName = h1->GetName();
  newName += "WgtAvg";
  TH1D *newH = new TH1D(newName, newName, finalBins,
			h1->GetXaxis()->GetXmin(),
			h1->GetXaxis()->GetXmax());
  newH->SetLineColor(kGreen);

  for(Int_t iNewBin = 0; iNewBin < finalBins; iNewBin++)
  {
    Double_t weightedSum = 0.;
    Double_t weights = 0.;
    
    for(Int_t iOldBin = 1; iOldBin < nCombine+1; iOldBin++)
    {
      Int_t binNum = iNewBin * nCombine + iOldBin;
      Double_t errorSq = pow(h1->GetBinError(binNum), 2.);
      Double_t content = h1->GetBinContent(binNum);

      weights += 1./errorSq;
      weightedSum += (1./errorSq) * content;
    }

    Double_t avg = weightedSum/weights;
    Double_t newErr = 1./sqrt(weights);

    newH->SetBinContent(iNewBin + 1, avg);
    newH->SetBinError(iNewBin + 1, newErr);
  }

  TDirectory *dir = out->GetDirectory("WeightedAvg");
  if(!dir) {
    dir = out->mkdir("WeightedAvg");
  }
  dir->cd();
  newH->SetDirectory(0);
  newH->Write(newH->GetName(), TObject::kOverwrite);
  
  // dir->Close();
  
}

void ecalSpectra()
{
  TH1::SetDefaultSumw2();

  //const Double_t PTBINS[] = {0,2.5,5,7.5,10,12.5,15,17.5, 20,25, 30,35, 40, 50, 60, 80, 120, 300};
  //const Int_t nPTBINS = 5;
  //const Int_t nPTBINS = sizeof(PTBINS)/sizeof(Double_t) -1;

  TChain *ecalTree = new TChain("ecalTree");
  ecalTree->Add("pA_promptRECO/*.root");
  //ecalTree->Add("pA_promptRECO/89.root");
  //TH1D *ecalSpectra = new TH1D("ecalSpectra","ecalSpectra", nPTBINS, PTBINS);
  TH1D *ecalSpectra = new TH1D("ecalSpectra","ecalSpectra",300,0,300);

  // const TCut ecalEtaCut = "abs(ecalTree.eta) < 1.44";
  // const TCut isoCut = "ecalRecHitSumEtConeDR04 < 4.2  &&  hcalTowerSumEtConeDR04 < 2.2  &&  trkSumPtHollowConeDR04 < 2 && hadronicOverEm<0.1";
  // const TCut genCut = "genCalIsoDR04<5 && abs(genMomId)<=22";
  // const TCut sbCut = "(cc4+cr4+ct4PtCut20>10) && (cc4+cr4+ct4PtCut20<20) && hadronicOverEm<0.1";
  // const TCut candidateCut = "sigmaIetaIeta<0.01";
  // const TCut decayCut = "(sigmaIetaIeta>0.011) && (sigmaIetaIeta<0.017)";

  //TCut ecalPtCut = "ecalTree.corrPt>40 && ecalTree.corrPt<300";

  //const TCut triggerCut = "HLT_PAPhoton10_NoCaloIdVL_v1 && HLT_PAPhoton10_NoCaloIdVL_v1_Prescl==1 && L1_SingleEG5_BptxAND_Prescl==1";
  const TCut triggerCut = "HLT_PAPhoton10_NoCaloIdVL_v1";
  const TCut showerCut = "sigmaIetaIeta > 0.011";
  const TCut etaCut = "abs(eta) < 1.44";
  //const TCut trackerIso = "trkSumPtHollowConeDR04 > 2";
  const TCut ecalIso = "ecalRecHitSumEtConeDR04 > 4.2";
  //const TCut hcalIso
  Int_t count = ecalTree->Project(ecalSpectra->GetName(),"pt", triggerCut && showerCut && etaCut && ecalIso);
  printf("Num events: %i\n", count);

  TCanvas *c1 = new TCanvas();
  ecalSpectra->Draw();
  // correctedEcalSpectra->Draw("same");
  c1->SetLogy();

  // TFile *outFile = new TFile("ecalSpectra.root", "RECREATE");
  // outFile->cd();
  // ecalSpectra->Write();
  // correctedEcalSpectra->Write();
  // outFile->Close();
}

void mc2ibd_PPX(char X, TChain *tMC, TFile *fBgnd, TCanvas *cv)
{
	char str[1024];
	TH1D *hExp;
	TH1D *hMC;
	TCut cut;

	sprintf(str, "hP%cA-diff", X);
	hExp = (TH1D *) fBgnd->Get(str);
	if (!hExp) {
		printf("Histogram hP%cA-diff not found in %s\n", X, fBgnd->GetName());
		return;
	}
	gROOT->cd();
	sprintf(str, "hMCPP%c", X);
	hMC = new TH1D(str, "MC XYZ", hExp->GetNbinsX(), 0, 100);
	sprintf(str, "PositronX[%d]+%4.1f", X - 'X', (X=='Z') ? 0.5 : 2.0);
	switch (X) {
	case 'X':
		cut = cY && "PositronX[0]>=0" && cZ && cR && c20 && cGamma && cGammaMax && cPe && cN;
		break;
	case 'Y':
		cut = cX && "PositronX[1]>=0" && cZ && cR && c20 && cGamma && cGammaMax && cPe && cN;
		break;
	default:
		cut = cX && cY && cR && c20 && cGamma && cGammaMax && cPe && cN;
	}
	tMC->Project(hMC->GetName(), str, cut);
	hMC->Sumw2();
	cv->Clear();
	hExp->SetLineColor(kBlack);
	hExp->SetLineWidth(3);
	hMC->SetLineColor(kBlue);
	hMC->SetLineWidth(1);
	hMC->Scale(hExp->Integral() / hMC->Integral());
	hExp->Draw();
	hMC->Draw("hist,same");
	TLegend *lg = new TLegend(0.5, 0.2, 0.65, 0.35);
	lg->AddEntry(hExp, "IBD", "LE");
	lg->AddEntry(hMC,  "MC",  "L");
	lg->Draw();
	cv->Update();
}

void RebinHist(TH1D* h1, TFile *out)
{
  // This rebinning method doesn't work
  // It needs to divide by the entries to make an average
  // Also, it doesn't include error weighting.
  cout<<"Rebinning histogram"<<endl;
  Int_t nBins = h1->GetNbinsX();
  TH1D *rebinHist = (TH1D*) h1->Clone();
  rebinHist->Rebin(nBins/2);

  // TFile out("Compare.root","update");
  TDirectory *dir = out->GetDirectory("Rebin");
  if(!dir) dir = out->mkdir("Rebin");
  dir->cd();

  TString histName = h1->GetName();
  histName += "Rebin";
  rebinHist->SetName(histName);
  rebinHist->SetTitle(histName);
  rebinHist->SetDirectory(0);
  rebinHist->Write(rebinHist->GetName(), TObject::kOverwrite);
  // out->Close();

}