-void MakeFinalSpectrum()\r
-{\r
- //-----------------------------------------------------------------------------\r
- // This macro takes the raw pi0 spectrum from LHC11a_FitResult_20130314.root,\r
- // correct it for feed down, then for efficiency,\r
- // adds all systematic errors and produces the production invariant spectrum\r
- //--\r
- // Last modification: 06.07.2012 by Yuri Kharlov\r
- //-----------------------------------------------------------------------------\r
-\r
- TFile *f = new TFile("LHC11a_FitResult_20130913.root");\r
- TH1D * nr1CB = (TH1D*)f->Get("Mix_CB_yr1") ;\r
- TH1D * nr1intCB = (TH1D*)f->Get("Mix_CB_yr1int") ;\r
- TH1D * nr2CB = (TH1D*)f->Get("Mix_CB_yr2") ;\r
- TH1D * nr2intCB = (TH1D*)f->Get("Mix_CB_yr2int") ;\r
-\r
- TH1D * nr1GS = (TH1D*)f->Get("Mix_yr1") ;\r
- TH1D * nr1intGS = (TH1D*)f->Get("Mix_yr1int") ;\r
- TH1D * nr2GS = (TH1D*)f->Get("Mix_yr2") ;\r
- TH1D * nr2intGS = (TH1D*)f->Get("Mix_yr2int") ;\r
-\r
- //Divide by bin width\r
- for(Int_t i=1;i<= nr1CB->GetNbinsX();i++){\r
- nr1CB ->SetBinContent(i,nr1CB->GetBinContent(i)/nr1CB->GetXaxis()->GetBinWidth(i)) ;\r
- nr1CB ->SetBinError (i,nr1CB->GetBinError(i)/nr1CB->GetXaxis()->GetBinWidth(i)) ;\r
- nr1intCB->SetBinContent(i,nr1intCB->GetBinContent(i)/nr1intCB->GetXaxis()->GetBinWidth(i)) ;\r
- nr1intCB->SetBinError (i,nr1intCB->GetBinError(i)/nr1intCB->GetXaxis()->GetBinWidth(i)) ;\r
- nr2CB ->SetBinContent(i,nr2CB->GetBinContent(i)/nr2CB->GetXaxis()->GetBinWidth(i)) ;\r
- nr2CB ->SetBinError (i,nr2CB->GetBinError(i)/nr2CB->GetXaxis()->GetBinWidth(i)) ;\r
- nr2intCB->SetBinContent(i,nr2intCB->GetBinContent(i)/nr2intCB->GetXaxis()->GetBinWidth(i)) ;\r
- nr2intCB->SetBinError (i,nr2intCB->GetBinError(i)/nr2intCB->GetXaxis()->GetBinWidth(i)) ;\r
-\r
- nr1GS ->SetBinContent(i,nr1GS->GetBinContent(i)/nr1GS->GetXaxis()->GetBinWidth(i)) ;\r
- nr1GS ->SetBinError (i,nr1GS->GetBinError(i)/nr1GS->GetXaxis()->GetBinWidth(i)) ;\r
- nr1intGS->SetBinContent(i,nr1intGS->GetBinContent(i)/nr1intGS->GetXaxis()->GetBinWidth(i)) ;\r
- nr1intGS->SetBinError (i,nr1intGS->GetBinError(i)/nr1intGS->GetXaxis()->GetBinWidth(i)) ;\r
- nr2GS ->SetBinContent(i,nr2GS->GetBinContent(i)/nr2GS->GetXaxis()->GetBinWidth(i)) ;\r
- nr2GS ->SetBinError (i,nr2GS->GetBinError(i)/nr2GS->GetXaxis()->GetBinWidth(i)) ;\r
- nr2intGS->SetBinContent(i,nr2intGS->GetBinContent(i)/nr2intGS->GetXaxis()->GetBinWidth(i)) ;\r
- nr2intGS->SetBinError (i,nr2intGS->GetBinError(i)/nr2intGS->GetXaxis()->GetBinWidth(i)) ;\r
- }\r
-\r
- // feed down correction\r
-\r
- TF1 *fKaonContaminationToPi0 = new TF1("kaonCont","1./(1.-1.33*1.2*exp(-2.95-0.16*x))",0.,30.);\r
- nr1GS ->Divide(fKaonContaminationToPi0);\r
- nr1intGS->Divide(fKaonContaminationToPi0);\r
- nr2GS ->Divide(fKaonContaminationToPi0);\r
- nr2intGS->Divide(fKaonContaminationToPi0);\r
- nr1CB ->Divide(fKaonContaminationToPi0);\r
- nr2CB ->Divide(fKaonContaminationToPi0);\r
- nr1intCB->Divide(fKaonContaminationToPi0);\r
- nr2intCB->Divide(fKaonContaminationToPi0);\r
-\r
- // SPD pileup correction\r
-\r
- TF1 *fSPDpileup = new TF1("SPDpileup","0.988",0.,30.);\r
- nr1GS ->Multiply(fSPDpileup);\r
- nr1intGS->Multiply(fSPDpileup);\r
- nr2GS ->Multiply(fSPDpileup);\r
- nr2intGS->Multiply(fSPDpileup);\r
- nr1CB ->Multiply(fSPDpileup);\r
- nr2CB ->Multiply(fSPDpileup);\r
- nr1intCB->Multiply(fSPDpileup);\r
- nr2intCB->Multiply(fSPDpileup);\r
-\r
- //correct for efficiency\r
- TFile *fEff = new TFile("Pi0_efficiency_LHC11a__20131029_Mall.root") ;\r
-\r
- TF1 * effGS=fEff->Get("eff_Pi0_Gaus_2760GeV") ;\r
- TF1 * effCB=fEff->Get("eff_Pi0_CB_2760GeV") ;\r
- effGS ->SetRange(0.,25.) ;\r
- effCB ->SetRange(0.,25.) ;\r
-\r
- nr1GS ->Divide(effGS) ;\r
- nr1intGS->Divide(effGS) ;\r
- nr2GS ->Divide(effGS) ;\r
- nr2intGS->Divide(effGS) ;\r
- nr1CB ->Divide(effCB) ;\r
- nr2CB ->Divide(effCB) ;\r
- nr1intCB->Divide(effCB) ;\r
- nr2intCB->Divide(effCB) ;\r
-\r
- //make 1/pt\r
- for(Int_t i=1;i<=nr1CB->GetNbinsX();i++){\r
- Double_t pt = TMath::TwoPi()*nr1CB->GetXaxis()->GetBinCenter(i);\r
- nr1CB ->SetBinContent(i,nr1CB ->GetBinContent(i)/pt) ;\r
- nr1CB ->SetBinError (i,nr1CB ->GetBinError(i) /pt) ;\r
- nr1intCB->SetBinContent(i,nr1intCB->GetBinContent(i)/pt) ;\r
- nr1intCB->SetBinError (i,nr1intCB->GetBinError(i) /pt) ;\r
- nr2CB ->SetBinContent(i,nr2CB ->GetBinContent(i)/pt) ;\r
- nr2CB ->SetBinError (i,nr2CB ->GetBinError(i) /pt) ;\r
- nr2intCB->SetBinContent(i,nr2intCB->GetBinContent(i)/pt) ;\r
- nr2intCB->SetBinError (i,nr2intCB->GetBinError(i) /pt) ;\r
-\r
- nr1GS ->SetBinContent(i,nr1GS ->GetBinContent(i)/pt) ;\r
- nr1GS ->SetBinError (i,nr1GS ->GetBinError(i) /pt) ;\r
- nr1intGS->SetBinContent(i,nr1intGS->GetBinContent(i)/pt) ;\r
- nr1intGS->SetBinError (i,nr1intGS->GetBinError(i) /pt) ;\r
- nr2GS ->SetBinContent(i,nr2GS ->GetBinContent(i)/pt) ;\r
- nr2GS ->SetBinError (i,nr2GS ->GetBinError(i) /pt) ;\r
- nr2intGS->SetBinContent(i,nr2intGS->GetBinContent(i)/pt) ;\r
- nr2intGS->SetBinError (i,nr2intGS->GetBinError(i) /pt) ;\r
- }\r
-\r
- //For the final spectrum we take average of fits\r
- //with numerical integration of entries in signal\r
- \r
- TH1D * hStat = (TH1D*)nr2intCB->Clone("hPi02760GeVStat") ;\r
- TH1D * hSys = (TH1D*)hStat ->Clone("hPi02760GeVSys") ;\r
- TH1D * hSys2 = (TH1D*)hStat ->Clone("hPi02760GeVSysTypeB") ;\r
- TH1D * hSys3 = (TH1D*)hStat ->Clone("hPi02760GeVSysTypeC") ;\r
- hStat->SetAxisRange(0.,14.9,"X");\r
- hSys ->SetAxisRange(0.,14.9,"X");\r
-\r
- //For systematic error estimate take largest deviation\r
- //of integrated yeilds (note, they are efficiency corrected)\r
- for(Int_t i=1;i<=nr1CB->GetNbinsX();i++){\r
- Double_t mean= hStat->GetBinContent(i) ;\r
- Double_t dev = TMath::Max(\r
- TMath::Max(TMath::Abs(nr1intCB->GetBinContent(i)-mean),\r
- TMath::Abs(nr2intCB->GetBinContent(i)-mean)),\r
- TMath::Max(TMath::Abs(nr1intGS->GetBinContent(i)-mean),\r
- TMath::Abs(nr2intGS->GetBinContent(i)-mean))\r
- );\r
- hSys ->SetBinError(i,dev) ;\r
- hSys2->SetBinError(i,dev) ;\r
- }\r
-\r
- //Add other sys errors\r
- TF1 * globalE = new TF1("globalE","1.-((x+1.354)/(x*1.002+1.354))^6.18 ",1.,30.) ; \r
- TF1 * conv = new TF1("conversion","0.035",0.,30.) ;\r
- TF1 * accept = new TF1("accept" ,"0.01" ,0.,30.) ;\r
- TF1 * pileup = new TF1("pileup" ,"0.004",0.,30.) ;\r
- TF1 * calib = new TF1("calib" ,"0.005",0.,30.) ;\r
- TF1 * modDiff = new TF1("modDiff" ,"0.04",0.,30.) ;\r
- // TF1 * modDiff = new TF1("modDiff" ,"16.9*exp(-4.5*x)+0.033",0.,30.) ;\r
- TF1 * tofCut = new TF1("tofCut" ,"0.0105" ,0.,30.) ;\r
-\r
- //Borya's estimate of non-linearity (found for pp @ 7 TeV)\r
- TF1 * nonlin= new TF1("nl","0.015+7.38*exp(-x/0.24)",0.,30.) ;\r
-\r
- //Draw sys errors\r
- TH1D * hRelSysRaw = (TH1D*)hSys->Clone("RelSysRaw") ;\r
- hRelSysRaw->SetTitle("Summary of systematic errors");\r
- for(Int_t i=1;i<=hSys->GetNbinsX();i++){\r
- Double_t mean= hSys->GetBinContent(i) ;\r
- Double_t a=hSys->GetBinError(i) ;\r
- if(mean>0)\r
- hRelSysRaw->SetBinContent(i,a/mean) ;\r
- else\r
- hRelSysRaw->SetBinContent(i,0.) ;\r
- hRelSysRaw->SetBinError(i,0.) ;\r
- }\r
-\r
- //Add errors in sys errors\r
-\r
- TH1D * hRelSysTot = (TH1D*)hSys->Clone("RelSys") ;\r
- hRelSysTot->SetTitle("Summary of systematic uncertainties");\r
-\r
- for(Int_t i=1;i<=hSys->GetNbinsX();i++){\r
- Double_t pt = hSys->GetXaxis()->GetBinCenter(i) ;\r
- Double_t mean = hSys->GetBinContent(i) ;\r
- Double_t a = hSys->GetBinError(i) ;\r
- // Double_t b = mean * hSysErrModules->GetBinContent(i) ;\r
- Double_t tot= mean*mean*nonlin ->Eval(pt)*nonlin ->Eval(pt) \r
- +mean*mean*conv ->Eval(pt)*conv ->Eval(pt)\r
- +mean*mean*accept ->Eval(pt)*accept ->Eval(pt)\r
- +mean*mean*pileup ->Eval(pt)*pileup ->Eval(pt)\r
- +mean*mean*calib ->Eval(pt)*calib ->Eval(pt)\r
- +mean*mean*modDiff->Eval(pt)*modDiff->Eval(pt)\r
- +mean*mean*tofCut ->Eval(pt)*tofCut ->Eval(pt)\r
- +mean*mean*globalE->Eval(pt)*globalE->Eval(pt); \r
- Double_t raa= mean*mean*nonlin->Eval(pt)*nonlin->Eval(pt) \r
- +mean*mean*pileup->Eval(pt)*pileup->Eval(pt)\r
- +mean*mean*calib ->Eval(pt)*calib ->Eval(pt); \r
- hSys3->SetBinError(i,TMath::Sqrt(tot)) ;\r
- // hSys->SetBinError(i,TMath::Sqrt(tot + a*a + b*b)) ;\r
- hSys2->SetBinError(i,TMath::Sqrt(raa + a*a)) ;\r
- hSys ->SetBinError(i,TMath::Sqrt(tot + a*a)) ;\r
- \r
- a = hSys->GetBinError(i) ;\r
- printf("i=%d, %g+-%g\n",i,mean,a);\r
- if(mean>0)\r
- hRelSysTot->SetBinContent(i,a/mean) ;\r
- else {\r
- hRelSysTot->SetBinContent(i,0.) ;\r
- hRelSysTot->SetBinError(i,0.) ;\r
- }\r
- }\r
-\r
- // TFile *fSysErrModules = TFile::Open("PHOS_sysErr_modules.root");\r
- // TH1D* hSysErrModules = (TH1D*)fSysErrModules->Get("hSyserr");\r
-\r
- gStyle->SetOptStat(0);\r
- TCanvas * c = new TCanvas("c","SysErrors") ;\r
- c->SetLogy();\r
- c->cd() ;\r
- gPad->SetRightMargin(0.02);\r
- gPad->SetTopMargin(0.07);\r
- hRelSysTot->SetAxisRange(0.8 ,11.9,"X");\r
- hRelSysTot->SetAxisRange(0.0031,0.45,"Y");\r
- hRelSysTot->GetYaxis()->SetMoreLogLabels();\r
- hRelSysTot->GetYaxis()->SetNoExponent();\r
- hRelSysTot->SetNdivisions(520,"X");\r
- hRelSysTot->SetLineColor(1) ;\r
- hRelSysTot->SetLineWidth(2) ;\r
- hRelSysRaw->SetLineColor(2) ;\r
- hRelSysRaw->SetLineWidth(2) ;\r
- globalE->SetLineColor(kGreen+2) ;\r
- nonlin ->SetLineColor(4) ;\r
- conv ->SetLineColor(6) ;\r
- accept ->SetLineColor(kOrange) ;\r
- pileup ->SetLineColor(42);\r
- calib ->SetLineColor(44);\r
- calib ->SetLineStyle(2);\r
- modDiff->SetLineColor(52);\r
- modDiff->SetLineStyle(2);\r
- tofCut ->SetLineColor(53);\r
- tofCut ->SetLineStyle(3);\r
- // hSysErrModules->SetLineColor(kOrange+2);\r
- // hSysErrModules->SetLineStyle(2);\r
- hRelSysTot->SetXTitle("p_{T} (GeV/c)") ;\r
- hRelSysTot->SetYTitle("Rel.syst.error") ;\r
- hRelSysTot->GetYaxis()->SetTitleOffset(1.2) ;\r
- hRelSysTot->Draw("hist") ;\r
- hRelSysRaw->Draw("h same") ;\r
- globalE->Draw("same") ; \r
- nonlin ->Draw("same") ;\r
- conv ->Draw("same") ;\r
- accept ->Draw("same") ;\r
- pileup ->Draw("same");\r
- calib ->Draw("same");\r
- modDiff->Draw("same");\r
- tofCut ->Draw("same");\r
- // hSysErrModules->Draw("same ][");\r
- TLegend * l = new TLegend(0.57,0.62,0.85,0.925) ;\r
- l->SetFillColor(kWhite);\r
- l->SetBorderSize(0);\r
- l->AddEntry(hRelSysTot,"Total uncertainty","l") ;\r
- l->AddEntry(hRelSysRaw,"Raw extraction","l") ;\r
- l->AddEntry(conv ,"Conversion","l") ;\r
- l->AddEntry(nonlin ,"Non-linearity","l") ;\r
- l->AddEntry(accept ,"Acceptance","l");\r
- l->AddEntry(pileup ,"Pileup","l");\r
- l->AddEntry(calib ,"Rel.calib.","l");\r
- l->AddEntry(modDiff,"Per-module yield","l");\r
- l->AddEntry(tofCut ,"Timing cut","l");\r
- // l->AddEntry(hSysErrModules,"Intermodule spectra","l");\r
- l->AddEntry(globalE,"Global E scale","l") ;\r
- l->Draw() ;\r
-\r
- c->Print("LHC11a_SysErrors.eps");\r
-\r
- hStat->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, stat.err.");\r
- hStat->SetXTitle("p_{T}, GeV/c");\r
- hStat->SetYTitle("1/N_{ev} 1/(2#pi p_{T}) d^{2}N/dydp_{T} (GeV^{-2}c^{2})");\r
- hSys ->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, total syst.err.");\r
- hSys ->SetXTitle("p_{T}, GeV/c");\r
- hSys ->SetYTitle("1/N_{ev} 1/(2#pi p_{T}) d^{2}N/dydp_{T} (GeV^{-2}c^{2})");\r
- hSys2->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, syst.err. for R_{AA}");\r
- hSys3->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, apparatus syst.err.");\r
- hSys ->SetFillColor(kBlue-10) ;\r
- hSys ->SetNdivisions(520,"X");\r
-\r
- hSysFinal = (TH1F*)hSys->Clone("hSysFinal");\r
- hSysFinal ->SetTitle("Production #pi^{0} yield, pp @ 2.76 TeV, 07.11.2013");\r
- hSysFinal ->GetYaxis()->SetTitleOffset(1.2);\r
- hStat->GetYaxis()->SetTitleOffset(1.2);\r
- hSysFinal ->SetAxisRange(0.8,9.9,"X");\r
-\r
- TCanvas *c2 = new TCanvas("c2","Production spectrum");\r
- gPad->SetRightMargin(0.02);\r
- gPad->SetTopMargin(0.07);\r
- gPad->SetGridx();\r
- gPad->SetGridy();\r
- c2->SetLogy();\r
- hSysFinal->Draw("E2") ;\r
- hStat->SetMarkerStyle(20) ;\r
- hStat->SetMarkerColor(4) ;\r
- hStat->SetLineColor(4) ;\r
- hStat->Draw("same") ;\r
- \r
- // //Apply bin width correction\r
- // TH1D * hBWcorr = BinWidthCorrection(hStat) ;\r
- // hSys->Divide(hBWcorr) ;\r
- // hSys2->Divide(hBWcorr) ;\r
- // hSys3->Divide(hBWcorr) ;\r
- // hStat->Divide(hBWcorr) ;\r
-\r
- c2->Print("LHC11a_pi0Spectrum.eps");\r
-\r
- TF1 *tsallis = FitTsallis(hStat,hSys);\r
- \r
- TFile fout("PHOS_pi0_2760eV_noBinWidthCorr_20131112.root","recreate") ;\r
- hSys ->Write() ;\r
- hSys2 ->Write() ;\r
- hSys3 ->Write() ;\r
- hStat ->Write() ;\r
- tsallis->Write() ;\r
- fout.Close() ;\r
-\r
-\r
-}\r
-\r
-//-----------------------------------------------------------------------------\r
-TF1 *FitTsallis(TH1D* hStat, TH1D* hSys)\r
-{\r
- TF1 * fit = new TF1("Tsalis","[0]/2./3.1415*([2]-1.)*([2]-2.)/([2]*[1]*([2]*[1]+0.135*([2]-2.)))*(1.+(sqrt(x*x+0.135*0.135)-0.135)/([2]*[1]))^-[2]",0.5,25.) ;\r
- fit->SetParameters(10.,0.2,8.) ;\r
- fit->SetLineColor(kBlack) ;\r
- fit->SetLineWidth(2) ;\r
- fit->SetParName(0,"dN/dy") ;\r
- fit->SetParName(1,"T") ;\r
- fit->SetParName(2,"n") ;\r
- \r
- TH1D * hSysRatio = (TH1D*)hSys ->Clone("RatioSys") ;\r
- TH1D * hStatRatio = (TH1D*)hStat->Clone("RatioStat") ;\r
-\r
- TH1D * hsum = (TH1D*)hStat->Clone("sum") ;\r
- for(Int_t i=1; i<=hsum->GetNbinsX();i++){\r
- Double_t a=hSys->GetBinError(i) ;\r
- Double_t b=hStat->GetBinError(i) ;\r
- hsum->SetBinError(i,TMath::Sqrt(a*a+b*b)) ;\r
- }\r
- hsum->SetStats(1) ;\r
- hsum->Fit(fit,"Q") ;\r
- Double_t meanPt=fit->Moment(1,1.,10.) ;\r
- printf("<pt>=%f \n",meanPt) ;\r
-\r
- hSysRatio ->Divide(fit) ;\r
- hStatRatio->Divide(fit) ;\r
-\r
- return fit;\r
-}\r
-//-----------------------------------------------------------------------------\r
-TH1D * BinWidthCorrection(TH1D * h){\r
- //We apply bin width a-la PHENIX \r
- //Use Tsalis fit to calculate shift in y direction\r
- \r
- TF1 * fit = new TF1("hag","[0]*(([2]*[1]+sqrt(x*x+0.135*0.135))/([2]*[1]+0.135))^-[2]",0.5,25.) ;\r
- fit->SetParameters(10.,0.2,8.) ;\r
- TCanvas * corr = new TCanvas("BWcorr","Bin width correction") ;\r
- Int_t col[6]={kRed,kOrange,kMagenta,kGreen,kCyan,kBlue} ;\r
- TH1D * hcorr[20] ;\r
- char key[55] ;\r
- Double_t rMax=10 ;\r
- Int_t iteration=0 ;\r
- TH1D * htmp = (TH1D*)h->Clone("tmp") ;\r
- while(iteration<6){\r
- printf(" Iteration %d: rMax=%f \n",iteration, rMax) ;\r
- htmp->Fit(fit,"N") ;\r
- sprintf(key,"Ineration%d",iteration) ;\r
- hcorr[iteration]=(TH1D*)h->Clone(key);\r
- rMax= 0; \r
- for(Int_t i=1;i<=h->GetNbinsX();i++){\r
- Double_t a=h->GetXaxis()->GetBinLowEdge(i) ;\r
- Double_t b=h->GetXaxis()->GetBinUpEdge(i) ;\r
- Double_t r=fit->Integral(a,b)/(b-a)/fit->Eval(0.5*(a+b)) ;\r
- hcorr[iteration]->SetBinContent(i,r) ;\r
- hcorr[iteration]->SetBinError(i,0.) ;\r
- if(rMax<r)rMax=r ;\r
- }\r
- delete htmp ;\r
- htmp = (TH1D*)h->Clone("tmp") ;\r
- htmp->Divide(hcorr[iteration]) ;\r
- corr->cd() ;\r
- hcorr[iteration]->SetLineColor(col[iteration]);\r
- if(iteration==0)\r
- hcorr[iteration]->Draw() ;\r
- else\r
- hcorr[iteration]->Draw("same") ;\r
- corr->Update() ;\r
- iteration++ ;\r
- } \r
-\r
- hcorr[5]->SetTitle("Bin-width correction for #pi^{0} spectrum");\r
- hcorr[5]->SetYTitle("Bin-width corrected / uncorrected");\r
- TFile fout("PHOS_pi0_7TeV_BinWidthCorrection.root","recreate") ;\r
- hcorr[5]->Write();\r
- fout.Close();\r
-\r
- return hcorr[5] ;\r
-}\r
+void MakeFinalSpectrum()
+{
+ //-----------------------------------------------------------------------------
+ // This macro takes the raw pi0 spectrum from LHC11a_FitResult_20130314.root,
+ // correct it for feed down, then for efficiency,
+ // adds all systematic errors and produces the production invariant spectrum
+ //--
+ // Last modification: 06.07.2012 by Yuri Kharlov
+ //-----------------------------------------------------------------------------
+
+ TFile *f = new TFile("LHC11a_FitResult_20130913.root");
+ TH1D * nr1CB = (TH1D*)f->Get("Mix_CB_yr1") ;
+ TH1D * nr1intCB = (TH1D*)f->Get("Mix_CB_yr1int") ;
+ TH1D * nr2CB = (TH1D*)f->Get("Mix_CB_yr2") ;
+ TH1D * nr2intCB = (TH1D*)f->Get("Mix_CB_yr2int") ;
+
+ TH1D * nr1GS = (TH1D*)f->Get("Mix_yr1") ;
+ TH1D * nr1intGS = (TH1D*)f->Get("Mix_yr1int") ;
+ TH1D * nr2GS = (TH1D*)f->Get("Mix_yr2") ;
+ TH1D * nr2intGS = (TH1D*)f->Get("Mix_yr2int") ;
+
+ //Divide by bin width
+ for(Int_t i=1;i<= nr1CB->GetNbinsX();i++){
+ nr1CB ->SetBinContent(i,nr1CB->GetBinContent(i)/nr1CB->GetXaxis()->GetBinWidth(i)) ;
+ nr1CB ->SetBinError (i,nr1CB->GetBinError(i)/nr1CB->GetXaxis()->GetBinWidth(i)) ;
+ nr1intCB->SetBinContent(i,nr1intCB->GetBinContent(i)/nr1intCB->GetXaxis()->GetBinWidth(i)) ;
+ nr1intCB->SetBinError (i,nr1intCB->GetBinError(i)/nr1intCB->GetXaxis()->GetBinWidth(i)) ;
+ nr2CB ->SetBinContent(i,nr2CB->GetBinContent(i)/nr2CB->GetXaxis()->GetBinWidth(i)) ;
+ nr2CB ->SetBinError (i,nr2CB->GetBinError(i)/nr2CB->GetXaxis()->GetBinWidth(i)) ;
+ nr2intCB->SetBinContent(i,nr2intCB->GetBinContent(i)/nr2intCB->GetXaxis()->GetBinWidth(i)) ;
+ nr2intCB->SetBinError (i,nr2intCB->GetBinError(i)/nr2intCB->GetXaxis()->GetBinWidth(i)) ;
+
+ nr1GS ->SetBinContent(i,nr1GS->GetBinContent(i)/nr1GS->GetXaxis()->GetBinWidth(i)) ;
+ nr1GS ->SetBinError (i,nr1GS->GetBinError(i)/nr1GS->GetXaxis()->GetBinWidth(i)) ;
+ nr1intGS->SetBinContent(i,nr1intGS->GetBinContent(i)/nr1intGS->GetXaxis()->GetBinWidth(i)) ;
+ nr1intGS->SetBinError (i,nr1intGS->GetBinError(i)/nr1intGS->GetXaxis()->GetBinWidth(i)) ;
+ nr2GS ->SetBinContent(i,nr2GS->GetBinContent(i)/nr2GS->GetXaxis()->GetBinWidth(i)) ;
+ nr2GS ->SetBinError (i,nr2GS->GetBinError(i)/nr2GS->GetXaxis()->GetBinWidth(i)) ;
+ nr2intGS->SetBinContent(i,nr2intGS->GetBinContent(i)/nr2intGS->GetXaxis()->GetBinWidth(i)) ;
+ nr2intGS->SetBinError (i,nr2intGS->GetBinError(i)/nr2intGS->GetXaxis()->GetBinWidth(i)) ;
+ }
+
+ // feed down correction
+
+ TF1 *fKaonContaminationToPi0 = new TF1("kaonCont","1./(1.-1.33*1.2*exp(-2.95-0.16*x))",0.,30.);
+ nr1GS ->Divide(fKaonContaminationToPi0);
+ nr1intGS->Divide(fKaonContaminationToPi0);
+ nr2GS ->Divide(fKaonContaminationToPi0);
+ nr2intGS->Divide(fKaonContaminationToPi0);
+ nr1CB ->Divide(fKaonContaminationToPi0);
+ nr2CB ->Divide(fKaonContaminationToPi0);
+ nr1intCB->Divide(fKaonContaminationToPi0);
+ nr2intCB->Divide(fKaonContaminationToPi0);
+
+ // SPD pileup correction
+
+ TF1 *fSPDpileup = new TF1("SPDpileup","0.988",0.,30.);
+ nr1GS ->Multiply(fSPDpileup);
+ nr1intGS->Multiply(fSPDpileup);
+ nr2GS ->Multiply(fSPDpileup);
+ nr2intGS->Multiply(fSPDpileup);
+ nr1CB ->Multiply(fSPDpileup);
+ nr2CB ->Multiply(fSPDpileup);
+ nr1intCB->Multiply(fSPDpileup);
+ nr2intCB->Multiply(fSPDpileup);
+
+ //correct for efficiency
+ TFile *fEff = new TFile("Pi0_efficiency_LHC11a__20131029_Mall.root") ;
+
+ TF1 * effGS=fEff->Get("eff_Pi0_Gaus_2760GeV") ;
+ TF1 * effCB=fEff->Get("eff_Pi0_CB_2760GeV") ;
+ effGS ->SetRange(0.,25.) ;
+ effCB ->SetRange(0.,25.) ;
+
+ nr1GS ->Divide(effGS) ;
+ nr1intGS->Divide(effGS) ;
+ nr2GS ->Divide(effGS) ;
+ nr2intGS->Divide(effGS) ;
+ nr1CB ->Divide(effCB) ;
+ nr2CB ->Divide(effCB) ;
+ nr1intCB->Divide(effCB) ;
+ nr2intCB->Divide(effCB) ;
+
+ //make 1/pt
+ for(Int_t i=1;i<=nr1CB->GetNbinsX();i++){
+ Double_t pt = TMath::TwoPi()*nr1CB->GetXaxis()->GetBinCenter(i);
+ nr1CB ->SetBinContent(i,nr1CB ->GetBinContent(i)/pt) ;
+ nr1CB ->SetBinError (i,nr1CB ->GetBinError(i) /pt) ;
+ nr1intCB->SetBinContent(i,nr1intCB->GetBinContent(i)/pt) ;
+ nr1intCB->SetBinError (i,nr1intCB->GetBinError(i) /pt) ;
+ nr2CB ->SetBinContent(i,nr2CB ->GetBinContent(i)/pt) ;
+ nr2CB ->SetBinError (i,nr2CB ->GetBinError(i) /pt) ;
+ nr2intCB->SetBinContent(i,nr2intCB->GetBinContent(i)/pt) ;
+ nr2intCB->SetBinError (i,nr2intCB->GetBinError(i) /pt) ;
+
+ nr1GS ->SetBinContent(i,nr1GS ->GetBinContent(i)/pt) ;
+ nr1GS ->SetBinError (i,nr1GS ->GetBinError(i) /pt) ;
+ nr1intGS->SetBinContent(i,nr1intGS->GetBinContent(i)/pt) ;
+ nr1intGS->SetBinError (i,nr1intGS->GetBinError(i) /pt) ;
+ nr2GS ->SetBinContent(i,nr2GS ->GetBinContent(i)/pt) ;
+ nr2GS ->SetBinError (i,nr2GS ->GetBinError(i) /pt) ;
+ nr2intGS->SetBinContent(i,nr2intGS->GetBinContent(i)/pt) ;
+ nr2intGS->SetBinError (i,nr2intGS->GetBinError(i) /pt) ;
+ }
+
+ //For the final spectrum we take average of fits
+ //with numerical integration of entries in signal
+
+ TH1D * hStat = (TH1D*)nr2intCB->Clone("hPi02760GeVStat") ;
+ TH1D * hSys = (TH1D*)hStat ->Clone("hPi02760GeVSys") ;
+ TH1D * hSys2 = (TH1D*)hStat ->Clone("hPi02760GeVSysTypeB") ;
+ TH1D * hSys3 = (TH1D*)hStat ->Clone("hPi02760GeVSysTypeC") ;
+ hStat->SetAxisRange(0.,14.9,"X");
+ hSys ->SetAxisRange(0.,14.9,"X");
+
+ //For systematic error estimate take largest deviation
+ //of integrated yeilds (note, they are efficiency corrected)
+ for(Int_t i=1;i<=nr1CB->GetNbinsX();i++){
+ Double_t mean= hStat->GetBinContent(i) ;
+ Double_t dev = TMath::Max(
+ TMath::Max(TMath::Abs(nr1intCB->GetBinContent(i)-mean),
+ TMath::Abs(nr2intCB->GetBinContent(i)-mean)),
+ TMath::Max(TMath::Abs(nr1intGS->GetBinContent(i)-mean),
+ TMath::Abs(nr2intGS->GetBinContent(i)-mean))
+ );
+ hSys ->SetBinError(i,dev) ;
+ hSys2->SetBinError(i,dev) ;
+ }
+
+ //Add other sys errors
+ TF1 * globalE = new TF1("globalE","1.-((x+1.354)/(x*1.002+1.354))^6.18 ",1.,30.) ;
+ TF1 * conv = new TF1("conversion","0.035",0.,30.) ;
+ TF1 * accept = new TF1("accept" ,"0.01" ,0.,30.) ;
+ TF1 * pileup = new TF1("pileup" ,"0.004",0.,30.) ;
+ TF1 * calib = new TF1("calib" ,"0.005",0.,30.) ;
+ TF1 * modDiff = new TF1("modDiff" ,"0.04",0.,30.) ;
+ // TF1 * modDiff = new TF1("modDiff" ,"16.9*exp(-4.5*x)+0.033",0.,30.) ;
+ TF1 * tofCut = new TF1("tofCut" ,"0.0105" ,0.,30.) ;
+
+ //Borya's estimate of non-linearity (found for pp @ 7 TeV)
+ TF1 * nonlin= new TF1("nl","0.015+7.38*exp(-x/0.24)",0.,30.) ;
+
+ //Draw sys errors
+ TH1D * hRelSysRaw = (TH1D*)hSys->Clone("RelSysRaw") ;
+ hRelSysRaw->SetTitle("Summary of systematic errors");
+ for(Int_t i=1;i<=hSys->GetNbinsX();i++){
+ Double_t mean= hSys->GetBinContent(i) ;
+ Double_t a=hSys->GetBinError(i) ;
+ if(mean>0)
+ hRelSysRaw->SetBinContent(i,a/mean) ;
+ else
+ hRelSysRaw->SetBinContent(i,0.) ;
+ hRelSysRaw->SetBinError(i,0.) ;
+ }
+
+ //Add errors in sys errors
+
+ TH1D * hRelSysTot = (TH1D*)hSys->Clone("RelSys") ;
+ hRelSysTot->SetTitle("Summary of systematic uncertainties");
+
+ for(Int_t i=1;i<=hSys->GetNbinsX();i++){
+ Double_t pt = hSys->GetXaxis()->GetBinCenter(i) ;
+ Double_t mean = hSys->GetBinContent(i) ;
+ Double_t a = hSys->GetBinError(i) ;
+ // Double_t b = mean * hSysErrModules->GetBinContent(i) ;
+ Double_t tot= mean*mean*nonlin ->Eval(pt)*nonlin ->Eval(pt)
+ +mean*mean*conv ->Eval(pt)*conv ->Eval(pt)
+ +mean*mean*accept ->Eval(pt)*accept ->Eval(pt)
+ +mean*mean*pileup ->Eval(pt)*pileup ->Eval(pt)
+ +mean*mean*calib ->Eval(pt)*calib ->Eval(pt)
+ +mean*mean*modDiff->Eval(pt)*modDiff->Eval(pt)
+ +mean*mean*tofCut ->Eval(pt)*tofCut ->Eval(pt)
+ +mean*mean*globalE->Eval(pt)*globalE->Eval(pt);
+ Double_t raa= mean*mean*nonlin->Eval(pt)*nonlin->Eval(pt)
+ +mean*mean*pileup->Eval(pt)*pileup->Eval(pt)
+ +mean*mean*calib ->Eval(pt)*calib ->Eval(pt);
+ hSys3->SetBinError(i,TMath::Sqrt(tot)) ;
+ // hSys->SetBinError(i,TMath::Sqrt(tot + a*a + b*b)) ;
+ hSys2->SetBinError(i,TMath::Sqrt(raa + a*a)) ;
+ hSys ->SetBinError(i,TMath::Sqrt(tot + a*a)) ;
+
+ a = hSys->GetBinError(i) ;
+ printf("i=%d, %g+-%g\n",i,mean,a);
+ if(mean>0)
+ hRelSysTot->SetBinContent(i,a/mean) ;
+ else {
+ hRelSysTot->SetBinContent(i,0.) ;
+ hRelSysTot->SetBinError(i,0.) ;
+ }
+ }
+
+ // TFile *fSysErrModules = TFile::Open("PHOS_sysErr_modules.root");
+ // TH1D* hSysErrModules = (TH1D*)fSysErrModules->Get("hSyserr");
+
+ gStyle->SetOptStat(0);
+ TCanvas * c = new TCanvas("c","SysErrors") ;
+ c->SetLogy();
+ c->cd() ;
+ gPad->SetRightMargin(0.02);
+ gPad->SetTopMargin(0.07);
+ hRelSysTot->SetAxisRange(0.8 ,11.9,"X");
+ hRelSysTot->SetAxisRange(0.0031,0.45,"Y");
+ hRelSysTot->GetYaxis()->SetMoreLogLabels();
+ hRelSysTot->GetYaxis()->SetNoExponent();
+ hRelSysTot->SetNdivisions(520,"X");
+ hRelSysTot->SetLineColor(1) ;
+ hRelSysTot->SetLineWidth(2) ;
+ hRelSysRaw->SetLineColor(2) ;
+ hRelSysRaw->SetLineWidth(2) ;
+ globalE->SetLineColor(kGreen+2) ;
+ nonlin ->SetLineColor(4) ;
+ conv ->SetLineColor(6) ;
+ accept ->SetLineColor(kOrange) ;
+ pileup ->SetLineColor(42);
+ calib ->SetLineColor(44);
+ calib ->SetLineStyle(2);
+ modDiff->SetLineColor(52);
+ modDiff->SetLineStyle(2);
+ tofCut ->SetLineColor(53);
+ tofCut ->SetLineStyle(3);
+ // hSysErrModules->SetLineColor(kOrange+2);
+ // hSysErrModules->SetLineStyle(2);
+ hRelSysTot->SetXTitle("p_{T} (GeV/c)") ;
+ hRelSysTot->SetYTitle("Rel.syst.error") ;
+ hRelSysTot->GetYaxis()->SetTitleOffset(1.2) ;
+ hRelSysTot->Draw("hist") ;
+ hRelSysRaw->Draw("h same") ;
+ globalE->Draw("same") ;
+ nonlin ->Draw("same") ;
+ conv ->Draw("same") ;
+ accept ->Draw("same") ;
+ pileup ->Draw("same");
+ calib ->Draw("same");
+ modDiff->Draw("same");
+ tofCut ->Draw("same");
+ // hSysErrModules->Draw("same ][");
+ TLegend * l = new TLegend(0.57,0.62,0.85,0.925) ;
+ l->SetFillColor(kWhite);
+ l->SetBorderSize(0);
+ l->AddEntry(hRelSysTot,"Total uncertainty","l") ;
+ l->AddEntry(hRelSysRaw,"Raw extraction","l") ;
+ l->AddEntry(conv ,"Conversion","l") ;
+ l->AddEntry(nonlin ,"Non-linearity","l") ;
+ l->AddEntry(accept ,"Acceptance","l");
+ l->AddEntry(pileup ,"Pileup","l");
+ l->AddEntry(calib ,"Rel.calib.","l");
+ l->AddEntry(modDiff,"Per-module yield","l");
+ l->AddEntry(tofCut ,"Timing cut","l");
+ // l->AddEntry(hSysErrModules,"Intermodule spectra","l");
+ l->AddEntry(globalE,"Global E scale","l") ;
+ l->Draw() ;
+
+ c->Print("LHC11a_SysErrors.eps");
+
+ hStat->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, stat.err.");
+ hStat->SetXTitle("p_{T}, GeV/c");
+ hStat->SetYTitle("1/N_{ev} 1/(2#pi p_{T}) d^{2}N/dydp_{T} (GeV^{-2}c^{2})");
+ hSys ->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, total syst.err.");
+ hSys ->SetXTitle("p_{T}, GeV/c");
+ hSys ->SetYTitle("1/N_{ev} 1/(2#pi p_{T}) d^{2}N/dydp_{T} (GeV^{-2}c^{2})");
+ hSys2->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, syst.err. for R_{AA}");
+ hSys3->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, apparatus syst.err.");
+ hSys ->SetFillColor(kBlue-10) ;
+ hSys ->SetNdivisions(520,"X");
+
+ hSysFinal = (TH1F*)hSys->Clone("hSysFinal");
+ hSysFinal ->SetTitle("Production #pi^{0} yield, pp @ 2.76 TeV, 07.11.2013");
+ hSysFinal ->GetYaxis()->SetTitleOffset(1.2);
+ hStat->GetYaxis()->SetTitleOffset(1.2);
+ hSysFinal ->SetAxisRange(0.8,9.9,"X");
+
+ TCanvas *c2 = new TCanvas("c2","Production spectrum");
+ gPad->SetRightMargin(0.02);
+ gPad->SetTopMargin(0.07);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ c2->SetLogy();
+ hSysFinal->Draw("E2") ;
+ hStat->SetMarkerStyle(20) ;
+ hStat->SetMarkerColor(4) ;
+ hStat->SetLineColor(4) ;
+ hStat->Draw("same") ;
+
+ // //Apply bin width correction
+ // TH1D * hBWcorr = BinWidthCorrection(hStat) ;
+ // hSys->Divide(hBWcorr) ;
+ // hSys2->Divide(hBWcorr) ;
+ // hSys3->Divide(hBWcorr) ;
+ // hStat->Divide(hBWcorr) ;
+
+ c2->Print("LHC11a_pi0Spectrum.eps");
+
+ TF1 *tsallis = FitTsallis(hStat,hSys);
+
+ TFile fout("PHOS_pi0_2760eV_noBinWidthCorr_20131112.root","recreate") ;
+ hSys ->Write() ;
+ hSys2 ->Write() ;
+ hSys3 ->Write() ;
+ hStat ->Write() ;
+ tsallis->Write() ;
+ fout.Close() ;
+
+
+}
+
+//-----------------------------------------------------------------------------
+TF1 *FitTsallis(TH1D* hStat, TH1D* hSys)
+{
+ TF1 * fit = new TF1("Tsalis","[0]/2./3.1415*([2]-1.)*([2]-2.)/([2]*[1]*([2]*[1]+0.135*([2]-2.)))*(1.+(sqrt(x*x+0.135*0.135)-0.135)/([2]*[1]))^-[2]",0.5,25.) ;
+ fit->SetParameters(10.,0.2,8.) ;
+ fit->SetLineColor(kBlack) ;
+ fit->SetLineWidth(2) ;
+ fit->SetParName(0,"dN/dy") ;
+ fit->SetParName(1,"T") ;
+ fit->SetParName(2,"n") ;
+
+ TH1D * hSysRatio = (TH1D*)hSys ->Clone("RatioSys") ;
+ TH1D * hStatRatio = (TH1D*)hStat->Clone("RatioStat") ;
+
+ TH1D * hsum = (TH1D*)hStat->Clone("sum") ;
+ for(Int_t i=1; i<=hsum->GetNbinsX();i++){
+ Double_t a=hSys->GetBinError(i) ;
+ Double_t b=hStat->GetBinError(i) ;
+ hsum->SetBinError(i,TMath::Sqrt(a*a+b*b)) ;
+ }
+ hsum->SetStats(1) ;
+ hsum->Fit(fit,"Q") ;
+ Double_t meanPt=fit->Moment(1,1.,10.) ;
+ printf("<pt>=%f \n",meanPt) ;
+
+ hSysRatio ->Divide(fit) ;
+ hStatRatio->Divide(fit) ;
+
+ return fit;
+}
+//-----------------------------------------------------------------------------
+TH1D * BinWidthCorrection(TH1D * h){
+ //We apply bin width a-la PHENIX
+ //Use Tsalis fit to calculate shift in y direction
+
+ TF1 * fit = new TF1("hag","[0]*(([2]*[1]+sqrt(x*x+0.135*0.135))/([2]*[1]+0.135))^-[2]",0.5,25.) ;
+ fit->SetParameters(10.,0.2,8.) ;
+ TCanvas * corr = new TCanvas("BWcorr","Bin width correction") ;
+ Int_t col[6]={kRed,kOrange,kMagenta,kGreen,kCyan,kBlue} ;
+ TH1D * hcorr[20] ;
+ char key[55] ;
+ Double_t rMax=10 ;
+ Int_t iteration=0 ;
+ TH1D * htmp = (TH1D*)h->Clone("tmp") ;
+ while(iteration<6){
+ printf(" Iteration %d: rMax=%f \n",iteration, rMax) ;
+ htmp->Fit(fit,"N") ;
+ sprintf(key,"Ineration%d",iteration) ;
+ hcorr[iteration]=(TH1D*)h->Clone(key);
+ rMax= 0;
+ for(Int_t i=1;i<=h->GetNbinsX();i++){
+ Double_t a=h->GetXaxis()->GetBinLowEdge(i) ;
+ Double_t b=h->GetXaxis()->GetBinUpEdge(i) ;
+ Double_t r=fit->Integral(a,b)/(b-a)/fit->Eval(0.5*(a+b)) ;
+ hcorr[iteration]->SetBinContent(i,r) ;
+ hcorr[iteration]->SetBinError(i,0.) ;
+ if(rMax<r)rMax=r ;
+ }
+ delete htmp ;
+ htmp = (TH1D*)h->Clone("tmp") ;
+ htmp->Divide(hcorr[iteration]) ;
+ corr->cd() ;
+ hcorr[iteration]->SetLineColor(col[iteration]);
+ if(iteration==0)
+ hcorr[iteration]->Draw() ;
+ else
+ hcorr[iteration]->Draw("same") ;
+ corr->Update() ;
+ iteration++ ;
+ }
+
+ hcorr[5]->SetTitle("Bin-width correction for #pi^{0} spectrum");
+ hcorr[5]->SetYTitle("Bin-width corrected / uncorrected");
+ TFile fout("PHOS_pi0_7TeV_BinWidthCorrection.root","recreate") ;
+ hcorr[5]->Write();
+ fout.Close();
+
+ return hcorr[5] ;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff