Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / lib / AliHighPtDeDxMc.cxx
1 #include "AliHighPtDeDxMc.h"
2
3 #ifndef __IOSTREAM__
4 #include <iostream>
5 #endif
6
7 using namespace std;
8
9 ClassImp(AliHighPtDeDxMc);
10
11 //
12 // AliHighPtDeDxMc class
13 //
14 // This class contains the AliHighPtDeDxMc information 
15 //
16
17 //_________________________________________________________
18 AliHighPtDeDxMc::AliHighPtDeDxMc():
19   AliHighPtDeDxBase(),
20   fTrackChargeMc(0),
21   fTrackEtaMc(-999),
22   fTrackPtMc(-1),
23   hVtxStatusMc(0x0),
24   hNeventsMc(0x0),
25   hNeventsMcTrig(0x0),
26   hPtMc(0x0),
27   hPtMcNeg(0x0),
28   hPtMcPos(0x0),
29   hPtPiMc(0x0),
30   hPtPiMcNeg(0x0),
31   hPtPiMcPos(0x0),
32   hPtKMc(0x0),
33   hPtKMcNeg(0x0),
34   hPtKMcPos(0x0),
35   hPtPMc(0x0),
36   hPtPMcNeg(0x0),
37   hPtPMcPos(0x0),
38   hMeanPtMc(0x0)
39 {
40   // default constructor - do not use
41 }
42
43 //_________________________________________________________
44 AliHighPtDeDxMc::AliHighPtDeDxMc(const char* name, const char* title):
45   AliHighPtDeDxBase(name, title),
46   fTrackChargeMc(0),
47   fTrackEtaMc(-999),
48   fTrackPtMc(-1),
49   hVtxStatusMc(0x0),
50   hNeventsMc(0x0),
51   hNeventsMcTrig(0x0),
52   hPtMc(0x0),
53   hPtMcNeg(0x0),
54   hPtMcPos(0x0),
55   hPtPiMc(0x0),
56   hPtPiMcNeg(0x0),
57   hPtPiMcPos(0x0),
58   hPtKMc(0x0),
59   hPtKMcNeg(0x0),
60   hPtKMcPos(0x0),
61   hPtPMc(0x0),
62   hPtPMcNeg(0x0),
63   hPtPMcPos(0x0),
64   hMeanPtMc(0x0)
65 {
66   // named constructor
67 }
68
69 //_________________________________________________________
70 AliHighPtDeDxMc::~AliHighPtDeDxMc()
71 {
72   delete hVtxStatusMc;
73   delete hNeventsMc;
74   delete hNeventsMcTrig;
75   delete hPtMc;
76   delete hPtMcNeg;
77   delete hPtMcPos;
78   delete hPtPiMc;
79   delete hPtPiMcNeg;
80   delete hPtPiMcPos;
81   delete hPtKMc;
82   delete hPtKMcNeg;
83   delete hPtKMcPos;
84   delete hPtPMc;
85   delete hPtPMcNeg;
86   delete hPtPMcPos;
87   delete hMeanPtMc;
88 }
89
90 //_________________________________________________________
91 void AliHighPtDeDxMc::Init(Int_t nPtBins, Double_t* ptBins)
92 {
93   //
94   // Create histograms
95   //
96
97   AliHighPtDeDxBase::Init(nPtBins, ptBins);
98   
99   hVtxStatusMc = new TH1D("hVtxStatus", "Vtx status (MC) - No Vtx = -1, Vtx outside cut = 0, Vtx inside = 1",
100                         3, -1.5, 1.5);
101   hVtxStatusMc->Sumw2();
102   hVtxStatusMc->SetDirectory(0);
103
104   hNeventsMc = new TH1D("hNeventsMc", "N events (Mc vtx) - No Vtx = 0, Vtx OK = 1",
105                         2, 0, 2);
106   hNeventsMc->Sumw2();
107   hNeventsMc->SetDirectory(0);
108
109   hNeventsMcTrig = new TH1D("hNeventsMcTrig", "N events (MC vtx+trigger) - No Vtx = 0, Vtx OK = 1",
110                            2, 0, 2);
111   hNeventsMcTrig->Sumw2();
112   hNeventsMcTrig->SetDirectory(0);
113
114   hPtMc = new TH1D("hPtMc", "p_{T} input spectrum (MC); p_{T} [GeV/c]; Counts",
115                    nPtBins, ptBins);
116   hPtMc->Sumw2();
117   hPtMc->SetDirectory(0);
118
119   hPtMcNeg = new TH1D("hPtMcNeg", "p_{T} input spectrum (MC) (q<0); p_{T} [GeV/c]; Counts",
120                    nPtBins, ptBins);
121   hPtMcNeg->Sumw2();
122   hPtMcNeg->SetDirectory(0);
123
124   hPtMcPos = new TH1D("hPtMcPos", "p_{T} input spectrum (MC) (q>0); p_{T} [GeV/c]; Counts",
125                    nPtBins, ptBins);
126   hPtMcPos->Sumw2();
127   hPtMcPos->SetDirectory(0);
128
129   hPtPiMc = new TH1D("hPtPiMc", "p_{T} input spectrum for pi (MC); p_{T} [GeV/c]; Counts",
130                      nPtBins, ptBins);
131   hPtPiMc->Sumw2();
132   hPtPiMc->SetDirectory(0);
133
134   hPtPiMcNeg = new TH1D("hPtPiMcNeg", "p_{T} input spectrum for pi (MC) (q<0); p_{T} [GeV/c]; Counts",
135                         nPtBins, ptBins);
136   hPtPiMcNeg->Sumw2();
137   hPtPiMcNeg->SetDirectory(0);
138   
139   hPtPiMcPos = new TH1D("hPtPiMcPos", "p_{T} input spectrum for pi (MC) (q>0); p_{T} [GeV/c]; Counts",
140                         nPtBins, ptBins);
141   hPtPiMcPos->Sumw2();
142   hPtPiMcPos->SetDirectory(0);
143
144   hPtKMc = new TH1D("hPtKMc", "p_{T} input spectrum for k (MC); p_{T} [GeV/c]; Counts",
145                      nPtBins, ptBins);
146   hPtKMc->Sumw2();
147   hPtKMc->SetDirectory(0);
148
149   hPtKMcNeg = new TH1D("hPtKMcNeg", "p_{T} input spectrum for k (MC) (q<0); p_{T} [GeV/c]; Counts",
150                         nPtBins, ptBins);
151   hPtKMcNeg->Sumw2();
152   hPtKMcNeg->SetDirectory(0);
153   
154   hPtKMcPos = new TH1D("hPtKMcPos", "p_{T} input spectrum for k (MC) (q>0); p_{T} [GeV/c]; Counts",
155                         nPtBins, ptBins);
156   hPtKMcPos->Sumw2();
157   hPtKMcPos->SetDirectory(0);
158   
159   hPtPMc = new TH1D("hPtPMc", "p_{T} input spectrum for p (MC); p_{T} [GeV/c]; Counts",
160                      nPtBins, ptBins);
161   hPtPMc->Sumw2();
162   hPtPMc->SetDirectory(0);
163
164   hPtPMcNeg = new TH1D("hPtPMcNeg", "p_{T} input spectrum for p (MC) (q<0); p_{T} [GeV/c]; Counts",
165                         nPtBins, ptBins);
166   hPtPMcNeg->Sumw2();
167   hPtPMcNeg->SetDirectory(0);
168   
169   hPtPMcPos = new TH1D("hPtPMcPos", "p_{T} input spectrum for p (MC) (q>0); p_{T} [GeV/c]; Counts",
170                         nPtBins, ptBins);
171   hPtPMcPos->Sumw2();
172   hPtPMcPos->SetDirectory(0);
173
174   hMeanPtMc = new TProfile("hMeanPtMc", "mean p_{T}; p_{T} [GeV/c]; mean p_{T}",
175                            nPtBins, ptBins);
176   hMeanPtMc->SetDirectory(0);
177 }
178
179 //_________________________________________________________
180 Bool_t AliHighPtDeDxMc::TrackAcceptedMc()
181 {
182   if(fUseEtaCut && (fTrackEtaMc<fEtaLow || fTrackEtaMc>fEtaHigh))
183     return kFALSE;
184   
185   // only accept hadrons = pi(1), K(2), p(3), other (Sigma+ etc. = 6)
186   if(fTrackPidMc < 1 || fTrackPidMc == 4 || fTrackPidMc == 5 || fTrackPidMc > 6)
187     return kFALSE;
188   
189   return kTRUE;
190 }
191
192 //_________________________________________________________
193 void AliHighPtDeDxMc::FillEventInfo() 
194 {
195   //
196   // We require that the MC vtx is withion our vertex range
197   // (fEventVtxStatusMc==1). In this way we have to accept both
198   // fEventVtxStatusMc==0 and fEventVtxStatusMc==1 to take into account
199   // migration effects. 
200   //
201   // The histogram associated with this is hNeventsMc (together with hNevents
202   // in the base class)
203   //
204   // To correct down the "no vtx" events in the data there are two possibilities
205   //
206   // 1) We can use the data fraction between accepted vtx and all vtx. Note
207   // that since the vtx rec efficiency is not constant vs the vtx position we
208   // make a small mistake here, but on the other hand we might be less
209   // sensitive to the MC.
210   //
211   // 2) We can use the MC prediction for estimating the fraction. In this way
212   // we might be have no problems with the slightly varying vtx efficiency,
213   // but we might be more sensitive to the MC.
214   //
215   // The histogram associated with this method is hNeventsMcTrig (together with hNevents
216   // in the base class)
217   //
218   // TO DO: We need to test which is better, e.g., by using PHOJET on PYTHIA
219   // and vice versa.
220   //
221   AliHighPtDeDxBase::FillEventInfo();
222  
223   hVtxStatusMc->Fill(fEventVtxStatus);
224
225   if(fEventVtxStatusMc==1) {
226     if(fEventVtxStatus==-1) { // no vtx class
227       hNeventsMc->Fill(0.5);      
228       if(fEventTrigger)
229         hNeventsMcTrig->Fill(0.5);
230     } else {                  // vtx class
231       hNeventsMc->Fill(1.5);
232       if(fEventTrigger)
233         hNeventsMcTrig->Fill(1.5);
234     }
235   }
236 }
237
238
239 //_________________________________________________________
240 TH1D* AliHighPtDeDxMc::GetPtSpectrum() 
241 {
242   TH1D* histSpectrum = dynamic_cast<TH1D*>(hPtMc->Clone("hPtSpectrum"));
243
244   const Double_t nEvents = hNeventsMc->GetBinContent(1) 
245     + hNeventsMc->GetBinContent(2);
246   const Double_t etaRange = fEtaHigh - fEtaLow;
247   
248   histSpectrum->Scale(1.0/etaRange);
249   NormalizePt(histSpectrum);
250   histSpectrum->Scale(1.0/nEvents);
251
252   return histSpectrum;
253 }
254
255 //_________________________________________________________
256 void AliHighPtDeDxMc::NormalizePt(TH1* hist)
257 {
258   const Int_t n = hist->GetNbinsX();
259
260   for(Int_t bin = 1; bin <= n; bin++) {
261     
262     const Float_t width = hist->GetXaxis()->GetBinWidth(bin);
263     hist->SetBinError(bin, hist->GetBinError(bin)/width);
264     hist->SetBinContent(bin, hist->GetBinContent(bin)/width);
265   }
266 }
267
268
269 // //_________________________________________________________
270 // void AliHighPtDeDxMc::FillTrackInfo(Float_t weight) 
271 // {
272   
273 //   AliHighPtDeDxBase::FillTrackInfo(weight);
274   
275 //   hPt->Fill(fTrackPt, weight);
276 // }
277
278 //_________________________________________________________
279 void AliHighPtDeDxMc::FillTrackInfoMc(Float_t weight) 
280 {
281   
282   hPtMc->Fill(fTrackPtMc, weight);
283   if(fTrackChargeMc<0)
284     hPtMcNeg->Fill(fTrackPtMc, weight);
285   else
286     hPtMcPos->Fill(fTrackPtMc, weight);
287
288   hMeanPtMc->Fill(fTrackPtMc, fTrackPtMc);
289
290   switch (fTrackPidMc) {
291     
292   case 1: // pion
293     hPtPiMc->Fill(fTrackPtMc, weight);
294     if(fTrackChargeMc<0)
295       hPtPiMcNeg->Fill(fTrackPtMc, weight);
296     else
297       hPtPiMcPos->Fill(fTrackPtMc, weight);
298     break;
299   case 2: // kaon
300     hPtKMc->Fill(fTrackPtMc, weight);
301     if(fTrackChargeMc<0)
302       hPtKMcNeg->Fill(fTrackPtMc, weight);
303     else
304       hPtKMcPos->Fill(fTrackPtMc, weight);
305     break;
306   case 3: // proton
307     hPtPMc->Fill(fTrackPtMc, weight);
308     if(fTrackChargeMc<0)
309       hPtPMcNeg->Fill(fTrackPtMc, weight);
310     else
311       hPtPMcPos->Fill(fTrackPtMc, weight);
312     break;
313   default:
314     break;
315   }
316 }
317
318 TH1D* AliHighPtDeDxMc::GetHistPtMc(Int_t pid, Int_t charge)
319 {
320   switch (pid) {
321         
322   case 0:
323     if(charge==0)
324       return hPtMc;
325     else if(charge<0)
326       return hPtMcNeg;
327     else
328       return hPtMcPos;
329     break;
330   case 1:
331     if(charge==0)
332       return hPtPiMc;
333     else if(charge<0)
334       return hPtPiMcNeg;
335     else
336       return hPtPiMcPos;
337     break;
338   case 2:
339     if(charge==0)
340       return hPtKMc;
341     else if(charge<0)
342       return hPtKMcNeg;
343     else
344       return hPtKMcPos;
345     break;
346   case 3:
347     if(charge==0)
348       return hPtPMc;
349     else if(charge<0)
350       return hPtPMcNeg;
351     else
352       return hPtPMcPos;
353     break;
354   default:
355     cout << "PID: " << pid << " not found" << endl;
356     break;
357   }
358   return 0;
359 }