From 47151f26b5f8625cfcf820b1061c949e99fc0387 Mon Sep 17 00:00:00 2001 From: cnattras Date: Wed, 15 May 2013 22:57:19 +0000 Subject: [PATCH] Implementing multiplicity dependent efficiency corrections --- PWGLF/totEt/AliAnalysisEt.cxx | 4 +- PWGLF/totEt/AliAnalysisEt.h | 2 +- PWGLF/totEt/AliAnalysisEtMonteCarlo.cxx | 4 +- PWGLF/totEt/AliAnalysisEtRecEffCorrection.cxx | 14 +- PWGLF/totEt/AliAnalysisEtRecEffCorrection.h | 13 +- PWGLF/totEt/AliAnalysisEtReconstructed.cxx | 10 +- PWGLF/totEt/AliAnalysisEtReconstructed.h | 2 +- .../totEt/AliAnalysisEtReconstructedEmcal.cxx | 4 +- PWGLF/totEt/AliAnalysisEtReconstructedEmcal.h | 2 +- PWGLF/totEt/macros/calculateCorrections.C | 102 ++++++++- PWGLF/totEt/macros/emEt/PlotGammaEfficiency.C | 55 ++++- .../macros/emEt/PlotGammaEfficiencyErrors.C | 216 ++++++++++++++++++ 12 files changed, 408 insertions(+), 20 deletions(-) create mode 100644 PWGLF/totEt/macros/emEt/PlotGammaEfficiencyErrors.C diff --git a/PWGLF/totEt/AliAnalysisEt.cxx b/PWGLF/totEt/AliAnalysisEt.cxx index 7e42d329ec3..fdd2bbfb178 100644 --- a/PWGLF/totEt/AliAnalysisEt.cxx +++ b/PWGLF/totEt/AliAnalysisEt.cxx @@ -701,12 +701,12 @@ Int_t AliAnalysisEt::ReadCorrections(TString filename) return -1; } -Double_t AliAnalysisEt::CorrectForReconstructionEfficiency(const AliESDCaloCluster& cluster) +Double_t AliAnalysisEt::CorrectForReconstructionEfficiency(const AliESDCaloCluster& cluster, int mult) { Float_t pos[3]; cluster.GetPosition(pos); TVector3 cp(pos); - Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E()); + Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E(), mult); //std::cout << "Original energy: " << cluster.E() << ", corrected energy: " << corrEnergy << std::endl; return TMath::Sin(cp.Theta())*corrEnergy; diff --git a/PWGLF/totEt/AliAnalysisEt.h b/PWGLF/totEt/AliAnalysisEt.h index f1758d7fd70..d579e52aa7a 100644 --- a/PWGLF/totEt/AliAnalysisEt.h +++ b/PWGLF/totEt/AliAnalysisEt.h @@ -136,7 +136,7 @@ protected: //AliAnalysisEtCuts *fCuts; // keeper of basic cuts // Return corrected cluster E_T - Double_t CorrectForReconstructionEfficiency(const AliESDCaloCluster &cluster); + Double_t CorrectForReconstructionEfficiency(const AliESDCaloCluster &cluster,int mult = 0); // Track matching (hadrdonic contamination) corrections AliAnalysisEtTrackMatchCorrections *fTmCorrections; diff --git a/PWGLF/totEt/AliAnalysisEtMonteCarlo.cxx b/PWGLF/totEt/AliAnalysisEtMonteCarlo.cxx index 32fa47f2094..80747c15e81 100644 --- a/PWGLF/totEt/AliAnalysisEtMonteCarlo.cxx +++ b/PWGLF/totEt/AliAnalysisEtMonteCarlo.cxx @@ -766,7 +766,7 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2) } fCutFlow->Fill(cf++); if(!fSelector->PassDistanceToBadChannelCut(*caloCluster)) continue; - Double_t clEt = CorrectForReconstructionEfficiency(*caloCluster); + Double_t clEt = CorrectForReconstructionEfficiency(*caloCluster,fClusterMult); // if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl; if(!fSelector->PassMinEnergyCut(*caloCluster)) continue; @@ -1076,7 +1076,7 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2) caloCluster->GetPosition(pos); TVector3 cp(pos); Double_t clEt = caloCluster->E()*TMath::Sin(cp.Theta()); - Double_t clEtCorr = CorrectForReconstructionEfficiency(*caloCluster); + Double_t clEtCorr = CorrectForReconstructionEfficiency(*caloCluster,fClusterMult); for(int l=0;l=etCuts[l]){ //cout<<", "<="<Eval(energy); } +Double_t AliAnalysisEtRecEffCorrection::CorrectedEnergy(Double_t energy, int mult) +{ + if(fRecoEff) return 1.0/ReconstructionEfficiency(energy, mult); + return 1.0; + +} + +Double_t AliAnalysisEtRecEffCorrection::ReconstructionEfficiency(float energy, int mult) const { + Double_t eff = 1.0; + if(fRecoEff) eff = fRecoEff->GetBinContent(fRecoEff->GetXaxis()->FindBin(energy),fRecoEff->GetXaxis()->FindBin(mult)); + //cout <<"eff "<TrackMatchingEfficiency(track->Pt(),fClusterMult); - Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster); + Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,fClusterMult); nChargedHadronsEtMeasured+= TMath::Sin(cp.Theta())*effCorrEt; nChargedHadronsEtTotal+= 1/fTmCorrections->TrackMatchingEfficiency(track->Pt(),fClusterMult) *effCorrEt; fHistMatchedTracksEvspTvsMult->Fill(track->P(),TMath::Sin(cp.Theta())*cluster->E(),fClusterMult); - fHistMatchedTracksEvspTvsMultEffCorr->Fill(track->P(),CorrectForReconstructionEfficiency(*cluster),fClusterMult); + fHistMatchedTracksEvspTvsMultEffCorr->Fill(track->P(),CorrectForReconstructionEfficiency(*cluster,fClusterMult),fClusterMult); const Double_t *pidWeights = track->PID(); Double_t maxpidweight = 0; @@ -283,7 +283,7 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev) fClusterEnergy->Fill(cluster->E()); fClusterEt->Fill(TMath::Sin(p2.Theta())*cluster->E()); - Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster); + Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,fClusterMult); fTotNeutralEt += effCorrEt; nominalRawEt += effCorrEt; nonlinHighRawEt += effCorrEt*GetCorrectionModification(*cluster,1,0); @@ -508,12 +508,12 @@ void AliAnalysisEtReconstructed::CreateHistograms() fHistNominalEffLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5); } -Double_t AliAnalysisEtReconstructed::ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr) +Double_t AliAnalysisEtReconstructed::ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t mult) { Float_t pos[3]; cluster.GetPosition(pos); TVector3 cp(pos); - Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E()); + Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E(),mult); Double_t factorNonLin = GetCorrectionModification(cluster, nonLinCorr,effCorr); diff --git a/PWGLF/totEt/AliAnalysisEtReconstructed.h b/PWGLF/totEt/AliAnalysisEtReconstructed.h index c7ba4268b80..096efe73580 100644 --- a/PWGLF/totEt/AliAnalysisEtReconstructed.h +++ b/PWGLF/totEt/AliAnalysisEtReconstructed.h @@ -84,7 +84,7 @@ protected: TH2D *fHistNominalEffHighEt;//Total ET from clusters with high bound on reconstruction efficiency and nominal nonlinearity correction vs centrality TH2D *fHistNominalEffLowEt;//Total ET from clusters with low bound on reconstruction efficiency and nominal nonlinearity correction vs centrality - Double_t ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr);//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr 0 = nominal 1 = high -1 = low + Double_t ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t mult);//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr 0 = nominal 1 = high -1 = low private: diff --git a/PWGLF/totEt/AliAnalysisEtReconstructedEmcal.cxx b/PWGLF/totEt/AliAnalysisEtReconstructedEmcal.cxx index 611bd36b37f..a10d93e6a0e 100644 --- a/PWGLF/totEt/AliAnalysisEtReconstructedEmcal.cxx +++ b/PWGLF/totEt/AliAnalysisEtReconstructedEmcal.cxx @@ -53,9 +53,9 @@ void AliAnalysisEtReconstructedEmcal::CreateHistograms() } -Double_t AliAnalysisEtReconstructedEmcal::GetCorrectionModification(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr){//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr 0 = nominal 1 = high -1 = low +Double_t AliAnalysisEtReconstructedEmcal::GetCorrectionModification(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t mult){//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr 0 = nominal 1 = high -1 = low Double_t factor = 1.0; - Double_t E = fReCorrections->CorrectedEnergy(cluster.E()); + Double_t E = fReCorrections->CorrectedEnergy(cluster.E(),mult); if(nonLinCorr!=0){ Double_t p0 = 9.90780e-01; Double_t p1 = 1.61503e-01; diff --git a/PWGLF/totEt/AliAnalysisEtReconstructedEmcal.h b/PWGLF/totEt/AliAnalysisEtReconstructedEmcal.h index 0cffc4d6a65..7591951eb07 100644 --- a/PWGLF/totEt/AliAnalysisEtReconstructedEmcal.h +++ b/PWGLF/totEt/AliAnalysisEtReconstructedEmcal.h @@ -27,7 +27,7 @@ public: virtual bool TrackHitsCalorimeter(AliVParticle *track, Double_t magField); - virtual Double_t GetCorrectionModification(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr);//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr 0 = nominal 1 = high -1 = low + virtual Double_t GetCorrectionModification(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t mult);//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr 0 = nominal 1 = high -1 = low private: ClassDef(AliAnalysisEtReconstructedEmcal, 1); diff --git a/PWGLF/totEt/macros/calculateCorrections.C b/PWGLF/totEt/macros/calculateCorrections.C index 4834561e243..c1a87007503 100644 --- a/PWGLF/totEt/macros/calculateCorrections.C +++ b/PWGLF/totEt/macros/calculateCorrections.C @@ -17,6 +17,8 @@ TCanvas *c1 = 0; +TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) ; +TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) ; //creates an empty set of correction factors for debugging purposes TF1 * generateRecEffFunction(Double_t p0, Double_t p1); @@ -201,11 +203,21 @@ int calculateCorrections(TString filename="rootFiles/LHC11a10a_bis/Et.ESD.simPbP Double_t meanGamma = 0.51/(p0 + p1*0.51); Double_t meanSecondary = meanGamma; - AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(("TmCorrections"+detector).Data(), fitcharged, fitneutral, fitgamma, fitsecondary, + TH2F *fHistHadronDepositsAllMult = l->FindObject("fHistHadronDepositsAllMult"); + TH2F *fHistHadronDepositsRecoMult = l->FindObject("fHistHadronDepositsRecoMult"); + TH2F *eff2D = (TH2F*) bayneseffdiv2D(fHistHadronDepositsRecoMult,fHistHadronDepositsAllMult,"eff2D"); + //cor->SetReconstructionEfficiency(eff2D); + + TH2F *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedMult"); + TH2F *fHistGammasFoundMult = l->FindObject("fHistGammasFoundMult"); + TH2F *gammaEff2D = (TH2F*) bayneseffdiv2D(fHistGammasFoundMult,fHistGammasGeneratedMult,"gammaEff2D"); + + AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(("TmCorrections"+detector).Data(), fitcharged, fitneutral, fitgamma, fitsecondary, *eff2D, meanCharged, meanNeutral, meanGamma, meanSecondary ); + TF1 *func = generateRecEffFunction(p0, p1); - AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(("ReCorrections"+detector).Data(), *func, 1000); + AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(("ReCorrections"+detector).Data(), *func,*gammaEff2D, 1000); TFile *outfile = TFile::Open("calocorrections.root","RECREATE"); cor->Write(); recor->Write(); @@ -221,3 +233,89 @@ TF1* generateRecEffFunction(Double_t p0, Double_t p1) f->SetParameters(params); return f; } + +TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) +{ + if(!numerator){ + cerr<<"Error: numerator does not exist!"<Clone(name); + Int_t nbins = numerator->GetNbinsX(); + for (Int_t ibin=0; ibin<= nbins+1; ++ibin) { + Double_t numeratorVal = numerator->GetBinContent(ibin); + Double_t denominatorVal = denominator->GetBinContent(ibin); + // Check if the errors are right or the thing is scaled + Double_t numeratorValErr = numerator->GetBinError(ibin); + if (!(numeratorValErr==0. || numeratorVal ==0.) ) { + Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal; + numeratorVal /= rescale; + } + Double_t denominatorValErr = denominator->GetBinError(ibin); + if (!(denominatorValErr==0. || denominatorVal==0. )) { + Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal; + denominatorVal /= rescale; + } + Double_t quotient = 0.; + if (denominatorVal!=0.) { + quotient = numeratorVal/denominatorVal; + } + Double_t quotientErr=0; + quotientErr = TMath::Sqrt( + (numeratorVal+1.0)/(denominatorVal+2.0)* + ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0))); + result->SetBinContent(ibin,quotient); + result->SetBinError(ibin,quotientErr); + //cout<<"Setting bin "<Clone(name); + Int_t nbinsX = numerator->GetNbinsX(); + Int_t nbinsY = numerator->GetNbinsY(); + for (Int_t ibin=0; ibin<= nbinsX+1; ++ibin) { + for (Int_t jbin=0; jbin<= nbinsY+1; ++jbin) { + Double_t numeratorVal = numerator->GetBinContent(ibin,jbin); + Double_t denominatorVal = denominator->GetBinContent(ibin,jbin); + // Check if the errors are right or the thing is scaled + Double_t numeratorValErr = numerator->GetBinError(ibin,jbin); + if (!(numeratorValErr==0. || numeratorVal ==0.) ) { + Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal; + numeratorVal /= rescale; + } + Double_t denominatorValErr = denominator->GetBinError(ibin,jbin); + if (!(denominatorValErr==0. || denominatorVal==0. )) { + Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal; + denominatorVal /= rescale; + } + Double_t quotient = 0.; + if (denominatorVal!=0.) { + quotient = numeratorVal/denominatorVal; + } + Double_t quotientErr=0; + quotientErr = TMath::Sqrt( + (numeratorVal+1.0)/(denominatorVal+2.0)* + ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0))); + result->SetBinContent(ibin,jbin,quotient); + result->SetBinError(ibin,jbin,quotientErr); + //cout<<"Setting bin "<SetOptTitle(0); gStyle->SetOptStat(0); @@ -10,6 +12,8 @@ void PlotGammaEfficiency(TString filename="rootFiles/LHC11a4_bis/Et.ESD.simPbPb. TH1F *fHistGammasGenerated = l->FindObject("fHistGammasGenerated"); TH1F *fHistGammasFound = l->FindObject("fHistGammasFound"); TH1F *hEfficiency = bayneseffdiv(fHistGammasFound,fHistGammasGenerated,"Efficiency"); + TH2F *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedMult"); + TH2F *fHistGammasFoundMult = l->FindObject("fHistGammasFoundMult"); hEfficiency->GetXaxis()->SetTitle("Energy"); hEfficiency->GetYaxis()->SetTitle("efficiency"); hEfficiency->GetXaxis()->SetRange(1,hEfficiency->FindBin(3)); @@ -27,10 +31,49 @@ void PlotGammaEfficiency(TString filename="rootFiles/LHC11a4_bis/Et.ESD.simPbPb. hEfficiency->Draw(); c->SaveAs("/tmp/GammaEfficiency.png"); + TCanvas *c1 = new TCanvas("c1","c1",600,400); + c1->SetTopMargin(0.02); + c1->SetRightMargin(0.02); + c1->SetBorderSize(0); + c1->SetFillColor(0); + c1->SetFillColor(0); + c1->SetBorderMode(0); + c1->SetFrameFillColor(0); + c1->SetFrameBorderMode(0); + TH1D *fHistGammasGeneratedCent[11] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL, NULL}; + TH1D *fHistGammasFoundCent[11] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL, NULL}; + TH1D *fEff[11] = {NULL,NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL,NULL, NULL}; + int colors[] = {TColor::kRed,TColor::kRed, TColor::kOrange-3, TColor::kOrange-3, TColor::kGreen+3, + TColor::kGreen+3, TColor::kBlue, TColor::kBlue, TColor::kViolet, TColor::kViolet, + TColor::kMagenta+3}; + int markers[] = {20,24,21,25,22, 26,23,32,33,27, 29}; + TLegend *leg1 = new TLegend(0.224832,0.153226,0.39094,0.647849); + leg1->SetFillStyle(0); + leg1->SetFillColor(0); + leg1->SetBorderSize(0); + leg1->SetTextSize(0.03); + leg1->SetTextSize(0.0456989); + for(int i=0;i<10;i++){ + fHistGammasGeneratedCent[i] = fHistGammasGeneratedMult->ProjectionX(Form("All%i",i),i+1,i+2); + fHistGammasFoundCent[i] = fHistGammasFoundMult->ProjectionX(Form("Reco%i",i),i+1,i+2); + fEff[i] = (TH1D*) bayneseffdiv(fHistGammasFoundCent[i],fHistGammasGeneratedCent[i],Form("eff%i",i)); + SetStyles(fEff[i],markers[i],colors[i],"energy","efficiency",1.0); + leg1->AddEntry(fEff[i],Form("%iSetMaximum(0.8); + fEff[i]->GetXaxis()->SetRange(1,fEff[i]->FindBin(3.0)); + fEff[i]->Draw(); + } + else{ fEff[i]->Draw("same");} + + } + //leg1->Draw(); +// fHistGammasFound->Draw(); + } -TH1F* bayneseffdiv(TH1F* numerator, TH1F* denominator,Char_t* name) +TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) { if(!numerator){ cerr<<"Error: numerator does not exist!"<Sumw2(); + histo->SetMarkerStyle(marker); + histo->SetMarkerColor(color); + histo->SetLineColor(color); + histo->GetXaxis()->SetTitle(xtitle); + histo->GetYaxis()->SetTitle(ytitle); +} diff --git a/PWGLF/totEt/macros/emEt/PlotGammaEfficiencyErrors.C b/PWGLF/totEt/macros/emEt/PlotGammaEfficiencyErrors.C new file mode 100644 index 00000000000..cf9c41cbb13 --- /dev/null +++ b/PWGLF/totEt/macros/emEt/PlotGammaEfficiencyErrors.C @@ -0,0 +1,216 @@ +void SetStyles(TH1 *histo,int marker, int color,char *xtitle = "p", char *ytitle = "E"); +TFile *junkfile = NULL; + +TFile *fNominal = NULL; +TFile *fLow = NULL; +TFile *fHigh = NULL; + + +TH1F* bayneseffdiv(TH1F* numerator, TH1F* denominator,Char_t* name); +TH1F *GetHisto(Bool_t isPHOS,int var){ +// TString myfilename = filename; +// if(var==1){myfilename = filenameLow;} +// if(var==2){myfilename = filenameHigh;} + TString listname = ""; + TString histoName = ""; + TFile *f = NULL; + switch(var){ + case 1: + histoName = "Efficiency1"; + listname = "List1"; + f = fLow; + break; + case 2: + histoName = "Efficiency2"; + listname = "List2"; + f = fHigh; + break; + default: + histoName = "Efficiency0"; + listname = "List0"; + f = fNominal; + } + //cout<<"My file name "<GetName()<cd(); + //TList *l = (TList*)f->Get("out1"); + f->ls(); + //l->SetName(listname.Data()); + //cout<<"List name "<GetName()<FindObject("fHistGammasGenerated")->Clone("GeneratedLHC11b1a"); + //TH1F *fHistGammasFound = l->FindObject("fHistGammasFound")->Clone("FoundLHC11b1a"); + TH1F *fHistGammasGeneratedLHC11b1a = (TH1F*)f->Get("fHistGammasGenerated"); + TH1F *fHistGammasFoundLHC11b1a = (TH1F*)f->Get("fHistGammasFound"); + junkfile->cd(); + TH1F *hEfficiency = bayneseffdiv(fHistGammasFoundLHC11b1a,fHistGammasGeneratedLHC11b1a,histoName.Data()); +// delete fHistGammasGenerated; +// delete fHistGammasFound; +// delete l; + hEfficiency->Rebin(2); + hEfficiency->Scale(0.5); + hEfficiency->SetName(histoName.Data()); + //cout<<"Histo name "<GetName()<<" entries "<GetEntries()<<" "<GetEntries()<GetXaxis()->SetRange(1,hEfficiency->FindBin(3)); + //f->Close(); + //delete f; + return hEfficiency; +} +void PlotGammaEfficiencyErrors(Bool_t isPHOS = kFALSE){ + + gStyle->SetOptTitle(0); + gStyle->SetOptStat(0); + gStyle->SetOptFit(0); + + +TString det = "EMCal"; +if(isPHOS) det = "PHOS"; +// TString run = ""; +// switch(runnum){ +// case 1: +// run = "118506";//bad file +// break; +// case 2: +// run = "121040"; +// break; +// case 3: +// run = "118518";//low stats +// break; +// case 4: +// run = "121039";//really low stats +// break; +// case 5: +// run = "118558";//bad file +// break; +// } +// TString run = "118558";//bad file - adding files didn't work. +//TString run = "121040"; +//TString run = "121039"; + TString outfilename = "junk.root"; + junkfile = new TFile(outfilename.Data(),"RECREATE"); + TString filename = "rootFiles/LHC11b1a/Efficiency."+det+".LHC11b1a.Combined.root"; + TString filenameHigh = "rootFiles/LHC11b1b/Efficiency."+det+".LHC11b1b.Combined.root"; + TString filenameLow = "rootFiles/LHC11b1c/Efficiency."+det+".LHC11b1c.Combined.root"; +// TString filenameHigh = "rootFiles/LHC11b1b/Et.ESD.sim."+det+".LHC11b1b.Run"+run+".root"; +// TString filename = "rootFiles/LHC11b1a/Et.ESD.sim."+det+".LHC11b1a.Run"+run+".root"; +// TString filenameLow = "rootFiles/LHC11b1c/Et.ESD.sim."+det+".LHC11b1c.Run"+run+".root"; + +fNominal = new TFile(filename.Data()); +fLow = new TFile(filenameLow.Data()); +fHigh = new TFile(filenameHigh.Data()); + + + hEfficiency = GetHisto(isPHOS,0); + SetStyles(hEfficiency,20,1,"Energy","efficiency"); + hEfficiencyLow = GetHisto(isPHOS,1); + SetStyles(hEfficiencyLow,23,2,"Energy","efficiency"); + hEfficiencyHigh = GetHisto(isPHOS,2); + SetStyles(hEfficiencyHigh,22,4,"Energy","efficiency"); + + junkfile->Write(); + //junkfile->Close(); + TCanvas *c = new TCanvas("c","c",600,400); + c->SetTopMargin(0.02); + c->SetRightMargin(0.02); + c->SetBorderSize(0); + c->SetFillColor(0); + c->SetFillColor(0); + c->SetBorderMode(0); + c->SetFrameFillColor(0); + c->SetFrameBorderMode(0); + + hEfficiency->Draw(); + hEfficiencyLow->Draw("same"); + hEfficiencyHigh->Draw("same"); + TString pic1name = "/tmp/GammaEfficiency"+det+".eps"; + c->SaveAs(pic1name.Data()); + + TH1F *hRatioLow = (TH1F*) hEfficiencyLow->Clone("hEfficiencyLowRatio"); + TH1F *hRatioHigh = (TH1F*) hEfficiencyHigh->Clone("hEfficiencyHighRatio"); + hRatioHigh->Divide(hEfficiency); + hRatioLow->Divide(hEfficiency); + + TCanvas *c1 = new TCanvas("c1","c1",600,400); + c1->SetTopMargin(0.02); + c1->SetRightMargin(0.02); + c1->SetBorderSize(0); + c1->SetFillColor(0); + c1->SetFillColor(0); + c1->SetBorderMode(0); + c1->SetFrameFillColor(0); + c1->SetFrameBorderMode(0); + hRatioLow->SetMaximum(1.2); + hRatioLow->SetMinimum(0.8); + hRatioLow->GetYaxis()->SetTitle("efficiency/nominal efficiency"); + TF1 *funcHigh = new TF1("funcHigh","[0]",0,3); + TF1 *funcLow = new TF1("funcLow","[0]",0,3); + funcHigh->SetLineColor(hRatioHigh->GetMarkerColor()); + funcLow->SetLineColor(hRatioLow->GetMarkerColor()); + TLegend *leg = new TLegend(0.144295,0.163978,0.234899,0.346774); + leg->SetFillStyle(0); + leg->SetFillColor(0); + leg->SetBorderSize(0); + leg->SetTextSize(0.03); + leg->AddEntry(hRatioLow,"Low material budget"); + leg->AddEntry(hRatioHigh,"High material budget"); + leg->SetTextSize(0.0537634); + hRatioLow->Fit(funcLow); + hRatioHigh->Fit(funcHigh); + hRatioLow->Draw(); + hRatioHigh->Draw("same"); + leg->Draw(); + pic1name = "/tmp/GammaEfficiencyError"+det+".eps"; + c1->SaveAs(pic1name.Data()); + +} + + +TH1F* bayneseffdiv(TH1F* numerator, TH1F* denominator,Char_t* name) +{ + if(!numerator){ + cerr<<"Error: numerator does not exist!"<Clone(name); + Int_t nbins = numerator->GetNbinsX(); + for (Int_t ibin=0; ibin<= nbins+1; ++ibin) { + Double_t numeratorVal = numerator->GetBinContent(ibin); + Double_t denominatorVal = denominator->GetBinContent(ibin); + // Check if the errors are right or the thing is scaled + Double_t numeratorValErr = numerator->GetBinError(ibin); + if (!(numeratorValErr==0. || numeratorVal ==0.) ) { + Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal; + numeratorVal /= rescale; + } + Double_t denominatorValErr = denominator->GetBinError(ibin); + if (!(denominatorValErr==0. || denominatorVal==0. )) { + Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal; + denominatorVal /= rescale; + } + Double_t quotient = 0.; + if (denominatorVal!=0.) { + quotient = numeratorVal/denominatorVal; + } + Double_t quotientErr=0; + quotientErr = TMath::Sqrt( + (numeratorVal+1.0)/(denominatorVal+2.0)* + ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0))); + result->SetBinContent(ibin,quotient); + result->SetBinError(ibin,quotientErr); + //cout<<"Setting bin "<SetName(name); + return result; +} + +void SetStyles(TH1 *histo,int marker, int color,char *xtitle, char *ytitle){ + histo->SetMarkerStyle(marker); + histo->SetMarkerColor(color); + histo->SetLineColor(color); + histo->GetXaxis()->SetTitle(xtitle); + histo->GetYaxis()->SetTitle(ytitle); + histo->Sumw2(); +} -- 2.43.0