4bb600af2eab4f2add574d131e2748cd4614f69b
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / PlotOutputMCCheck.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TCanvas.h>
3 #include <TH1.h>
4 #include <TH2.h>
5 #include <TStyle.h>
6 #include <TFile.h>
7 #endif
8
9 /* $Id$ */ 
10
11 // Macro to plot the output of AliAnalysisTaskCheckHFMCProd
12 // Author: F. Prino, prino@to.infn.it
13
14
15 void PlotOutputMCCheck(){
16   TFile *fil=new TFile("AnalysisResults.root");
17   TDirectoryFile* df=(TDirectoryFile*)fil->Get("HFMCCheck");
18   TList* l=(TList*)df->Get("clistHFMCCheck");
19   l->ls();
20
21   TH1F* hNEvents=(TH1F*)l->FindObject("hNEvents");
22   Int_t nAnalEv=hNEvents->GetBinContent(1);
23   printf("Number of events= %d\n",nAnalEv);
24
25
26   TCanvas* cv=new TCanvas("cv","Vertex");
27   cv->Divide(3,3);
28   cv->cd(1);
29   TH1F* hSPD3DvX=(TH1F*)l->FindObject("hSPD3DvX");
30   hSPD3DvX->Draw();
31   cv->cd(2);
32   TH1F* hSPD3DvY=(TH1F*)l->FindObject("hSPD3DvY");
33   hSPD3DvY->Draw();
34   cv->cd(3);
35   TH1F* hSPD3DvZ=(TH1F*)l->FindObject("hSPD3DvZ");
36   hSPD3DvZ->Draw();
37   cv->cd(4);
38   TH1F* hSPDZvX=(TH1F*)l->FindObject("hSPDZvX");
39   hSPDZvX->Draw();
40   cv->cd(5);
41   TH1F* hSPDZvY=(TH1F*)l->FindObject("hSPDZvY");
42   hSPDZvY->Draw();
43   cv->cd(6);
44   TH1F* hSPDZvZ=(TH1F*)l->FindObject("hSPDZvZ");
45   hSPDZvZ->Draw();
46   cv->cd(7);
47   TH1F* hTRKvX=(TH1F*)l->FindObject("hTRKvX");
48   hTRKvX->Draw();
49   cv->cd(8);
50   TH1F* hTRKvY=(TH1F*)l->FindObject("hTRKvY");
51   hTRKvY->Draw();
52   cv->cd(9);
53   TH1F* hTRKvZ=(TH1F*)l->FindObject("hTRKvZ");
54   hTRKvZ->Draw();
55
56   TCanvas* c1=new TCanvas("c1","Multipl");
57   c1->Divide(2,2);
58   c1->cd(1);
59   gPad->SetLogy();
60   TH1F* hPhysPrim=(TH1F*)l->FindObject("hPhysPrim");
61   hPhysPrim->Draw();
62   c1->cd(2);
63   gPad->SetLogy();
64   TH1F* hTracklets=(TH1F*)l->FindObject("hTracklets");
65   hTracklets->Draw();
66   c1->cd(3);
67   gPad->SetLogy();
68   TH1F* hTracks=(TH1F*)l->FindObject("hTracks");
69   hTracks->Draw();
70   c1->cd(4);
71   gPad->SetLogy();
72   TH1F* hSelTracks=(TH1F*)l->FindObject("hSelTracks");
73   hSelTracks->Draw();
74
75   TH1F* hncharmed=(TH1F*)l->FindObject("hncharmed");
76   TCanvas* cn=new TCanvas("cn","ncharm");
77   hncharmed->Draw("box");
78   hncharmed->GetXaxis()->SetTitle("dNch/dy");
79   hncharmed->GetYaxis()->SetTitle("N Charm hadrons in golden channels");
80   cn->Update();
81
82   TH1F* hnbvsnc=(TH1F*)l->FindObject("hnbvsnc");
83   TCanvas* cnhf=new TCanvas("cnhf","nb/c");
84   hnbvsnc->Draw("colztext");
85   hnbvsnc->GetXaxis()->SetTitle("Nc");
86   hnbvsnc->GetYaxis()->SetTitle("Nb");
87   cnhf->Update();
88   TH2F* hyptD0all=(TH2F*)l->FindObject("hyptD0AllDecay");
89   TH2F* hyptD0promptall=(TH2F*)l->FindObject("hyptD0promptAllDecay");
90   TH2F* hyptD0feeddownall=(TH2F*)l->FindObject("hyptD0feeddownAllDecay");
91   TH2F* hyptD0prompt=(TH2F*)l->FindObject("hyptD0prompt");
92   TH2F* hyptD0feeddown=(TH2F*)l->FindObject("hyptD0feeddown");
93   TH2F* hyptD02=(TH2F*)l->FindObject("hyptD02");
94   TH2F* hyptD04=(TH2F*)l->FindObject("hyptD04");
95   TH1D* hptD0all=hyptD0all->ProjectionX();
96   TH1D* hptD0promptall=hyptD0promptall->ProjectionX();
97   TH1D* hptD0feeddownall=hyptD0feeddownall->ProjectionX();
98   TH2F* hyptDplusall=(TH2F*)l->FindObject("hyptDplusAllDecay");
99   TH2F* hyptDpluspromptall=(TH2F*)l->FindObject("hyptDpluspromptAllDecay");
100   TH2F* hyptDplusfeeddownall=(TH2F*)l->FindObject("hyptDplusfeeddownAllDecay");
101   TH2F* hyptDplusprompt=(TH2F*)l->FindObject("hyptDplusprompt");
102   TH2F* hyptDplusfeeddown=(TH2F*)l->FindObject("hyptDplusfeeddown");
103   TH2F* hyptDplusnonreson=(TH2F*)l->FindObject("hyptDplusnonreson");
104   TH2F* hyptDplusreson=(TH2F*)l->FindObject("hyptDplusreson");
105   TH2F* hyptDsall=(TH2F*)l->FindObject("hyptDsAllDecay");
106   TH2F* hyptDsprompt=(TH2F*)l->FindObject("hyptDsprompt");
107   TH2F* hyptDsfeeddown=(TH2F*)l->FindObject("hyptDsfeedown");
108   TH2F* hyptDsphi=(TH2F*)l->FindObject("hyptDsphi");
109   TH2F* hyptDsK0st=(TH2F*)l->FindObject("hyptDsk0st");
110   TH2F* hyptDstarall=(TH2F*)l->FindObject("hyptDstarAllDecay");
111   TH2F* hyptDstarprompt=(TH2F*)l->FindObject("hyptDstarprompt");
112   TH2F* hyptDstarfeedown=(TH2F*)l->FindObject("hyptDstarfeedown");
113   TH2F* hyptLcprompt=(TH2F*)l->FindObject("hyptLcprompt");
114   TH2F* hyptLcfeedown=(TH2F*)l->FindObject("hyptLcfeedown");
115   TH2F* hyptLcall=(TH2F*)l->FindObject("hyptLcAllDecay");
116
117   TH2F* hyptB0all=(TH2F*)l->FindObject("hyptB0AllDecay");
118   TH2F* hyptBplusall=(TH2F*)l->FindObject("hyptBplusAllDecay");
119   TH2F* hyptBsall=(TH2F*)l->FindObject("hyptBsAllDecay");
120   TH2F* hyptBstarall=(TH2F*)l->FindObject("hyptBstarAllDecay");
121   TH2F* hyptLball=(TH2F*)l->FindObject("hyptLbAllDecay");
122
123
124
125   TH1F* hD0fonll7=HistoFONLL7TeV();
126   hD0fonll7->Scale(hptD0all->GetMaximum()/hD0fonll7->GetMaximum());  
127   hD0fonll7->SetLineColor(kGreen+1);
128   hD0fonll7->SetLineWidth(2);
129
130   TH1F* hD0fonll2=HistoFONLL2_76TeV();
131   hD0fonll2->Scale(hptD0all->GetMaximum()/hD0fonll2->GetMaximum());  
132   hD0fonll2->SetLineColor(kBlue+1);
133   hD0fonll2->SetLineWidth(2);
134
135   TH1F* hptD0pythia=HistoPYTHIA7();
136   hptD0pythia->Scale(hptD0all->GetMaximum()/hptD0pythia->GetMaximum());
137   hptD0pythia->SetLineColor(kRed+1);
138   hptD0pythia->SetLineWidth(2);
139
140   hptD0all->SetLineWidth(2);
141   hptD0all->SetMarkerStyle(20);
142   hptD0all->SetTitle("");
143   hptD0all->GetYaxis()->SetTitle("D^{0} dN/dp_{T} (a.u.)");
144   hptD0all->GetXaxis()->SetTitle("p_{T} (GeV/c)");
145   hptD0all->SetStats(0);  
146   hptD0all->GetYaxis()->SetTitleOffset(1.2);
147   hptD0all->GetXaxis()->SetTitleOffset(1.2);
148
149   TCanvas* cd0a=new TCanvas("cd0a","D0 spectra",700,700);
150   cd0a->SetLogy();
151   cd0a->SetLeftMargin(0.13);
152   cd0a->SetRightMargin(0.07);
153   hptD0all->Draw();
154   hD0fonll7->Draw("lsame");
155   hD0fonll2->Draw("lsame");
156   hptD0pythia->Draw("lsame");
157   TLegend* leg=new TLegend(0.45,0.6,0.89,0.85);
158   leg->SetFillStyle(0);
159   leg->SetBorderSize(0);
160   leg->AddEntry(hptD0all,"MC production","LP");
161   leg->AddEntry(hD0fonll7,"FONLL, #sqrt{s}=7 TeV","L");
162   leg->AddEntry(hD0fonll2,"FONLL, #sqrt{s}=2.76 TeV","L");
163   leg->AddEntry(hptD0pythia,"PYTHIA Perugia0, #sqrt{s}=7 TeV","L");
164   leg->Draw();
165
166   // Prompt and Feeddown
167   TH1F* hOriginPrompt=(TH1F*)l->FindObject("hOriginPrompt");
168   TH1F* hOriginFeeddown=(TH1F*)l->FindObject("hOriginFeeddown");
169
170   hptD0promptall->SetLineColor(4);
171   hptD0promptall->SetMarkerColor(4);
172   hptD0promptall->SetMarkerStyle(26);
173   hptD0feeddownall->SetLineColor(2);
174   hptD0feeddownall->SetMarkerColor(2);
175   hptD0feeddownall->SetMarkerStyle(23)
176 ;
177   TCanvas* cprf1=new TCanvas("cprf1","Prompt/Feeddown",700,700);
178   gPad->SetLogy();
179   gPad->SetLeftMargin(0.13);
180   gPad->SetRightMargin(0.07);
181   hptD0all->Draw();
182   hptD0promptall->Draw("same");
183   hptD0feeddownall->Draw("same");
184   TLegend* leg2=new TLegend(0.4,0.5,0.89,0.85);
185   leg2->SetFillStyle(0);
186   leg2->SetBorderSize(0);
187   leg2->AddEntry(hptD0all,"All D0","LP");
188   leg2->AddEntry("",Form("Entries=%.0f  <pt>=%.2f GeV/c",hptD0all->Integral(),hptD0all->GetMean()),"");
189   leg2->AddEntry(hptD0promptall,"Prompt D0","LP");
190   leg2->AddEntry("",Form("Entries=%.0f  <pt>=%.2f GeV/c",hptD0promptall->Integral(),hptD0promptall->GetMean()),"");
191   leg2->AddEntry(hptD0feeddownall,"Feeddown D0","LP");
192   leg2->AddEntry("",Form("Entries=%.0f  <pt>=%.2f GeV/c",hptD0feeddownall->Integral(),hptD0feeddownall->GetMean()),"");
193   leg2->Draw();
194   hOriginPrompt->GetXaxis()->SetTitle("Distance of prompt D meson origin to vertex (cm)");
195   hOriginFeeddown->GetXaxis()->SetTitle("Distance of feed-down D meson origin to vertex (cm)");
196
197   TCanvas* cprf2=new TCanvas("cprf2","Origin",1000,600);
198   cprf2->Divide(2,1);
199   cprf2->cd(1);
200   hOriginPrompt->Draw();
201   cprf2->cd(2);
202   hOriginFeeddown->Draw();
203
204   // Hadrons
205   Double_t nev=hNEvents->GetBinContent(1);
206   Double_t nD0=hyptD0all->GetEntries();
207   Double_t nDp=hyptDplusall->GetEntries();
208   Double_t nDs=hyptDsall->GetEntries();
209   Double_t nDst=hyptDstarall->GetEntries();
210   Double_t nLc=hyptLcall->GetEntries();
211   Double_t nB0=hyptB0all->GetEntries();
212   Double_t nBp=hyptBplusall->GetEntries();
213   Double_t nBs=hyptBsall->GetEntries();
214   Double_t nBst=hyptBstarall->GetEntries();
215   Double_t nLb=hyptLball->GetEntries();
216
217
218   Double_t nccbar=0;
219   Double_t nbbbar=0;
220   for(Int_t iBinx=1; iBinx<=hnbvsnc->GetNbinsX(); iBinx++){
221     for(Int_t iBiny=1; iBiny<=hnbvsnc->GetNbinsY(); iBiny++){
222       Double_t bincentx=hnbvsnc->GetXaxis()->GetBinCenter(iBinx);
223       Double_t bincenty=hnbvsnc->GetYaxis()->GetBinCenter(iBiny);
224       Double_t bincont=hnbvsnc->GetBinContent(iBinx,iBiny);
225       nccbar+=(bincentx*bincont);
226       nbbbar+=(bincenty*bincont);
227     }
228   }
229
230
231   printf("Events =%f\n",nev);
232   printf("c+cbar =%f\n",nccbar);
233   printf("D0 =%f\n",nD0);
234   printf("D+ =%f\n",nDp);
235   printf("D*+ =%f\n",nDst);
236   printf("Ds =%f\n",nDs);
237   printf("Lc =%f\n",nLc);
238   printf("----------\n");
239   printf("b+bbar =%f\n",nbbbar);
240   printf("B0 =%f\n",nB0);
241   printf("B+ =%f\n",nBp);
242   printf("B*0 =%f\n",nBst);
243   printf("Bs =%f\n",nBs);
244   printf("Lb =%f\n",nLb);
245
246   if(nccbar==0) nccbar=nD0+nDp+nDs+nLc;
247   if(nbbbar==0) nbbbar=nB0+nBp+nBs+nLb;
248
249
250   TH1F* hCharmHad=new TH1F("hCharmHad","",5,-0.5,4.5);
251   hCharmHad->GetXaxis()->SetBinLabel(1,"D0");
252   hCharmHad->SetBinContent(1,nD0/nccbar);
253   hCharmHad->GetXaxis()->SetBinLabel(2,"D+");
254   hCharmHad->SetBinContent(2,nDp/nccbar);
255   hCharmHad->GetXaxis()->SetBinLabel(3,"Ds");
256   hCharmHad->SetBinContent(3,nDs/nccbar);
257   hCharmHad->GetXaxis()->SetBinLabel(4,"Lc");
258   hCharmHad->SetBinContent(4,nLc/nccbar);
259   hCharmHad->GetXaxis()->SetBinLabel(5,"D*+");
260   hCharmHad->SetBinContent(5,nDst/nccbar);
261   hCharmHad->SetMinimum(0);
262   hCharmHad->SetStats(0);
263   hCharmHad->GetYaxis()->SetTitle("N(species)/N(c+#bar{c})");
264   hCharmHad->GetYaxis()->SetTitleOffset(1.3);
265
266   TH1F* hBeautyHad=new TH1F("hBeautyHad","",5,-0.5,4.5);
267   hBeautyHad->GetXaxis()->SetBinLabel(1,"B0");
268   hBeautyHad->SetBinContent(1,nB0/nbbbar);
269   hBeautyHad->GetXaxis()->SetBinLabel(2,"B+");
270   hBeautyHad->SetBinContent(2,nBp/nbbbar);
271   hBeautyHad->GetXaxis()->SetBinLabel(3,"Bs");
272   hBeautyHad->SetBinContent(3,nBs/nbbbar);
273   hBeautyHad->GetXaxis()->SetBinLabel(4,"Lb");
274   hBeautyHad->SetBinContent(4,nLb/nbbbar);
275   hBeautyHad->GetXaxis()->SetBinLabel(5,"B*0");
276   hBeautyHad->SetBinContent(5,nBst/nbbbar);
277   hBeautyHad->SetMinimum(0);
278   hBeautyHad->SetStats(0);
279   hBeautyHad->GetYaxis()->SetTitle("N(species)/N(b+#bar{b})");
280   hBeautyHad->GetYaxis()->SetTitleOffset(1.3);
281
282
283
284   TCanvas* chad=new TCanvas("chad","Hadrons",1200,600);
285   chad->Divide(2,1);
286   chad->cd(1);
287   hCharmHad->Draw();
288   chad->cd(2);
289   hBeautyHad->Draw();
290   
291
292
293   TCanvas* cd0=new TCanvas("cd0","D0");
294   cd0->Divide(2,2);
295   cd0->cd(1);
296   hyptD0prompt->Draw("colz");
297   cd0->cd(2);
298   hyptD0feeddown->Draw("colz");
299   cd0->cd(3);
300   hyptD02->Draw("colz");
301   cd0->cd(4);
302   hyptD04->Draw("colz");
303
304
305   TCanvas* cdplus=new TCanvas("cdplus","Dplus");
306   cdplus->Divide(2,2);
307   cdplus->cd(1);
308   hyptDplusprompt->Draw("colz");
309   cdplus->cd(2);
310   hyptDplusfeeddown->Draw("colz");
311   cdplus->cd(3);
312   hyptDplusnonreson->Draw("colz");
313   cdplus->cd(4);
314   hyptDplusreson->Draw("colz");
315
316    
317   TCanvas* cds=new TCanvas("cds","Ds");
318   cds->Divide(2,2);
319   cds->cd(1);
320   hyptDsprompt->Draw("colz");
321   cds->cd(2);
322   hyptDsfeeddown->Draw("colz");
323   cds->cd(3);
324   hyptDsphi->Draw("colz");
325   cds->cd(4);
326   hyptDsK0st->Draw("colz");
327
328
329   TCanvas* cdstlc=new TCanvas("cdstls","Dstar LambdaC");
330   cdstlc->Divide(2,2);
331   cdstlc->cd(1);
332   hyptDstarprompt->Draw("colz");
333   cdstlc->cd(2);
334   hyptDstarfeedown->Draw("colz");
335   cdstlc->cd(3);
336   hyptLcprompt->Draw("colz");
337   cdstlc->cd(4);
338   hyptLcfeedown->Draw("colz");
339
340 }
341
342
343
344 TH1F* HistoFONLL7TeV(){
345
346   TH1F* hFONLL7=new TH1F("hFONLLD07TeV","",61,0.,30.5);
347   Float_t val[61]={
348     1390542.31,3512269.33,4572233.65,4116353.65,3104057.40,
349     2185147.21,1507632.40,1038687.03,721889.43,509311.55,
350     365094.01,265684.43,196609.08,147415.64,112019.94,
351     86170.85,66997.46,52651.54,41787.55,33458.64,
352     27012.62,21981.48,18020.69,14873.73,12354.91,
353     10324.90,8677.32,7331.56,6225.67,5311.74,
354     4552.30,3918.10,3385.73,2936.79,2556.52,
355     2233.25,1957.23,1720.64,1517.12,1341.44,
356     1189.32,1057.17,942.05,841.47,753.33,
357     675.89,607.68,547.45,494.14,446.87,
358     404.83,367.39,333.96,304.08,277.30,
359     253.25,231.62,212.13,194.55,178.66,
360     164.27};
361   for(Int_t i=0; i<61; i++) hFONLL7->SetBinContent(i+1,val[i]);
362   return hFONLL7;
363 }
364
365 TH1F* HistoFONLL2_76TeV(){
366   TH1F* hFONLL2=new TH1F("hFONLLD02_76TeV","",61,0.,30.5);
367   Float_t val[61]={
368     1154043.73,2596175.89,2995937.57,2442988.08,1701598.07,
369     1122452.81,732896.42,482314.10,322062.75,218493.05,
370     151653.77,107053.80,76999.08,56239.64,41687.89,
371     31322.25,23818.86,18326.64,14252.10,11192.03,
372     8869.29,7089.11,5711.50,4635.50,3788.31,
373     3116.19,2578.84,2146.51,1796.18,1510.67,
374     1276.72,1083.93,924.20,791.18,679.88,
375     586.36,507.49,440.73,383.99,335.54,
376     294.07,258.43,227.68,201.11,178.07,
377     158.04,140.58,125.32,111.95,100.20,
378     89.86,80.73,72.66,65.51,59.16,
379     53.51,48.48,43.98,39.96,36.36,
380     33.12};
381   for(Int_t i=0; i<61; i++) hFONLL2->SetBinContent(i+1,val[i]);
382   return hFONLL2;
383 }
384
385 TH1F* HistoPYTHIA7(){
386   TH1F* hPYTHIA7=new TH1F("hPYTHIAD07TeV","",40,0.,20.);
387   Float_t val[40]={
388     826307.00,1753264.00,1877464.00,1634664.00,1278586.00,
389     959137.00,713389.00,535745.00,410250.00,316284.00,
390     249107.00,198235.00,158832.00,129936.00,106957.00,
391     89098.00,73690.00,62752.00,53247.00,45004.00,
392     38475.00,32691.00,28510.00,24516.00,21204.00,
393     18276.00,15890.00,13702.00,12326.00,10612.00,
394     9184.00,8028.00,7194.00,6384.00,5767.00,
395     5102.00,4505.00,3939.00,3578.00,3288.00
396   };
397   for(Int_t i=0; i<40; i++) hPYTHIA7->SetBinContent(i+1,val[i]);
398   return hPYTHIA7;
399 }