First update for AOD compatibility (distributions and efficiency), still missing...
[u/mrichter/AliRoot.git] / PWGCF / EBYE / NetParticle / AliAnalysisNetParticleEffCont.cxx
1 //-*- Mode: C++ -*-
2
3 #include "TMath.h"
4 #include "TAxis.h"
5 #include "TDatabasePDG.h"
6
7 #include "AliESDEvent.h"
8 #include "AliESDInputHandler.h"
9 #include "AliStack.h"
10 #include "AliMCEvent.h"
11 #include "AliESDtrackCuts.h"
12 #include "AliAODEvent.h"
13 #include "AliAODInputHandler.h"
14 #include "AliAODMCParticle.h"
15
16 #include "AliAnalysisNetParticleEffCont.h"
17
18 using namespace std;
19
20 // Task for NetParticle checks
21 // Author: Jochen Thaeder <jochen@thaeder.de>
22
23 ClassImp(AliAnalysisNetParticleEffCont)
24
25 /*
26  * ---------------------------------------------------------------------------------
27  *                            Constructor / Destructor
28  * ---------------------------------------------------------------------------------
29  */
30
31 //________________________________________________________________________
32 AliAnalysisNetParticleEffCont::AliAnalysisNetParticleEffCont() :
33   fHelper(NULL),
34   fPdgCode(2212),
35
36   fESD(NULL), 
37   fESDTrackCuts(NULL),
38
39   fAOD(NULL), 
40   fArrayMC(NULL),
41
42   fCentralityBin(-1.),
43   fNTracks(0),
44
45   fAODtrackCutBit(1024),
46
47   fStack(NULL),
48   fMCEvent(NULL),
49
50   fLabelsRec(NULL),
51
52   fHnEff(NULL),
53   fHnCont(NULL) {
54   // Constructor   
55
56   AliLog::SetClassDebugLevel("AliAnalysisNetParticleEffCont",10);
57 }
58
59 //________________________________________________________________________
60 AliAnalysisNetParticleEffCont::~AliAnalysisNetParticleEffCont() {
61   // Destructor
62
63   if (fLabelsRec[0])
64     delete[] (fLabelsRec[0]);
65   if (fLabelsRec[1])
66     delete[] (fLabelsRec[1]);
67   if (fLabelsRec)
68     delete[] fLabelsRec;
69 }
70
71 /*
72  * ---------------------------------------------------------------------------------
73  *                                 Public Methods
74  * ---------------------------------------------------------------------------------
75  */
76
77 //________________________________________________________________________
78 void AliAnalysisNetParticleEffCont::Initialize(AliESDtrackCuts *cuts , AliAnalysisNetParticleHelper* helper, Int_t trackCutBit) {
79
80   // -- ESD track cuts
81   // -------------------
82   fESDTrackCuts     = cuts;
83
84   // -- Get Helper
85   // ---------------
86   fHelper           = helper;
87
88   // -- Get particle species / pdgCode
89   // -------------------------
90   fPdgCode          = AliPID::ParticleCode(fHelper->GetParticleSpecies());
91
92   // -- MC Labels for efficiency
93   // -----------------------------
94   fLabelsRec        = new Int_t*[2];
95   for (Int_t ii = 0; ii < 2 ; ++ii)
96     fLabelsRec[ii] = NULL;
97
98   // -- AOD track filter bit
99   // -------------------------
100   fAODtrackCutBit = trackCutBit;
101
102
103   // -- Create THnSparse Histograms
104   // --------------------------------
105   CreateHistograms();
106   
107   return;
108 }
109
110 /*
111  * ---------------------------------------------------------------------------------
112  *                            Setup/Reset Methods - private
113  * ---------------------------------------------------------------------------------
114  */
115
116 //________________________________________________________________________
117 Int_t AliAnalysisNetParticleEffCont::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
118   // -- Setup event
119
120   // -- Reset of event
121   // -------------------
122   ResetEvent();
123
124   // -- Setup of event
125   // -------------------
126   fESD           = esdHandler->GetEvent();
127   fNTracks       = fESD->GetNumberOfTracks();
128
129   fMCEvent       = mcEvent;
130   fStack         = fMCEvent->Stack();
131
132   fCentralityBin = fHelper->GetCentralityBin();
133
134   // -- Create label arrays
135   // ------------------------
136   fLabelsRec[0] = new Int_t[fNTracks];
137   if(!fLabelsRec[0]) {
138     AliError("Cannot create fLabelsRec[0]");
139     return -1;
140   }
141
142   fLabelsRec[1] = new Int_t[fNTracks];
143   if(!fLabelsRec[1]) {
144     AliError("Cannot create fLabelsRec[1] for TPC");
145     return -1;
146   }
147
148   for(Int_t ii = 0; ii < fNTracks; ++ii) {
149     fLabelsRec[0][ii] = 0;
150     fLabelsRec[1][ii] = 0;
151   }
152
153   return 0;
154 }
155
156 //________________________________________________________________________
157 Int_t AliAnalysisNetParticleEffCont::SetupEvent(AliAODInputHandler *aodHandler) {
158   // -- Setup event
159
160   // -- Reset of event
161   // -------------------
162   ResetEvent();
163
164   // -- Setup of event
165   // -------------------
166   fAOD           = aodHandler->GetEvent();
167   fNTracks       = fAOD->GetNumberOfTracks();
168
169   fArrayMC       = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
170   if (!fArrayMC)
171     AliFatal("No array of MC particles found !!!");
172
173   fCentralityBin = fHelper->GetCentralityBin();
174
175   // -- Create label arrays
176   // ------------------------
177   fLabelsRec[0] = new Int_t[fNTracks];
178   if(!fLabelsRec[0]) {
179     AliError("Cannot create fLabelsRec[0]");
180     return -1;
181   }
182
183   fLabelsRec[1] = new Int_t[fNTracks];
184   if(!fLabelsRec[1]) {
185     AliError("Cannot create fLabelsRec[1] for TPC");
186     return -1;
187   }
188
189   for(Int_t ii = 0; ii < fNTracks; ++ii) {
190     fLabelsRec[0][ii] = 0;
191     fLabelsRec[1][ii] = 0;
192   }
193
194   return 0;
195 }
196
197 //________________________________________________________________________
198 void AliAnalysisNetParticleEffCont::ResetEvent() {
199   // -- Reset event
200   
201   for (Int_t ii = 0; ii < 2 ; ++ii) {
202     if (fLabelsRec[ii])
203       delete[] fLabelsRec[ii];
204     fLabelsRec[ii] = NULL;
205   }
206
207   return;
208 }
209
210 //________________________________________________________________________
211 void AliAnalysisNetParticleEffCont::Process() {
212   // -- Process event
213
214   // -- Setup (clean, create and fill) MC labels
215   // ---------------------------------------------
216   if(fESD)
217     FillMCLabels();
218   else if(fAOD)
219     FillMCLabelsAOD();
220   
221   // -- Fill  MC histograms for efficiency studies
222   // -----------------------------------------------
223   if(fESD)
224     FillMCEffHist();
225   else if(fAOD)
226     FillMCEffHistAOD();
227
228   return;
229 }      
230
231 /*
232  * ---------------------------------------------------------------------------------
233  *                                 Private Methods
234  * ---------------------------------------------------------------------------------
235  */
236
237 //________________________________________________________________________
238 void AliAnalysisNetParticleEffCont::CreateHistograms() {
239   // -- Create histograms
240
241   Double_t dCent[2] = {-0.5, 8.5};
242   Int_t iCent       = 9;
243   
244   Double_t dEta[2]  = {-0.9, 0.9}; 
245   Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ; 
246
247   Double_t dRap[2]  = {-0.5, 0.5}; 
248   Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ; 
249
250   Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
251   Int_t iPhi        = 42;
252
253   Double_t dPt[2]   = {0.1, 3.0};
254   Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
255
256   Double_t dSign[2] = {-1.5, 1.5};
257   Int_t iSign       = 3;
258
259   // ------------------------------------------------------------------
260   // -- Create THnSparseF - Eff
261   // ------------------------------------------------------------------
262
263   //                           cent:   etaMC:     yMC:  phiMC:   ptMC:     sign:findable:recStatus:pidStatus:   etaRec:  phiRec:  ptRec:deltaEta: deltaPhi: deltaPt
264   Int_t    binHnEff[15] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       2,      2  ,      2  ,    iEta,    iPhi,    iPt,    iEta, 2*iPhi+1, 2*iPt+1};
265   Double_t minHnEff[15] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0],    -0.5,     -0.5,     -0.5, dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
266   Double_t maxHnEff[15] = {dCent[1], dEta[1], dRap[0], dPhi[1], dPt[1], dSign[1],     1.5,      1.5,      1.5, dEta[1], dPhi[1], dPt[1], dEta[1],  dPhi[1],  dPt[1]};
267
268   fHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt", 
269                           15, binHnEff, minHnEff, maxHnEff);
270
271   fHnEff->Sumw2();    
272
273   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
274   fHnEff->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9, 0.9]
275   fHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
276   fHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. , 2Pi]
277   fHnEff->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pt   [ 0.1, 3.0]
278   
279   fHnEff->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1 
280   fHnEff->GetAxis(6)->SetTitle("findable");                     //  0 not findable      |  1 findable
281   fHnEff->GetAxis(7)->SetTitle("recStatus");                    //  0 not reconstructed |  1 reconstructed
282   fHnEff->GetAxis(8)->SetTitle("recPid");                       //  0 not accepted      |  1 accepted
283
284   fHnEff->GetAxis(9)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
285   fHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)");         //  phi  [ 0. , 2Pi]
286   fHnEff->GetAxis(11)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 0.1, 3.0]
287
288   fHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}");                      //  eta  [-0.9,0.9]
289   fHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");          //  phi  [ 0. ,2Pi]
290   fHnEff->GetAxis(14)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 3.0,3.0]
291
292   fHelper->BinLogAxis(fHnEff,  4);
293   fHelper->BinLogAxis(fHnEff, 11);
294
295   /* 
296      >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
297      
298      creation -> findable -> reconstructed -> pid_TPC+TOF
299         (1)         (2)          (3)            (4)      
300                                                                       ||   findable | recStatus | recPid 
301      1) all primary probeParticles_MC                                 ||      -           -         -
302      2) all findable primary probeParticles_MC                        ||      x           -         -
303      3) all reconstructed primary probeParticles_MC                   ||      x           x         -
304      4) all reconstructed primary probeParticles_MC & recPid_TPC+TOF  ||      x           x         x
305      
306      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
307   */ 
308
309   // ------------------------------------------------------------------
310   // -- Create THnSparseF - Cont
311   // ------------------------------------------------------------------
312
313   //                            cent:   etaMC:    yMC:   phiMC:   ptMC:     sign:contPart: contSign:  etaRec:  phiRec:  ptRec:deltaEta: deltaPhi: deltaPt
314   Int_t    binHnCont[14] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       8,    iSign,    iEta,    iPhi,    iPt,    iEta, 2*iPhi+1, 2*iPt+1};
315   Double_t minHnCont[14] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0],     0.5, dSign[0], dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
316   Double_t maxHnCont[14] = {dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1],     8.5, dSign[1], dEta[1], dPhi[1], dPt[1], dEta[1],  dPhi[1],  dPt[1]};
317   fHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",
318                            14, binHnCont, minHnCont, maxHnCont);
319
320   fHnCont->Sumw2();    
321   
322   fHnCont->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
323   fHnCont->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9,0.9]
324   fHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
325   fHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. ,2Pi]
326   fHnCont->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pT   [ 0.1,1.3]
327   fHnCont->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1 
328   
329   fHnCont->GetAxis(6)->SetTitle("contPart");                     //  1 pi | 2 K | 3 p | 4 e | 5 mu | 6 other | 7 p from WeakDecay | 8 p from Material
330   fHnCont->GetAxis(7)->SetTitle("contSign");                     //  -1 | 0 | +1 
331    
332   fHnCont->GetAxis(8)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
333   fHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)");          //  phi  [ 0. , 2Pi]
334   fHnCont->GetAxis(10)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 0.1, 3.0]
335
336   fHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}");                      //  eta  [-0.9,0.9]
337   fHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");          //  phi  [ 0. ,2Pi]
338   fHnCont->GetAxis(13)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 3.0,3.0]
339
340   fHelper->BinLogAxis(fHnCont,  4);
341   fHelper->BinLogAxis(fHnCont, 10);
342
343   // ------------------------------------------------------------------
344   
345   return;
346 }
347
348 //________________________________________________________________________
349 void AliAnalysisNetParticleEffCont::FillMCLabels() {
350   // Fill MC labels
351   // Loop over ESD tracks and fill arrays with MC lables
352   //  fLabelsRec[0] : all Tracks
353   //  fLabelsRec[1] : all Tracks accepted by PID of TPC
354   // Check every accepted track if correctly identified
355   //  otherwise check for contamination
356
357   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
358     AliESDtrack *track = fESD->GetTrack(idxTrack); 
359     
360     // -- Check if track is accepted for basic parameters
361     if (!fHelper->IsTrackAcceptedBasicCharged(track))
362       continue;
363     
364     // -- Check if accepted
365     if (!fESDTrackCuts->AcceptTrack(track)) 
366       continue;
367
368     // -- Check if accepted in rapidity window
369     Double_t yP;
370     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
371       continue;
372
373     // -- Check if accepted with thighter DCA cuts
374     if (!fHelper->IsTrackAcceptedDCA(track))
375       continue;
376
377     Int_t label  = TMath::Abs(track->GetLabel()); 
378     
379     // -- Fill Label of all reconstructed
380     fLabelsRec[0][idxTrack] = label;
381
382     // -- Check if accepted by PID from TPC or TPC+TOF
383     Double_t pid[2];
384     if (!fHelper->IsTrackAcceptedPID(track, pid))
385       continue;
386
387     // -- Fill Label of all reconstructed && recPid_TPC+TOF    
388     fLabelsRec[1][idxTrack] = label;    
389     
390     // -- Check for contamination and fill contamination THnSparse
391     CheckContTrack(label, track->GetSign(), idxTrack);
392
393   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
394
395   return;
396 }
397
398 //________________________________________________________________________
399 void AliAnalysisNetParticleEffCont::FillMCLabelsAOD() {
400   // Fill MC labels for AOD analysis
401   // Loop over ESD tracks and fill arrays with MC lables
402   //  fLabelsRec[0] : all Tracks
403   //  fLabelsRec[1] : all Tracks accepted by PID of TPC
404   // Check every accepted track if correctly identified
405   //  otherwise check for contamination
406
407   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
408     AliAODTrack *track = fAOD->GetTrack(idxTrack); 
409     
410     // -- Check if track is accepted for basic parameters
411     if (!fHelper->IsTrackAcceptedBasicCharged(track))
412       continue;
413     
414     // -- Check if accepted
415     if (!track->TestFilterBit(fAODtrackCutBit)) 
416       continue;
417
418     // -- Check if accepted in rapidity window
419     Double_t yP;
420     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
421       continue;
422
423     // -- Check if accepted with thighter DCA cuts (not yet done for AODs XXX)
424     //if (!fHelper->IsTrackAcceptedDCA(track))
425     //  continue;
426
427     Int_t label  = TMath::Abs(track->GetLabel()); 
428     
429     // -- Fill Label of all reconstructed
430     fLabelsRec[0][idxTrack] = label;
431
432     // -- Check if accepted by PID from TPC or TPC+TOF
433     Double_t pid[2];
434     if (!fHelper->IsTrackAcceptedPID(track, pid))
435       continue;
436
437     // -- Fill Label of all reconstructed && recPid_TPC+TOF    
438     fLabelsRec[1][idxTrack] = label;    
439
440     // -- Check for contamination and fill contamination THnSparse
441     CheckContTrackAOD(label, track->Charge(), idxTrack);
442
443   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
444
445   return;
446 }
447
448 //________________________________________________________________________
449 void AliAnalysisNetParticleEffCont::CheckContTrack(Int_t label, Float_t sign, Int_t idxTrack) {
450   // Check if particle is contamination or correctly identified for ESDs
451   // Fill contamination THnSparse
452   
453   TParticle* particle = fStack->Particle(label);
454   if (!particle)
455     return;
456
457   Bool_t isSecondaryFromWeakDecay = kFALSE;
458   Bool_t isSecondaryFromMaterial  = kFALSE;
459   
460   // -- Check if correctly identified 
461   if (particle->GetPdgCode() == (sign*fPdgCode)) {
462
463     // Check if is physical primary -> all ok 
464     if (fStack->IsPhysicalPrimary(label))
465       return;
466     
467     // -- Check if secondaries from material or weak decay
468     isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label);
469     isSecondaryFromMaterial  = fStack->IsSecondaryFromMaterial(label);
470   } 
471   
472   // -- Get MC pdg
473   Float_t contSign = 0.;
474   if      (particle->GetPDG()->Charge() == 0.) contSign =  0.;
475   else if (particle->GetPDG()->Charge() <  0.) contSign = -1.;  
476   else if (particle->GetPDG()->Charge() >  0.) contSign =  1.;  
477
478   // -- Get contaminating particle
479   Float_t contPart = 0;
480   if        (isSecondaryFromWeakDecay)                   contPart = 7; // probeParticle from WeakDecay
481   else if   (isSecondaryFromMaterial)                    contPart = 8; // probeParticle from Material
482   else {
483     if      (TMath::Abs(particle->GetPdgCode()) ==  211) contPart = 1; // pion
484     else if (TMath::Abs(particle->GetPdgCode()) ==  321) contPart = 2; // kaon
485     else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
486     else if (TMath::Abs(particle->GetPdgCode()) ==   11) contPart = 4; // electron
487     else if (TMath::Abs(particle->GetPdgCode()) ==   13) contPart = 5; // muon
488     else                                                 contPart = 6; // other
489   }
490   
491   // -- Get Reconstructed values 
492   Float_t etaRec = 0.;
493   Float_t phiRec = 0.;
494   Float_t ptRec  = 0.;
495
496   AliESDtrack *track = fESD->GetTrack(idxTrack);
497   if (track) {
498     // if no track present (which should not happen)
499     // -> pt = 0. , which is not in the looked at range
500     
501     // -- Get Reconstructed eta/phi/pt
502     etaRec = track->Eta();
503     phiRec = track->Phi();        
504     ptRec  = track->Pt();
505   }     
506
507   // -- Fill THnSparse
508   Double_t hnCont[14] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, 
509                          contPart, contSign, etaRec, phiRec, ptRec, 
510                          particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
511   fHnCont->Fill(hnCont);
512 }
513
514 //________________________________________________________________________
515 void AliAnalysisNetParticleEffCont::CheckContTrackAOD(Int_t label, Float_t sign, Int_t idxTrack) {
516   // Check if particle is contamination or correctly identified for AODs
517   // Fill contamination THnSparse
518   
519   AliAODMCParticle* particle = (AliAODMCParticle*)fArrayMC->At(label);
520
521   if (!particle)
522     return;
523
524   Bool_t isSecondaryFromWeakDecay = kFALSE;
525   Bool_t isSecondaryFromMaterial  = kFALSE;
526   
527   // -- Check if correctly identified 
528   if (particle->GetPdgCode() == (sign*fPdgCode)) {
529     
530     // Check if is physical primary -> all ok 
531     if(particle->IsPhysicalPrimary())
532       return;
533     
534     // -- Check if secondaries from material or weak decay
535     isSecondaryFromWeakDecay = particle->IsSecondaryFromWeakDecay();
536     isSecondaryFromMaterial  = particle->IsSecondaryFromMaterial();
537     
538   } 
539   
540   // -- Get MC pdg 
541   Float_t contSign = 0.;
542   if      ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() == 0.) contSign =  0.;
543   else if ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() <  0.) contSign = -1.;        
544   else if ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() >  0.) contSign =  1.;        
545
546   // -- Get contaminating particle
547   Float_t contPart = 0;
548   if        (isSecondaryFromWeakDecay)                   contPart = 7; // probeParticle from WeakDecay
549   else if   (isSecondaryFromMaterial)                    contPart = 8; // probeParticle from Material
550   else {
551     if      (TMath::Abs(particle->GetPdgCode()) ==  211) contPart = 1; // pion
552     else if (TMath::Abs(particle->GetPdgCode()) ==  321) contPart = 2; // kaon
553     else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
554     else if (TMath::Abs(particle->GetPdgCode()) ==   11) contPart = 4; // electron
555     else if (TMath::Abs(particle->GetPdgCode()) ==   13) contPart = 5; // muon
556     else                                                 contPart = 6; // other
557   }
558   
559   // -- Get Reconstructed values 
560   Float_t etaRec = 0.;
561   Float_t phiRec = 0.;
562   Float_t ptRec  = 0.;
563
564   AliAODTrack *track = fAOD->GetTrack(idxTrack);
565   
566   if (track) {
567     // if no track present (which should not happen)
568     // -> pt = 0. , which is not in the looked at range
569     
570     // -- Get Reconstructed eta/phi/pt
571     etaRec = track->Eta();
572     phiRec = track->Phi();        
573     ptRec  = track->Pt();
574   }     
575
576   // -- Fill THnSparse
577   Double_t hnCont[14] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, 
578                          contPart, contSign, etaRec, phiRec, ptRec, 
579                          particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
580   fHnCont->Fill(hnCont);
581 }
582
583 //________________________________________________________________________
584 void AliAnalysisNetParticleEffCont::FillMCEffHist() {
585   // Fill efficiency THnSparse for ESDs
586   
587   Int_t nPart  = fStack->GetNprimary();
588
589   for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
590     TParticle* particle = fStack->Particle(idxMC);
591
592     // -- Check basic MC properties -> charged physical primary
593     if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
594       continue;
595
596     // -- Check rapidity window
597     Double_t yMC;
598     if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
599       continue;
600
601     // -- Check if probeParticle / anti-probeParticle
602     if (TMath::Abs(particle->GetPdgCode()) != fPdgCode)
603       continue;
604     
605     // -- Get sign of particle
606     Float_t sign      = (particle->GetPdgCode() < 0) ? -1. : 1.;
607
608     // -- Get if particle is findable 
609     Float_t findable  = Float_t(fHelper->IsParticleFindable(idxMC));
610
611     // -- Get recStatus and pidStatus
612     Float_t recStatus = 0.;
613     Float_t recPid    = 0.;
614
615     // -- Get Reconstructed values 
616     Float_t etaRec = 0.;
617     Float_t phiRec = 0.;
618     Float_t ptRec  = 0.;
619
620     // -- Loop over all labels
621     for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
622       if (idxMC == fLabelsRec[0][idxRec]) {
623         recStatus = 1.;
624         
625         if (idxMC == fLabelsRec[1][idxRec])
626           recPid = 1.;
627
628         AliESDtrack *track = fESD->GetTrack(idxRec);
629         if (track) {
630           // if no track present (which should not happen)
631           // -> pt = 0. , which is not in the looked at range
632
633           // -- Get Reconstructed eta/phi/pt
634           etaRec = track->Eta();
635           phiRec = track->Phi();          
636           ptRec  = track->Pt();
637         }       
638         break;
639       }
640     } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {  
641     
642     // -- Fill THnSparse
643     Double_t hnEff[15] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, 
644                           findable, recStatus, recPid, etaRec, phiRec, ptRec, 
645                           particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
646     fHnEff->Fill(hnEff);
647
648   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
649   
650   return;
651 }
652
653 //________________________________________________________________________
654 void AliAnalysisNetParticleEffCont::FillMCEffHistAOD() {
655   // Fill efficiency THnSparse for AODs
656   
657   Int_t nPart  = fArrayMC->GetEntriesFast();
658
659   for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
660     
661     AliAODMCParticle* particle = (AliAODMCParticle*)fArrayMC->At(idxMC);
662              
663     // -- Check basic MC properties -> charged physical primary
664     if (!fHelper->IsParticleAcceptedBasicCharged(particle))
665       continue;
666
667     // -- Check rapidity window
668     Double_t yMC;
669     if (!fHelper->IsParticleAcceptedRapidity((TParticle*)particle, yMC))
670       continue;
671
672     // -- Check if probeParticle / anti-probeParticle
673     if (TMath::Abs(particle->GetPdgCode()) != fPdgCode)
674       continue;
675     
676     // -- Get sign of particle
677     Float_t sign      = (particle->GetPdgCode() < 0) ? -1. : 1.;
678
679     // -- Get if particle is findable 
680     Float_t findable  = 1.;// Float_t(fHelper->IsParticleFindable(idxMC)); XXX
681
682     // -- Get recStatus and pidStatus
683     Float_t recStatus = 0.;
684     Float_t recPid    = 0.;
685
686     // -- Get Reconstructed values 
687     Float_t etaRec = 0.;
688     Float_t phiRec = 0.;
689     Float_t ptRec  = 0.;
690
691     // -- Loop over all labels
692     for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
693       if (idxMC == fLabelsRec[0][idxRec]) {
694         recStatus = 1.;
695         
696         if (idxMC == fLabelsRec[1][idxRec])
697           recPid = 1.;
698
699         AliVTrack *track = NULL;
700         if(fESD)
701           track = fESD->GetTrack(idxRec);
702         else if(fAOD)
703           track = fAOD->GetTrack(idxRec);
704
705         if (track) {
706           // if no track present (which should not happen)
707           // -> pt = 0. , which is not in the looked at range
708
709           // -- Get Reconstructed eta/phi/pt
710           etaRec = track->Eta();
711           phiRec = track->Phi();          
712           ptRec  = track->Pt();
713         }       
714         break;
715       }
716     } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {  
717     
718     // -- Fill THnSparse
719     Double_t hnEff[15] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, 
720                           findable, recStatus, recPid, etaRec, phiRec, ptRec, 
721                           particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
722     fHnEff->Fill(hnEff);
723
724   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
725   
726   return;
727 }