Update PR task: drathee
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidRatioEffCont.cxx
CommitLineData
0a28d543 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: ALICE Offline. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17//=========================================================================//
18// AliEbyE Analysis for Particle Ratio Fluctuation //
19// Deepika Rathee | Satyajit Jena //
20// drathee@cern.ch | sjena@cern.ch //
21// Date: Wed Jul 9 18:38:30 CEST 2014 //
22// New approch to find particle ratio to reduce memory //
23// (Test Only) //
24//=========================================================================//
25
26#include "TMath.h"
27#include "TAxis.h"
28
29#include "AliESDEvent.h"
30#include "AliStack.h"
31#include "AliMCEvent.h"
32#include "AliESDtrackCuts.h"
33#include "AliAODEvent.h"
34#include "AliAODMCParticle.h"
35
36#include "AliEbyEPidRatioEffCont.h"
37
38using namespace std;
39
40
41ClassImp(AliEbyEPidRatioEffCont)
42
43//________________________________________________________________________
44AliEbyEPidRatioEffCont::AliEbyEPidRatioEffCont() :
45 AliEbyEPidRatioBase("EffCont", "EffCont"),
46 fLabelsRec(NULL),
f7ea34d2 47 fHnEffMc(NULL),
48 fHnContMc(NULL),
49 fHnEffRec(NULL),
50 fHnContRec(NULL) {
0a28d543 51 // Constructor
52
53 AliLog::SetClassDebugLevel("AliEbyEPidRatioEffCont",10);
54}
55
56//________________________________________________________________________
57AliEbyEPidRatioEffCont::~AliEbyEPidRatioEffCont() {
58 // Destructor
59
60 for (Int_t ii = 0; ii < 2; ++ii) {
61 for (Int_t kk = 0; kk < 4; ++kk)
62 if (fLabelsRec[ii][kk]) delete[] fLabelsRec[ii][kk];
63 if (fLabelsRec[ii]) delete[] fLabelsRec[ii];
64 }
65 if (fLabelsRec) delete[] fLabelsRec;
66
67
68}
69
70//________________________________________________________________________
71void AliEbyEPidRatioEffCont::Process() {
72 // -- Process event
73
74 // -- Setup (clean, create and fill) MC labels
75 FillMCLabels();
76
77 // -- Fill MC histograms for efficiency studies
78 FillMCEffHist();
79
80 return;
81}
82
83//________________________________________________________________________
84void AliEbyEPidRatioEffCont::Init() {
85 // -- Init eventwise
86
87 fLabelsRec = new Int_t**[2];
88 for (Int_t ii = 0 ; ii < 2; ++ii) {
89 fLabelsRec[ii] = new Int_t*[4];
90 for (Int_t kk = 0 ; kk < 4; ++kk)
91 fLabelsRec[ii][kk] = NULL;
92 }
93}
94
95//________________________________________________________________________
96void AliEbyEPidRatioEffCont::CreateHistograms() {
97 // Copied from NetParticle class
f7ea34d2 98 Int_t binHnEff[10] = { AliEbyEPidRatioHelper::fgkfHistNBinsCent, 4,
99 AliEbyEPidRatioHelper::fgkfHistNBinsSign, 2, 2 , 2 ,
100 AliEbyEPidRatioHelper::fgkfHistNBinsEta,
101 AliEbyEPidRatioHelper::fgkfHistNBinsRap,
102 AliEbyEPidRatioHelper::fgkfHistNBinsPhi,
103 AliEbyEPidRatioHelper::fgkfHistNBinsPt};
104
105 Double_t minHnEff[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], -0.5,
106 AliEbyEPidRatioHelper::fgkfHistRangeSign[0], -0.5, -0.5, -0.5,
107 AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
108 AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
109 AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
110 AliEbyEPidRatioHelper::fgkfHistRangePt[0]};
111
112
113 Double_t maxHnEff[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], 3.5,
114 AliEbyEPidRatioHelper::fgkfHistRangeSign[1], 1.5, 1.5, 1.5,
115 AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
116 AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
117 AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
118 AliEbyEPidRatioHelper::fgkfHistRangePt[1]};
119
120 fHnEffMc = new THnSparseF("hnEffMc", "cent:pid:SignMC:findable:recStatus:pidStatus:etaMC:yMC:phiMC:ptMC", 10, binHnEff, minHnEff, maxHnEff);
121 fHnEffRec = new THnSparseF("hnEffRec", "cent:pid:SignMC:findable:recStatus:pidStatus:etaRec:yRec:phiRec:ptRec", 10, binHnEff, minHnEff, maxHnEff);
122 //fHnEffRecMc = new THnSparseF("hnEffRecMMc", "cent:pid:SignMC:findable:recStatus:pidStatus:deltaEta:deltaY:deltaPhi:deltaPt" 10, binHnEff, minHnEff, maxHnEff);
123
124 fHnEffMc->Sumw2();
125 fHnEffRec->Sumw2();
126
127
128
129 fHnEffMc->GetAxis(0)->SetTitle("centrality"); // 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
130 fHnEffMc->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}"); // 0 | 1 | 2 | 3
131 fHnEffMc->GetAxis(2)->SetTitle("sign"); // -1 | 0 | +1
132 fHnEffMc->GetAxis(3)->SetTitle("findable"); // 0 not findable | 1 findable
133 fHnEffMc->GetAxis(4)->SetTitle("recStatus"); // 0 not reconstructed | 1 reconstructed
134 fHnEffMc->GetAxis(5)->SetTitle("recPid"); // 0 not accepted | 1 accepted
135 fHnEffMc->GetAxis(6)->SetTitle("#eta_{MC}"); // eta [-0.9, 0.9]
136 fHnEffMc->GetAxis(7)->SetTitle("#it{y}_{MC}"); // rapidity [-0.5, 0.5]
137 fHnEffMc->GetAxis(8)->SetTitle("#varphi_{MC} (rad)"); // phi [ 0. , 2Pi]
138 fHnEffMc->GetAxis(9)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})"); // pT [ 0.2, 2.3]
139
140
141 fHnEffRec->GetAxis(0)->SetTitle("centrality"); // 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
142 fHnEffRec->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}"); // 0 | 1 | 2 | 3
143 fHnEffRec->GetAxis(2)->SetTitle("sign"); // -1 | 0 | +1
144 fHnEffRec->GetAxis(3)->SetTitle("findable"); // 0 not findable | 1 findable
145 fHnEffRec->GetAxis(4)->SetTitle("recStatus"); // 0 not reconstructed | 1 reconstructed
146 fHnEffRec->GetAxis(5)->SetTitle("recPid"); // 0 not accepted | 1 accepted
147 fHnEffRec->GetAxis(6)->SetTitle("#eta_{Rec}"); // eta [-0.9, 0.9]
148 fHnEffRec->GetAxis(7)->SetTitle("#it{y}_{Rec}"); // rapidity [-0.5, 0.5]
149 fHnEffRec->GetAxis(8)->SetTitle("#varphi_{Rec} (rad)"); // phi [ 0. , 2Pi]
150 fHnEffRec->GetAxis(9)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 0.2, 2.3]
151
0a28d543 152
f7ea34d2 153 // fHnEff->GetAxis(0)->SetTitle("centrality"); // 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
154 // fHnEff->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}"); // 0 | 1 | 2 | 3
155 // fHnEff->GetAxis(2)->SetTitle("sign"); // -1 | 0 | +1
156 // fHnEff->GetAxis(3)->SetTitle("findable"); // 0 not findable | 1 findable
157 // fHnEff->GetAxis(4)->SetTitle("recStatus"); // 0 not reconstructed | 1 reconstructed
158 // fHnEff->GetAxis(5)->SetTitle("recPid"); // 0 not accepted | 1 accepted
159 // fHnEff->GetAxis(6)->SetTitle("#eta_{MC}-#eta_{Rec}"); // eta [-0.9, 0.9]
160 // fHnEff->GetAxis(7)->SetTitle("#it{y}_{MC}-#it{y}_{Rec}"); // rapidity [-0.5, 0.5]
161 // fHnEff->GetAxis(8)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); // phi [ -2Pi , 2Pi]
162 // fHnEff->GetAxis(9)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ -2.3, 2.3]
163 // fHnEff->GetAxis(10)->SetTitle("sign_{MC}-sign_{Rec}"); // -2 | 0 | +2
0a28d543 164
0a28d543 165
f7ea34d2 166
6ce4ad92 167
f7ea34d2 168 fHelper->BinLogAxis(fHnEffMc, 9);
f7ea34d2 169 fHelper->BinLogAxis(fHnEffRec, 9);
0a28d543 170
171 /*
172 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Copied from NetParticle
173
174 creation -> findable -> reconstructed -> pid_TPC+TOF
175 (1) (2) (3) (4)
176 || findable | recStatus | recPid
177 1) all primary probeParticles_MC || - - -
178 2) all findable primary probeParticles_MC || x - -
179 3) all reconstructed primary probeParticles_MC || x x -
180 4) all reconstructed primary probeParticles_MC & recPid_TPC+TOF || x x x
181
182 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
183 */
184
185 // ------------------------------------------------------------------
186 // -- Create THnSparse - Cont
187 // ------------------------------------------------------------------
188
f7ea34d2 189 Int_t binHnCont[8] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, 4,
190 AliEbyEPidRatioHelper::fgkfHistNBinsSign, 8,
191 AliEbyEPidRatioHelper::fgkfHistNBinsEta,
192 AliEbyEPidRatioHelper::fgkfHistNBinsRap,
193 AliEbyEPidRatioHelper::fgkfHistNBinsPhi,
194 AliEbyEPidRatioHelper::fgkfHistNBinsPt};
0a28d543 195
f7ea34d2 196 Double_t minHnCont[8] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], -0.5,
197 AliEbyEPidRatioHelper::fgkfHistRangeSign[0], 0.5,
198 AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
199 AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
200 AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
201 AliEbyEPidRatioHelper::fgkfHistRangePt[0]};
202
203 Double_t maxHnCont[8] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], 3.5,
204 AliEbyEPidRatioHelper::fgkfHistRangeSign[1], 8.5,
205 AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
206 AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
207 AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
208 AliEbyEPidRatioHelper::fgkfHistRangePt[1]};
209
210 fHnContMc = new THnSparseF("hnContMc", "cent:pid:SignMC:contPart:etaMC:yMC:phiMC:ptMC",8, binHnCont, minHnCont, maxHnCont);
211 fHnContRec = new THnSparseF("hnContRec", "cent:pid:SignRec:contPart:etaRec:yRec:phiRec:ptRec",8, binHnCont, minHnCont, maxHnCont);
212
213 fHnContMc->Sumw2();
214 fHnContRec->Sumw2();
215
216 fHnContMc->GetAxis(0)->SetTitle("centrality"); // 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
217 fHnContMc->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}"); // 0 | 1 | 2 | 3
218 fHnContMc->GetAxis(2)->SetTitle("sign"); // -1 | 0 | +1
219 fHnContMc->GetAxis(3)->SetTitle("contPart"); // 1 pi | 2 K | 3 p | 4 e | 5 mu | 6 other | 7 p from WeakDecay | 8 p from Material
220 fHnContMc->GetAxis(4)->SetTitle("#eta_{MC}"); // eta [-0.9,0.9]
221 fHnContMc->GetAxis(5)->SetTitle("#it{y}_{MC}"); // rapidity [-0.5, 0.5]
222 fHnContMc->GetAxis(6)->SetTitle("#varphi_{MC} (rad)"); // phi [ 0. ,2Pi]
223 fHnContMc->GetAxis(7)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})"); // pT [ 0.2,2.3]
0a28d543 224
0a28d543 225
f7ea34d2 226 fHnContRec->GetAxis(0)->SetTitle("centrality"); // 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
227 fHnContRec->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}"); // 0 | 1 | 2 | 3
228 fHnContRec->GetAxis(2)->SetTitle("sign"); // -1 | 0 | +1
229 fHnContRec->GetAxis(3)->SetTitle("contPart"); // 1 pi | 2 K | 3 p | 4 e | 5 mu | 6 other | 7 p from WeakDecay | 8 p from Material
230 fHnContRec->GetAxis(4)->SetTitle("#eta_{Rec}"); // eta [-0.9, 0.9]
231 fHnContRec->GetAxis(5)->SetTitle("#it{y}_{Rec}"); // rapidity [-0.5, 0.5]
232 fHnContRec->GetAxis(6)->SetTitle("#varphi_{Rec} (rad)"); // phi [ 0. , 2Pi]
233 fHnContRec->GetAxis(7)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 0.2, 2.3]
0a28d543 234
f7ea34d2 235
236 // fHnCont->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}"); // eta [-0.9, 0.9]
237 // fHnCont->GetAxis(13)->SetTitle("#it{y}_{MC}-#it{y}_{Rec}"); // rapidity [-0.5, 0.5]
238 // fHnCont->GetAxis(14)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); // phi [ -2Pi , 2Pi]
239 // fHnCont->GetAxis(15)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ -2.3, 2.3]
240 // fHnCont->GetAxis(16)->SetTitle("sign_{MC}-sign_{Rec}"); // -2 | 0 | +2
241 // fHnCont->GetAxis(17)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}"); // 0 | 1 | 2 | 3
242
6ce4ad92 243 fHelper->BinLogAxis(fHnContMc, 7);
244 fHelper->BinLogAxis(fHnContRec, 7);
0a28d543 245
246 return;
247}
248
249//________________________________________________________________________
250Int_t AliEbyEPidRatioEffCont::Setup() {
251 // -- Setup eventwise
252
253 // -- Create label arrays
254 for(Int_t i = 0; i < 4; i++) {
255 fLabelsRec[0][i] = new Int_t[fNTracks];
256 if(!fLabelsRec[0][i]) {
257 AliError("Cannot create fLabelsRec[0]");
258 return -1;
259 }
260
261 fLabelsRec[1][i] = new Int_t[fNTracks];
262 if(!fLabelsRec[1][i]) {
263 AliError("Cannot create fLabelsRec[1] for PID");
264 return -1;
265 }
266
267 for(Int_t ii = 0; ii < fNTracks; ++ii) {
268 fLabelsRec[0][i][ii] = 0;
269 fLabelsRec[1][i][ii] = 0;
270 }
271 }
272 return 0;
273}
274
275//________________________________________________________________________
276void AliEbyEPidRatioEffCont::Reset() {
277 // -- Reset eventwise
278 for(Int_t i = 0; i < 4; i++) {
279 for (Int_t ii = 0; ii < 2 ; ++ii) {
280 if (fLabelsRec[ii][i])
281 delete[] fLabelsRec[ii][i];
282 fLabelsRec[ii][i] = NULL;
283 }
284 }
285}
286
287//________________________________________________________________________
288void AliEbyEPidRatioEffCont::FillMCLabels() {
289 Float_t etaRange[2];
290 fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
291
292 Float_t ptRange[2];
293 fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
294
295 // -- Track Loop
296 for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
297
298 AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
299
300 // -- Check if track is accepted for basic parameters
301 if (!fHelper->IsTrackAcceptedBasicCharged(track))
302 continue;
303
304 // -- Check if accepted - ESD
305 if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
306 continue;
307
308 // -- Check if accepted - AOD
309 if (fAOD){
310 AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
311
312 if (!trackAOD) {
313 AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
314 continue;
315 }
316 if (!trackAOD->TestFilterBit(fAODtrackCutBit))
317 continue;
318
319 // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
320 if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
321 continue;
322 }
323
324 Int_t gPdgCode = 0;
325
326 Int_t iPid = 0;
327 Double_t pid[3];
328 if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion))) { iPid = 1; gPdgCode = 211;}
329 else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon))) { iPid = 2; gPdgCode = 321;}
330 else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))){ iPid = 3; gPdgCode = 2212;}
331 else iPid = 0;
332
15bb9247 333 // cout << " --- EFF ---- " << iPid << " " << gPdgCode << endl;
0a28d543 334
335 Double_t yP;
336 if (!fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
337 continue;
338
339 if (!fHelper->IsTrackAcceptedDCA(track))
340 continue;
341
342 Int_t label = TMath::Abs(track->GetLabel());
343
344 // -- Fill Label of all reconstructed
345 if(iPid != 0) fLabelsRec[0][0][idxTrack] = label;
346 fLabelsRec[0][iPid][idxTrack] = label;
347
348 // -- Fill Label of all reconstructed && recPid_TPC+TOF
349 if(iPid != 0) fLabelsRec[1][0][idxTrack] = label;
350 fLabelsRec[1][iPid][idxTrack] = label;
351
352 // -- Check for contamination and fill contamination THnSparse
353 CheckContTrack(track, iPid, gPdgCode);
354 if(iPid != 0) CheckContTrack(track, 0, 0);
355
356 } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
357
358 return;
359}
360
361//________________________________________________________________________
362void AliEbyEPidRatioEffCont::CheckContTrack(AliVTrack *track, Int_t iPid, Int_t gPdgCode) {
363 Int_t label = TMath::Abs(track->GetLabel());
364 Float_t signRec = track->Charge();
365
366
367 AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
368 if (!particle)
369 return;
370
371 Bool_t isPhysicalPrimary = (fESD) ? fStack->IsPhysicalPrimary(label): (static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary();
372 if (iPid == 0) {
373 if (particle->PdgCode() == (signRec*gPdgCode))
374 if (isPhysicalPrimary)
375 return;
376 }
377 else {
378 if (isPhysicalPrimary)
379 return;
380 }
381
382 // -- Check if secondaries from material or weak decay
383 Bool_t isSecondaryFromWeakDecay = (fESD) ? fStack->IsSecondaryFromWeakDecay(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromWeakDecay();
384 Bool_t isSecondaryFromMaterial = (fESD) ? fStack->IsSecondaryFromMaterial(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromMaterial();
385
386 // -- Get PDG Charge of contaminating particle
387 Float_t signMC = 0.;
388 if (particle->Charge() == 0.) signMC = 0.;
389 else if (particle->Charge() < 0.) signMC = -1.;
390 else if (particle->Charge() > 0.) signMC = 1.;
391
392 // -- Get contaminating particle
393 Float_t contPart = 0;
394 if (isSecondaryFromWeakDecay) contPart = 7; // probeParticle from WeakDecay
395 else if (isSecondaryFromMaterial) contPart = 8; // probeParticle from Material
396 else {
397 if (TMath::Abs(particle->PdgCode()) == 211) contPart = 1; // pion
398 else if (TMath::Abs(particle->PdgCode()) == 321) contPart = 2; // kaon
399 else if (TMath::Abs(particle->PdgCode()) == 2212) contPart = 3; // proton
400 else if (TMath::Abs(particle->PdgCode()) == 11) contPart = 4; // electron
401 else if (TMath::Abs(particle->PdgCode()) == 13) contPart = 5; // muon
402 else contPart = 6; // other
403 }
404
15bb9247 405 // cout << " --- CONT ---- " << iPid << " " << gPdgCode << endl;
406
0a28d543 407 // -- Get Reconstructed y
408 // yRec = y for identified particles | yRec = eta for charged particles
409 Double_t yRec = 0.;
410 fHelper->IsTrackAcceptedRapidity(track, yRec, iPid);
411
412 Double_t deltaPhi = particle->Phi()-track->Phi();
413 if (TMath::Abs(deltaPhi) > TMath::TwoPi()) {
414 if (deltaPhi < 0)
415 deltaPhi += TMath::TwoPi();
416 else
417 deltaPhi -= TMath::TwoPi();
418 }
419
f7ea34d2 420 Double_t hnContMc[8] = {fCentralityBin,iPid,signMC,contPart,particle->Eta(),particle->Y(),particle->Phi(),particle->Pt()};
421 Double_t hnContRec[8] = {fCentralityBin,iPid,signRec,contPart, track->Eta(),yRec,track->Phi(),track->Pt()};
422 fHnContMc->Fill(hnContMc);
423 fHnContRec->Fill(hnContRec);
0a28d543 424
425}
426
427//________________________________________________________________________
428void AliEbyEPidRatioEffCont::FillMCEffHist() {
429 // Fill efficiency THnSparse for ESDs
430
431 Float_t etaRange[2];
432 fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
433
434 Int_t nPart = (fESD) ? fStack->GetNprimary() : fArrayMC->GetEntriesFast();
435
436 for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
437 AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : static_cast<AliVParticle*>(fArrayMC->At(idxMC));
438
439 // -- Check basic MC properties -> charged physical primary
440 if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
441 continue;
442
443 Int_t iPid = 0;
444 Int_t gPdgCode = 0;
eda5b3c1 445 if ( particle->PdgCode() == 211 ) { iPid = 1; gPdgCode = 211;}
446 else if ( particle->PdgCode() == 321 ) { iPid = 1; gPdgCode = 321;}
447 else if ( particle->PdgCode() == 2212 ) { iPid = 1; gPdgCode = 2212;}
0a28d543 448 else {iPid = 0; gPdgCode = 0;}
449
450 // -- Check if accepted in rapidity window -- for identified particles
451 Double_t yMC;
452 if (iPid != 0 && !fHelper->IsParticleAcceptedRapidity(particle, yMC, iPid))
453 continue;
454
455 // -- Check if accepted in eta window -- for charged particles
456 if (iPid == 0 && TMath::Abs(particle->Eta()) > etaRange[1])
457 continue;
458
459 // -- Check if probeParticle / anti-probeParticle
460 // > skip check if PID is not required
461 if (iPid == 0 && TMath::Abs(particle->PdgCode()) != gPdgCode)
462 continue;
463
464 // -- Get sign of particle
465 Float_t signMC = (particle->PdgCode() < 0) ? -1. : 1.;
466
467 // -- Get if particle is findable --- not availible for AODs yet
468 Float_t findable = (fESD) ? Float_t(fHelper->IsParticleFindable(idxMC)) : 1.;
469
470 // -- Get recStatus and pidStatus
471 Float_t recStatus = 0.;
472 Float_t recPid = 0.;
473
474 // -- Get Reconstructed values
475 Float_t etaRec = 0.;
476 Float_t phiRec = 0.;
477 Float_t ptRec = 0.;
478 Double_t yRec = 0.;
479 Float_t signRec = 0.;
480
481 // -- Loop over all labels
482 for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
483 if (idxMC == fLabelsRec[0][iPid][idxRec]) {
484 recStatus = 1.;
485
486 if (idxMC == fLabelsRec[1][iPid][idxRec])
487 recPid = 1.;
488
489 AliVTrack *track = NULL;
490 if(fESD)
491 track = fESD->GetTrack(idxRec);
492 else if(fAOD)
493 track = fAOD->GetTrack(idxRec);
494
495 if (track) {
496 // if no track present (which should not happen)
497 // -> pt = 0. , which is not in the looked at range
498
499 // -- Get Reconstructed values
500 etaRec = track->Eta();
501 phiRec = track->Phi();
502 ptRec = track->Pt();
503 signRec = track->Charge();
504 fHelper->IsTrackAcceptedRapidity(track, yRec, iPid); // yRec = y for identified particles | yRec = eta for charged particles
505 }
506 break;
507 }
508 } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
509
510 Double_t deltaPhi = particle->Phi()-phiRec;
511 if (TMath::Abs(deltaPhi) > TMath::TwoPi()) {
512 if (deltaPhi < 0)
513 deltaPhi += TMath::TwoPi();
514 else
515 deltaPhi -= TMath::TwoPi();
516 }
517
518 if(iPid != 0) {
f7ea34d2 519 Double_t hnEffMc[10] = {fCentralityBin,0,signMC,findable, recStatus,recPid,particle->Eta(), particle->Y(), particle->Phi(),particle->Pt()};
520 Double_t hnEffRec[10] = {fCentralityBin,0,signRec,findable, recStatus,recPid,etaRec, yRec, phiRec, ptRec};
521 fHnEffMc->Fill(hnEffMc);
522 fHnEffRec->Fill(hnEffRec);
0a28d543 523 }
f7ea34d2 524 Double_t hnEffMc[10] = {fCentralityBin,iPid,signMC,findable, recStatus,recPid,particle->Eta(), particle->Y(), particle->Phi(),particle->Pt()};
525 Double_t hnEffRec[10] = {fCentralityBin,iPid,signRec,findable, recStatus,recPid,etaRec, yRec, phiRec, ptRec};
526 fHnEffMc->Fill(hnEffMc);
527 fHnEffRec->Fill(hnEffRec);
0a28d543 528
f7ea34d2 529
0a28d543 530
531 } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
532
533 return;
534}