Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / lib / AliHighPtDeDxMc.cxx
CommitLineData
4ebdd20e 1#include "AliHighPtDeDxMc.h"
2
3#ifndef __IOSTREAM__
4#include <iostream>
5#endif
6
7using namespace std;
8
9ClassImp(AliHighPtDeDxMc);
10
11//
12// AliHighPtDeDxMc class
13//
14// This class contains the AliHighPtDeDxMc information
15//
16
17//_________________________________________________________
18AliHighPtDeDxMc::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//_________________________________________________________
44AliHighPtDeDxMc::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//_________________________________________________________
70AliHighPtDeDxMc::~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//_________________________________________________________
91void 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//_________________________________________________________
180Bool_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//_________________________________________________________
193void 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//_________________________________________________________
240TH1D* 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//_________________________________________________________
256void 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//_________________________________________________________
279void 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
318TH1D* 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}