8d9cb18035f1e1f72faf59ad2ad349f0bcc737ad
[u/mrichter/AliRoot.git] / PWG4 / macros / electrons / plotMCRates.C
1 /////////////////////////////////////////////////\r
2 //\r
3 // Macro for plotting MC rates of electrons\r
4 // for the EMCAL PPR\r
5 //\r
6 // J.L. Klay (Cal Poly)\r
7 //\r
8 /////////////////////////////////////////////////\r
9 \r
10 TH1F* all;\r
11 TH1F* bele;\r
12 TH1F* cele;\r
13 TH1F* candb;\r
14 TH1F* conv;\r
15 TH1F* dalitz;\r
16 TH1F* wz;\r
17 TH1F* other;\r
18 TH1F* mchad;\r
19 TH1F* sige;\r
20 TH1F* bkge;\r
21 TH1F* walle;\r
22 TH1F* hije;\r
23 \r
24 TLegend* leg;\r
25 \r
26 void plotMCRates(char* hijfname = "scale/histosLHC08d6.root",\r
27                  char* jjfname = "scale/histosscaledLHC09b2ESD.root",\r
28                  char* bfname = "scale/histosscaledLHC09b4.root",\r
29                  char* wfname = "scale/histosWBoson.root") {\r
30 \r
31   //For HIJING need to divide by the number of events, which we\r
32   //can get from the file and do when we perform scaling\r
33   double hijscale = 0.05*(1.E6)*0.5*7700; //0-5% * seconds*lumi*PbPb x-section\r
34   //For bjet and jet-jet events\r
35   double pyscale = (1.E6)*0.5*208*208*100/360; //seconds*lumi*Pb*Pb*acceptance\r
36   double bscale = pyscale*0.10; //Branching ratio for forced\r
37                                 //semi-leptonic decays\r
38   double wscale = pyscale*6.29e-05; //JLK: This is temporary X-sec\r
39                                      //info from 2-29 GeV bin until we\r
40                                      //get pyxsec files; also need to\r
41                                      //divide by nTrials.  For now,\r
42                                      //use nEvt\r
43   \r
44   TFile* hijfile = new TFile(hijfname);\r
45   if(!hijfile) { printf("NO HIJING FILE\n"); return; }\r
46   TList* hijlist = (TList*)hijfile->Get("histos");\r
47   TH2F* hijmcele = (TH2F*)histos->FindObject("AnaElectron_hPtMCElectron");\r
48   TH1F* hijmchad = (TH1F*)histos->FindObject("AnaElectron_hPtMCHadron");\r
49   TH1F* hijmult = (TH1F*)histos->FindObject("AnaElectron_hRefMult");\r
50   Int_t nEvt = hijmult->GetEntries();\r
51   if(nEvt == 0) { printf("NO HIJING EVENTS\n"); return; }\r
52   hijmcele->Scale(hijscale/nEvt);\r
53   hijmchad->Scale(hijscale/nEvt);\r
54 \r
55   TFile* jjfile = new TFile(jjfname);\r
56   if(!jjfile) { printf("NO JET-JET FILE\n"); return; }\r
57   TH2F* jjmcele = (TH2F*)jjfile->Get("AnaElectron_hPtMCElectronScaled");\r
58   TH1F* jjmchad = (TH1F*)jjfile->Get("AnaElectron_hPtMCHadronScaled");\r
59   TH1F* jjmult = (TH1F*)jjfile->Get("AnaElectron_hRefMultScaled");\r
60   Int_t nEvtJJ = jjmult->GetEntries();\r
61   jjmcele->Scale(pyscale);\r
62   jjmchad->Scale(pyscale);\r
63 \r
64   TFile* bfile = new TFile(bfname);\r
65   if(!bfile) { printf("NO B-JET FILE\n"); return; }\r
66   TH2F* bmcele = (TH2F*)bfile->Get("AnaElectron_hPtMCElectronScaled");\r
67   TH1F* bmchad = (TH1F*)bfile->Get("AnaElectron_hPtMCHadronScaled");\r
68   TH1F* bmult = (TH1F*)bfile->Get("AnaElectron_hRefMultScaled");\r
69   Int_t nEvtB = bmult->GetEntries();\r
70   bmcele->Scale(bscale);\r
71   bmchad->Scale(bscale);\r
72 \r
73   TFile* wfile = new TFile(wfname);\r
74   if(!wfile) { printf("NO W-BOSON FILE\n"); return; }\r
75   TList* wlist = (TList*)wfile->Get("histos");\r
76   TH2F* wmcele = (TH2F*)histos->FindObject("AnaElectron_hPtMCElectron");\r
77   TH1F* wmchad = (TH1F*)histos->FindObject("AnaElectron_hPtMCHadron");\r
78   TH1F* wmult = (TH1F*)histos->FindObject("AnaElectron_hRefMult");\r
79   Int_t nEvtW = wmult->GetEntries();\r
80   wmcele->Scale(wscale/nEvtW);\r
81   wmchad->Scale(wscale/nEvtW);\r
82 \r
83   printf("Event statistics: %d (HIJING)  %d (JET-JET)  %d (B-JET)  %d (W-Boson)\n",nEvt,nEvtJJ,nEvtB,nEvtW);\r
84 \r
85   TH2F* combined = (TH2F*)hijmcele->Clone();\r
86   combined->Add(jjmcele);\r
87   combined->Add(bmcele);\r
88   combined->SetTitle("MC electrons in Pb+Pb in EMCAL acceptance");\r
89   combined->SetName("CombinedMCEle");\r
90   combined->SetXTitle("p_T (GeV/c)");\r
91 \r
92   mchad = (TH1F*)hijmchad->Clone();\r
93   mchad->Add(jjmchad);\r
94   mchad->Add(bmchad);\r
95   mchad->Add(wmchad);\r
96   mchad->SetTitle("MC hadrons in Pb+Pb in EMCAL acceptance");\r
97   mchad->SetName("CombinedMCHad");\r
98   mchad->SetXTitle("p_T (GeV/c)");\r
99 \r
100   all = (TH1F*)combined->ProjectionX("all",1,1);\r
101   bele = (TH1F*)combined->ProjectionX("b",2,2);\r
102   cele = (TH1F*)combined->ProjectionX("c",3,3);\r
103   candb = (TH1F*)combined->ProjectionX("candb",4,4);\r
104   conv = (TH1F*)combined->ProjectionX("conv",5,5);\r
105   dalitz = (TH1F*)combined->ProjectionX("dalitz",6,6);\r
106   other = (TH1F*)combined->ProjectionX("other",8,8);\r
107 \r
108   wz = (TH1F*)wmcele->ProjectionX("wz",7,7);\r
109   \r
110   all->Add(wz); //because it had to be done separately\r
111 \r
112   //For comparing contributions\r
113   walle = (TH1F*)wmcele->ProjectionX("walle",7,7);\r
114   sige = (TH1F*)bmcele->ProjectionX("sige",1,1);\r
115   bkge = (TH1F*)jjmcele->ProjectionX("bkge",1,1);\r
116   hije = (TH1F*)hijmcele->ProjectionX("hije",1,1);\r
117 \r
118   double myscale = 1.; //we already scaled them\r
119   ScaleAndConfigure(all,myscale,kBlack,kFALSE);\r
120   ScaleAndConfigure(bele,myscale,kRed,kFALSE);\r
121   ScaleAndConfigure(sige,myscale,kRed,kFALSE);\r
122   ScaleAndConfigure(cele,myscale,kBlue,kFALSE);\r
123   ScaleAndConfigure(candb,myscale,kViolet,kFALSE);\r
124   ScaleAndConfigure(conv,myscale,kOrange-3,kFALSE);\r
125   ScaleAndConfigure(bkge,myscale,kOrange-3,kFALSE);\r
126   ScaleAndConfigure(dalitz,myscale,kGreen-3,kFALSE);\r
127   ScaleAndConfigure(wz,myscale,kOrange-7,kFALSE);\r
128   ScaleAndConfigure(walle,myscale,kOrange-7,kFALSE);\r
129   ScaleAndConfigure(mchad,myscale,kGreen+2,kFALSE);\r
130   ScaleAndConfigure(hije,myscale,kGreen+2,kFALSE);\r
131 \r
132   gStyle->SetOptStat(0);\r
133   //drawXSRates();\r
134   drawAnnualYields();\r
135   drawPtCutRates();\r
136   drawHadEleRatios();\r
137   drawSigBkg();\r
138 \r
139 }\r
140 \r
141 void ScaleAndConfigure(TH1F* hist,Double_t scale, Int_t color,Bool_t keepErr)\r
142 {\r
143   hist->Scale(scale);\r
144   hist->SetLineColor(color);\r
145   hist->SetLineWidth(2);\r
146   if(keepErr == kFALSE) {\r
147     //remove the error bars - useful for MC rates\r
148     for(Int_t i = 1; i < hist->GetNbinsX(); i++) {\r
149       hist->SetBinError(i,0.);\r
150     }\r
151   }\r
152 }\r
153 \r
154 void drawAnnualYields() {\r
155 \r
156   TCanvas* crates = new TCanvas();\r
157   crates->cd();\r
158   gPad->SetLogy();\r
159   all->SetXTitle("p_{T} (GeV/c)");\r
160   all->SetYTitle("Annual yield");\r
161   all->GetYaxis()->SetRangeUser(1.E1,3.E9);\r
162   all->GetXaxis()->SetRangeUser(0.,50.);\r
163   all->Draw();\r
164   bele->Draw("same");  \r
165   cele->Draw("same");  \r
166   candb->Draw("same");  \r
167   conv->Draw("same");  \r
168   dalitz->Draw("same");  \r
169   wz->Draw("same");\r
170 \r
171   leg = new TLegend(0.6,0.6,0.9,0.9);\r
172   leg->SetTextSize(leg->GetTextSize()*1.2);\r
173   leg->AddEntry(all,"All MC electrons","l");\r
174   leg->AddEntry(bele,"Bottom e","l");\r
175   leg->AddEntry(cele,"Charm e","l");\r
176   leg->AddEntry(candb,"B-->C e","l");\r
177   leg->AddEntry(dalitz,"Dalitz e","l");\r
178   leg->AddEntry(conv,"Conversion e","l");\r
179   leg->AddEntry(wz,"W Boson e","l");\r
180   leg->Draw();\r
181 \r
182 }\r
183 \r
184 void drawSigBkg() {\r
185 \r
186   TCanvas* csigbkg = new TCanvas();\r
187   csigbkg->cd();\r
188   gPad->SetLogy();\r
189   all->SetXTitle("p_{T} (GeV/c)");\r
190   all->SetYTitle("Annual yield");\r
191   all->GetYaxis()->SetRangeUser(1.E1,3.E9);\r
192   all->GetXaxis()->SetRangeUser(0.,50.);\r
193   all->Draw();\r
194   sige->Draw("same");  \r
195   bkge->Draw("same");  \r
196   hije->Draw("same");  \r
197   walle->Draw("same");\r
198 \r
199   TLegend* leg1 = new TLegend(0.6,0.6,0.9,0.9);\r
200   leg1->SetTextSize(leg->GetTextSize()*1.2);\r
201   leg1->AddEntry(all,"All MC electrons","l");\r
202   leg1->AddEntry(sige,"B-Jet Events","l");\r
203   leg1->AddEntry(hije,"Pb+Pb Underlying Event","l");\r
204   leg1->AddEntry(bkge,"Jet-Jet Events","l");\r
205   leg1->AddEntry(walle,"W-decay Events","l");\r
206   leg1->Draw();\r
207 \r
208 }\r
209 \r
210 void drawPtCutRates() {\r
211 \r
212   TCanvas* cptcut = new TCanvas();\r
213   cptcut->cd();\r
214   gPad->SetLogy();\r
215   TH1F* alleptcut = GetPtCutHisto(all);\r
216   TH1F* beleptcut = GetPtCutHisto(bele);\r
217   TH1F* celeptcut = GetPtCutHisto(cele);\r
218   TH1F* cbeleptcut = GetPtCutHisto(candb);\r
219   TH1F* dalitzptcut = GetPtCutHisto(dalitz);\r
220   TH1F* convptcut = GetPtCutHisto(conv);\r
221   TH1F* wzptcut = GetPtCutHisto(wz);\r
222   alleptcut->GetXaxis()->SetRangeUser(0,50);\r
223   alleptcut->GetYaxis()->SetRangeUser(100,alleptcut->GetMaximum()*2);\r
224   alleptcut->SetXTitle("p_{T}^{cut} (GeV/c)");\r
225   alleptcut->SetYTitle("Annual Yield in EMCAL for p_{T}>p_{T}^{cut}");\r
226   alleptcut->SetTitle("Annual electron yield in Pb+Pb for p_{T}>p_{T}^{cut}");\r
227   alleptcut->Draw();\r
228   beleptcut->Draw("same");\r
229   celeptcut->Draw("same");\r
230   cbeleptcut->Draw("same");\r
231   dalitzptcut->Draw("same");\r
232   convptcut->Draw("same");\r
233   wzptcut->Draw("same");\r
234   leg->Draw();\r
235 \r
236 }\r
237 \r
238 void drawHadEleRatios() {\r
239 \r
240   TCanvas* ceh = new TCanvas();\r
241   ceh->cd();\r
242   gPad->SetLogy();\r
243   gStyle->SetOptStat(0);\r
244   TH1F* allratio = (TH1F*)all->Clone();\r
245   TH1F* behratio = (TH1F*)bele->Clone();\r
246   allratio->SetTitle("Pb+Pb, 5.5 TeV");\r
247   allratio->SetXTitle("p_{T} (GeV/c)");\r
248   allratio->SetYTitle("Hadrons/Electrons");\r
249   for(Int_t i = 1; i < all->GetNbinsX(); i++) {\r
250     Double_t vale = all->GetBinContent(i);\r
251     Double_t valb = bele->GetBinContent(i);\r
252     Double_t valh = mchad->GetBinContent(i);\r
253     //printf("pT %.2f, Hadron %.1f, Electron %.1f, B-electron %.1f\n",all->GetBinCenter(i),valh,vale,valb);\r
254     if(vale>0) allratio->SetBinContent(i,valh/vale);\r
255     else allratio->SetBinContent(i,0.);\r
256 \r
257     if(valb>0) behratio->SetBinContent(i,valh/valb);\r
258     else behratio->SetBinContent(i,0.);\r
259 \r
260     allratio->SetBinError(i,0.);\r
261     behratio->SetBinError(i,0.);\r
262   }\r
263   allratio->Rebin();\r
264   behratio->Rebin();\r
265   allratio->GetYaxis()->SetRangeUser(1,10000);\r
266   allratio->GetXaxis()->SetRangeUser(0,50);\r
267   behratio->GetXaxis()->SetRangeUser(0,50);\r
268   allratio->SetMarkerStyle(20);\r
269   behratio->SetMarkerStyle(24);\r
270   allratio->Fit("pol0");\r
271   allratio->Draw();\r
272   behratio->Draw("psame");\r
273 \r
274   TLegend *heleg = new TLegend(0.4,0.75,0.75,0.9);\r
275   heleg->SetTextSize(heleg->GetTextSize()*1.5);\r
276   heleg->AddEntry(allratio,"All electrons","l");\r
277   heleg->AddEntry(behratio,"Bottom electrons","p");\r
278   heleg->Draw();\r
279 \r
280 }\r
281 \r
282 TH1F* GetPtCutHisto(TH1F* input) \r
283 {\r
284   //Given a rate histogram vs pt, return the histogram with yield\r
285   //above a given pTcut\r
286 \r
287   TH1F* result = (TH1F*)input->Clone();\r
288   char name[100];\r
289   sprintf(name,"%s_ptCut",result->GetName());\r
290   result->SetNameTitle(name,name);\r
291   for(Int_t i = 1; i <= result->GetNbinsX(); i++) {\r
292     Double_t val = input->Integral(i,result->GetNbinsX());\r
293     result->SetBinContent(i,val);\r
294     result->SetBinError(i,0.);\r
295   }\r
296 \r
297   return result;\r
298 \r
299 }\r
300 \r