Net Particle fluctuation task including efficiency correction (by Jochen Thaeder)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / NetParticle / AliAnalysisNetParticleEffCont.cxx
1 //-*- Mode: C++ -*-
2
3 #include "TMath.h"
4 #include "TAxis.h"
5
6 #include "AliESDEvent.h"
7 #include "AliESDInputHandler.h"
8 #include "AliStack.h"
9 #include "AliMCEvent.h"
10
11 #include "AliESDtrackCuts.h"
12
13 #include "AliAnalysisNetParticleEffCont.h"
14
15 using namespace std;
16
17 // Task for NetParticle checks
18 // Author: Jochen Thaeder <jochen@thaeder.de>
19
20 ClassImp(AliAnalysisNetParticleEffCont)
21
22 /*
23  * ---------------------------------------------------------------------------------
24  *                            Constructor / Destructor
25  * ---------------------------------------------------------------------------------
26  */
27
28 //________________________________________________________________________
29 AliAnalysisNetParticleEffCont::AliAnalysisNetParticleEffCont() :
30   fHelper(NULL),
31   fPdgCode(2212),
32
33   fESD(NULL), 
34   fESDTrackCuts(NULL),
35
36   fCentralityBin(-1.),
37   fNTracks(0),
38
39   fStack(NULL),
40   fMCEvent(NULL),
41
42   fLabelsRec(NULL),
43
44   fHnEff(NULL),
45   fHnCont(NULL) {
46   // Constructor   
47
48   AliLog::SetClassDebugLevel("AliAnalysisNetParticleEffCont",10);
49 }
50
51 //________________________________________________________________________
52 AliAnalysisNetParticleEffCont::~AliAnalysisNetParticleEffCont() {
53   // Destructor
54
55   if (fLabelsRec[0])
56     delete[] (fLabelsRec[0]);
57   if (fLabelsRec[1])
58     delete[] (fLabelsRec[1]);
59   if (fLabelsRec)
60     delete[] fLabelsRec;
61 }
62
63 /*
64  * ---------------------------------------------------------------------------------
65  *                                 Public Methods
66  * ---------------------------------------------------------------------------------
67  */
68
69 //________________________________________________________________________
70 void AliAnalysisNetParticleEffCont::Initialize(AliESDtrackCuts *cuts , AliAnalysisNetParticleHelper* helper) {
71
72   // -- ESD track cuts
73   // -------------------
74   fESDTrackCuts     = cuts;
75
76   // -- Get Helper
77   // ---------------
78   fHelper           = helper;
79
80   // -- Get particle species / pdgCode
81   // -------------------------
82   fPdgCode          = AliPID::ParticleCode(fHelper->GetParticleSpecies());
83
84   // -- MC Labels for efficiency
85   // -----------------------------
86   fLabelsRec        = new Int_t*[2];
87   for (Int_t ii = 0; ii < 2 ; ++ii)
88     fLabelsRec[ii] = NULL;
89
90   // -- Create THnSparse Histograms
91   // --------------------------------
92   CreateHistograms();
93   
94   return;
95 }
96
97 /*
98  * ---------------------------------------------------------------------------------
99  *                            Setup/Reset Methods - private
100  * ---------------------------------------------------------------------------------
101  */
102
103 //________________________________________________________________________
104 Int_t AliAnalysisNetParticleEffCont::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
105   // -- Setup event
106
107   // -- Reset of event
108   // -------------------
109   ResetEvent();
110
111   // -- Setup of event
112   // -------------------
113   fESD           = esdHandler->GetEvent();
114   fNTracks       = fESD->GetNumberOfTracks();
115
116   fMCEvent       = mcEvent;
117   fStack         = fMCEvent->Stack();
118
119   fCentralityBin = fHelper->GetCentralityBin();
120
121   // -- Create label arrays
122   // ------------------------
123   fLabelsRec[0] = new Int_t[fNTracks];
124   if(!fLabelsRec[0]) {
125     AliError("Cannot create fLabelsRec[0]");
126     return -1;
127   }
128
129   fLabelsRec[1] = new Int_t[fNTracks];
130   if(!fLabelsRec[1]) {
131     AliError("Cannot create fLabelsRec[1] for TPC");
132     return -1;
133   }
134
135   for(Int_t ii = 0; ii < fNTracks; ++ii) {
136     fLabelsRec[0][ii] = 0;
137     fLabelsRec[1][ii] = 0;
138   }
139
140   return 0;
141 }
142
143 //________________________________________________________________________
144 void AliAnalysisNetParticleEffCont::ResetEvent() {
145   // -- Reset event
146   
147   for (Int_t ii = 0; ii < 2 ; ++ii) {
148     if (fLabelsRec[ii])
149       delete[] fLabelsRec[ii];
150     fLabelsRec[ii] = NULL;
151   }
152
153   return;
154 }
155
156 //________________________________________________________________________
157 void AliAnalysisNetParticleEffCont::Process() {
158   // -- Process event
159
160   // -- Setup (clean, create and fill) MC labels
161   // ---------------------------------------------
162   FillMCLabels();
163
164   // -- Fill  MC histograms for efficiency studies
165   // -----------------------------------------------
166   FillMCEffHist();
167
168   return;
169 }      
170
171 /*
172  * ---------------------------------------------------------------------------------
173  *                                 Private Methods
174  * ---------------------------------------------------------------------------------
175  */
176
177 //________________________________________________________________________
178 void AliAnalysisNetParticleEffCont::CreateHistograms() {
179   // -- Create histograms
180
181   Double_t dCent[2] = {-0.5, 8.5};
182   Int_t iCent       = 9;
183   
184   Double_t dEta[2]  = {-0.9, 0.9}; 
185   Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ; 
186
187   Double_t dRap[2]  = {-0.5, 0.5}; 
188   Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ; 
189
190   Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
191   Int_t iPhi        = 42;
192
193   Double_t dPt[2]   = {0.1, 3.0};
194   Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
195
196   Double_t dSign[2] = {-1.5, 1.5};
197   Int_t iSign       = 3;
198
199   // ------------------------------------------------------------------
200   // -- Create THnSparseF - Eff
201   // ------------------------------------------------------------------
202
203   //                           cent:   etaMC:     yMC:  phiMC:   ptMC:     sign:findable:recStatus:pidStatus:   etaRec:  phiRec:  ptRec:deltaEta: deltaPhi: deltaPt
204   Int_t    binHnEff[15] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       2,      2  ,      2  ,    iEta,    iPhi,    iPt,    iEta, 2*iPhi+1, 2*iPt+1};
205   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]};
206   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]};
207
208   fHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt", 
209                           15, binHnEff, minHnEff, maxHnEff);
210
211   fHnEff->Sumw2();    
212
213   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
214   fHnEff->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9, 0.9]
215   fHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
216   fHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. , 2Pi]
217   fHnEff->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pt   [ 0.1, 3.0]
218   
219   fHnEff->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1 
220   fHnEff->GetAxis(6)->SetTitle("findable");                     //  0 not findable      |  1 findable
221   fHnEff->GetAxis(7)->SetTitle("recStatus");                    //  0 not reconstructed |  1 reconstructed
222   fHnEff->GetAxis(8)->SetTitle("recPid");                       //  0 not accepted      |  1 accepted
223
224   fHnEff->GetAxis(9)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
225   fHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)");         //  phi  [ 0. , 2Pi]
226   fHnEff->GetAxis(11)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 0.1, 3.0]
227
228   fHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}");                      //  eta  [-0.9,0.9]
229   fHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");          //  phi  [ 0. ,2Pi]
230   fHnEff->GetAxis(14)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 3.0,3.0]
231
232   fHelper->BinLogAxis(fHnEff,  4);
233   fHelper->BinLogAxis(fHnEff, 11);
234
235   /* 
236      >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
237      
238      creation -> findable -> reconstructed -> pid_TPC+TOF
239         (1)         (2)          (3)            (4)      
240                                                                       ||   findable | recStatus | recPid 
241      1) all primary probeParticles_MC                                 ||      -           -         -
242      2) all findable primary probeParticles_MC                        ||      x           -         -
243      3) all reconstructed primary probeParticles_MC                   ||      x           x         -
244      4) all reconstructed primary probeParticles_MC & recPid_TPC+TOF  ||      x           x         x
245      
246      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
247   */ 
248
249   // ------------------------------------------------------------------
250   // -- Create THnSparseF - Cont
251   // ------------------------------------------------------------------
252
253   //                            cent:   etaMC:    yMC:   phiMC:   ptMC:     sign:contPart: contSign:  etaRec:  phiRec:  ptRec:deltaEta: deltaPhi: deltaPt
254   Int_t    binHnCont[14] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       8,    iSign,    iEta,    iPhi,    iPt,    iEta, 2*iPhi+1, 2*iPt+1};
255   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]};
256   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]};
257   fHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",
258                            14, binHnCont, minHnCont, maxHnCont);
259
260   fHnCont->Sumw2();    
261   
262   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
263   fHnCont->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9,0.9]
264   fHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
265   fHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. ,2Pi]
266   fHnCont->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pT   [ 0.1,1.3]
267   fHnCont->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1 
268   
269   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
270   fHnCont->GetAxis(7)->SetTitle("contSign");                     //  -1 | 0 | +1 
271    
272   fHnCont->GetAxis(8)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
273   fHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)");          //  phi  [ 0. , 2Pi]
274   fHnCont->GetAxis(10)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 0.1, 3.0]
275
276   fHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}");                      //  eta  [-0.9,0.9]
277   fHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");          //  phi  [ 0. ,2Pi]
278   fHnCont->GetAxis(13)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); //  pt   [ 3.0,3.0]
279
280   fHelper->BinLogAxis(fHnCont,  4);
281   fHelper->BinLogAxis(fHnCont, 10);
282
283   // ------------------------------------------------------------------
284   
285   return;
286 }
287
288 //________________________________________________________________________
289 void AliAnalysisNetParticleEffCont::FillMCLabels() {
290   // Fill MC labels
291   // Loop over ESD tracks and fill arrays with MC lables
292   //  fLabelsRec[0] : all Tracks
293   //  fLabelsRec[1] : all Tracks accepted by PID of TPC
294   // Check every accepted track if correctly identified
295   //  otherwise check for contamination
296
297   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
298     AliESDtrack *track = fESD->GetTrack(idxTrack); 
299     
300     // -- Check if track is accepted for basic parameters
301     if (!fHelper->IsTrackAcceptedBasicCharged(track))
302       continue;
303     
304     // -- Check if accepted
305     if (!fESDTrackCuts->AcceptTrack(track)) 
306       continue;
307
308     // -- Check if accepted in rapidity window
309     Double_t yP;
310     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
311       continue;
312
313     // -- Check if accepted with thighter DCA cuts
314     if (!fHelper->IsTrackAcceptedDCA(track))
315       continue;
316
317     Int_t label  = TMath::Abs(track->GetLabel()); 
318     
319     // -- Fill Label of all reconstructed
320     fLabelsRec[0][idxTrack] = label;
321
322     // -- Check if accepted by PID from TPC or TPC+TOF
323     Double_t pid[2];
324     if (!fHelper->IsTrackAcceptedPID(track, pid))
325       continue;
326
327     // -- Fill Label of all reconstructed && recPid_TPC+TOF    
328     fLabelsRec[1][idxTrack] = label;    
329     
330     // -- Check for contamination and fill contamination THnSparse
331     CheckContTrack(label, track->GetSign(), idxTrack);
332
333   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
334
335   return;
336 }
337
338 //________________________________________________________________________
339 void AliAnalysisNetParticleEffCont::CheckContTrack(Int_t label, Float_t sign, Int_t idxTrack) {
340   // Check if particle is contamination or correctly identified
341   // Fill contamination THnSparse
342
343   TParticle* particle = fStack->Particle(label);
344   if (!particle)
345     return;
346
347   Bool_t isSecondaryFromWeakDecay = kFALSE;
348   Bool_t isSecondaryFromMaterial  = kFALSE;
349   
350   // -- Check if correctly identified 
351   if (particle->GetPdgCode() == (sign*fPdgCode)) {
352
353     // Check if is physical primary -> all ok 
354     if (fStack->IsPhysicalPrimary(label))
355       return;
356     
357     // -- Check if secondaries from material or weak decay
358     isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label);
359     isSecondaryFromMaterial  = fStack->IsSecondaryFromMaterial(label);
360   } 
361   
362   // -- Get MC pdg
363   Float_t contSign = 0.;
364   if      (particle->GetPDG()->Charge() == 0.) contSign =  0.;
365   else if (particle->GetPDG()->Charge() <  0.) contSign = -1.;  
366   else if (particle->GetPDG()->Charge() >  0.) contSign =  1.;  
367
368   // -- Get contaminating particle
369   Float_t contPart = 0;
370   if        (isSecondaryFromWeakDecay)                   contPart = 7; // probeParticle from WeakDecay
371   else if   (isSecondaryFromMaterial)                    contPart = 8; // probeParticle from Material
372   else {
373     if      (TMath::Abs(particle->GetPdgCode()) ==  211) contPart = 1; // pion
374     else if (TMath::Abs(particle->GetPdgCode()) ==  321) contPart = 2; // kaon
375     else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
376     else if (TMath::Abs(particle->GetPdgCode()) ==   11) contPart = 4; // electron
377     else if (TMath::Abs(particle->GetPdgCode()) ==   13) contPart = 5; // muon
378     else                                                 contPart = 6; // other
379   }
380   
381   // -- Get Reconstructed values 
382   Float_t etaRec = 0.;
383   Float_t phiRec = 0.;
384   Float_t ptRec  = 0.;
385
386   AliESDtrack *track = fESD->GetTrack(idxTrack);
387   if (track) {
388     // if no track present (which should not happen)
389     // -> pt = 0. , which is not in the looked at range
390     
391     // -- Get Reconstructed eta/phi/pt
392     etaRec = track->Eta();
393     phiRec = track->Phi();        
394     ptRec  = track->Pt();
395   }     
396
397   // -- Fill THnSparse
398   Double_t hnCont[14] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, 
399                          contPart, contSign, etaRec, phiRec, ptRec, 
400                          particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
401   fHnCont->Fill(hnCont);
402 }
403
404 //________________________________________________________________________
405 void AliAnalysisNetParticleEffCont::FillMCEffHist() {
406   // Fill efficiency THnSparse
407   
408   Int_t nPart  = fStack->GetNprimary();
409
410   for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
411     TParticle* particle = fStack->Particle(idxMC);
412
413     // -- Check basic MC properties -> charged physical primary
414     if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
415       continue;
416
417     // -- Check rapidity window
418     Double_t yMC;
419     if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
420       continue;
421
422     // -- Check if probeParticle / anti-probeParticle
423     if (TMath::Abs(particle->GetPdgCode()) != fPdgCode)
424       continue;
425     
426     // -- Get sign of particle
427     Float_t sign      = (particle->GetPdgCode() < 0) ? -1. : 1.;
428
429     // -- Get if particle is findable 
430     Float_t findable  = Float_t(fHelper->IsParticleFindable(idxMC));
431
432     // -- Get recStatus and pidStatus
433     Float_t recStatus = 0.;
434     Float_t recPid    = 0.;
435
436     // -- Get Reconstructed values 
437     Float_t etaRec = 0.;
438     Float_t phiRec = 0.;
439     Float_t ptRec  = 0.;
440
441     // -- Loop over all labels
442     for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
443       if (idxMC == fLabelsRec[0][idxRec]) {
444         recStatus = 1.;
445         
446         if (idxMC == fLabelsRec[1][idxRec])
447           recPid = 1.;
448
449         AliESDtrack *track = fESD->GetTrack(idxRec);
450         if (track) {
451           // if no track present (which should not happen)
452           // -> pt = 0. , which is not in the looked at range
453
454           // -- Get Reconstructed eta/phi/pt
455           etaRec = track->Eta();
456           phiRec = track->Phi();          
457           ptRec  = track->Pt();
458         }       
459         break;
460       }
461     } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {  
462     
463     // -- Fill THnSparse
464     Double_t hnEff[15] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, 
465                           findable, recStatus, recPid, etaRec, phiRec, ptRec, 
466                           particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
467     fHnEff->Fill(hnEff);
468
469   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
470   
471   return;
472 }