Net Particle fluctuation task including efficiency correction (by Jochen Thaeder)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / NetParticle / AliAnalysisNetParticleDCA.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 "AliAnalysisNetParticleDCA.h"
14
15 using namespace std;
16
17 // Task for NetParticle checks
18 // Author: Jochen Thaeder <jochen@thaeder.de>
19
20 ClassImp(AliAnalysisNetParticleDCA)
21
22 /*
23  * ---------------------------------------------------------------------------------
24  *                            Constructor / Destructor
25  * ---------------------------------------------------------------------------------
26  */
27
28 //________________________________________________________________________
29 AliAnalysisNetParticleDCA::AliAnalysisNetParticleDCA() :
30   fHelper(NULL),
31
32   fESD(NULL), 
33   fESDTrackCuts(NULL),
34   fESDTrackCutsBkg(NULL),
35
36   fStack(NULL),
37   fMCEvent(NULL),
38
39   fHnDCA(NULL) {
40   // Constructor   
41
42   AliLog::SetClassDebugLevel("AliAnalysisNetParticleDCA",10);
43 }
44
45 //________________________________________________________________________
46 AliAnalysisNetParticleDCA::~AliAnalysisNetParticleDCA() {
47   // Destructor
48
49 }
50
51
52
53 /*
54  * ---------------------------------------------------------------------------------
55  *                                 Public Methods
56  * ---------------------------------------------------------------------------------
57  */
58
59 //________________________________________________________________________
60 void AliAnalysisNetParticleDCA::Initialize(AliESDtrackCuts *cuts, AliESDtrackCuts *cutsBkg, AliAnalysisNetParticleHelper* helper) {
61
62   // -- Get Helper
63   // ---------------
64   fHelper          = helper;
65
66   // -- ESD track cuts
67   // -------------------
68   fESDTrackCuts    = cuts;
69   fESDTrackCutsBkg = cutsBkg;
70
71 #if 0
72   // -- Get particle species / pdgCode
73   // -------------------------
74   fPdgCode          = AliPID::ParticleCode(fHelper->GetParticleSpecies());
75
76   // -- MC Labels for efficiency
77   // -----------------------------
78   fLabelsRec        = new Int_t*[2];
79   for (Int_t ii = 0; ii < 2 ; ++ii)
80     fLabelsRec[ii] = NULL;
81
82   // -- Create THnSparse Histograms
83   // --------------------------------
84   CreateHistograms();
85 #endif  
86   return;
87 }
88
89 /*
90  * ---------------------------------------------------------------------------------
91  *                            Setup/Reset Methods - private
92  * ---------------------------------------------------------------------------------
93  */
94
95 //________________________________________________________________________
96 Int_t AliAnalysisNetParticleDCA::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
97   // -- Setup event
98
99   // -- Reset of event
100   // -------------------
101   ResetEvent();
102
103   // -- Setup of event
104   // -------------------
105   fESD           = esdHandler->GetEvent();
106
107   fMCEvent       = mcEvent;
108   fStack         = fMCEvent->Stack();
109 #if 0
110   fCentralityBin = centBin;
111
112   // -- Create label arrays
113   // ------------------------
114   fLabelsRec[0] = new Int_t[fNTracks];
115   if(!fLabelsRec[0]) {
116     AliError("Cannot create fLabelsRec[0]");
117     return -1;
118   }
119
120   fLabelsRec[1] = new Int_t[fNTracks];
121   if(!fLabelsRec[1]) {
122     AliError("Cannot create fLabelsRec[1] for TPC");
123     return -1;
124   }
125
126   for(Int_t ii = 0; ii < fNTracks; ++ii) {
127     fLabelsRec[0][ii] = 0;
128     fLabelsRec[1][ii] = 0;
129   }
130 #endif
131   return 0;
132 }
133
134 //________________________________________________________________________
135 void AliAnalysisNetParticleDCA::ResetEvent() {
136   // -- Reset event
137   /*
138   for (Int_t ii = 0; ii < 2 ; ++ii) {
139     if (fLabelsRec[ii])
140       delete[] fLabelsRec[ii];
141     fLabelsRec[ii] = NULL;
142   }
143   */
144   return;
145 }
146
147 //________________________________________________________________________
148 void AliAnalysisNetParticleDCA::Process() {
149   // -- Process event
150
151   // -- Setup (clean, create and fill) MC labels
152   // ---------------------------------------------
153   //  FillMCLabels();
154
155   // -- Fill  MC histograms for efficiency studies
156   // -----------------------------------------------
157   //  FillMCEffHist();
158
159   return;
160 }      
161
162 /*
163  * ---------------------------------------------------------------------------------
164  *                                 Private Methods
165  * ---------------------------------------------------------------------------------
166  */
167 #if 0
168 //________________________________________________________________________
169 void AliAnalysisNetParticleDCA::CreateHistograms() {
170   // -- Create histograms
171
172   Double_t dCent[2] = {-0.5, 8.5};
173   Int_t iCent       = 9;
174   
175   Double_t dEta[2]  = {-0.9, 0.9}; 
176   Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ; 
177
178   Double_t dRap[2]  = {-0.5, 0.5}; 
179   Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ; 
180
181   Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
182   Int_t iPhi        = 42;
183
184   Double_t dPt[2]   = {0.1, 3.0};
185   Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
186
187   Double_t dSign[2] = {-1.5, 1.5};
188   Int_t iSign       = 3;
189
190   // ------------------------------------------------------------------
191   // -- Create THnSparseF - DCA
192   // ------------------------------------------------------------------
193
194   //                              cent:   etaMC:     yMC:   phiMC:   ptMC:     sign:contPart:DCArAccepted:DCAr
195   Int_t    binHnDCA[9] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       4,           2,  50 };     
196   Double_t minHnDCA[9] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0],     0.5,        -0.5, -10.};
197   Double_t maxHnDCA[9] = {dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1],     4.5,         1.5,  10.};
198   fHnDCA = new THnSparseF("fHnDCA", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:DCArAccepted:DCAr", 9, binHnDCA, minHnDCA, maxHnDCA);
199   
200   fHnDCA->Sumw2();
201   fHnDCA->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
202   fHnDCA->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9,0.9]
203   fHnDCA->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
204   fHnDCA->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. ,2Pi]
205   fHnDCA->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pT   [ 0.1,1.3]
206   fHnDCA->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1 
207   fHnDCA->GetAxis(6)->SetTitle("contPart");                     //  1  primary | 2 missId | 3 from WeakDecay | 4 p from Material
208   fHnDCA->GetAxis(7)->SetTitle("DCArAccepted");                 //  0 not accepted | 1 accepted 
209   fHnDCA->GetAxis(8)->SetTitle("DCAr");                         //  DCAr [-10,10]
210
211   // ------------------------------------------------------------------
212   
213   return;
214 }
215
216 //________________________________________________________________________
217 void AliAnalysisNetParticleDCA::FillMCDCA() {
218   // Fill MC labels
219   // Loop over ESD tracks and fill arrays with MC lables
220   //  fLabelsRec[0] : all Tracks
221   //  fLabelsRec[1] : all Tracks accepted by PID of TPC
222   // Check every accepted track if correctly identified
223   //  otherwise check for contamination
224
225   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
226     AliESDtrack *track = fESD->GetTrack(idxTrack); 
227     
228     // -- Check if track is accepted for basic parameters
229     if (!fHelper->IsTrackAcceptedBasicCharged(track))
230       continue;
231     
232     // -- Check if accepted
233     if (!fESDTrackCuts->AcceptTrack(track)) 
234       continue;
235
236     // -- Check if accepted in rapidity window
237     Double_t yP;
238     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
239       continue;
240
241     // -- Check if accepted with thighter DCA cuts
242     if (!fHelper->IsTrackAcceptedDCA(track))
243       continue;
244
245     Int_t label  = TMath::Abs(track->GetLabel()); 
246     
247     // -- Fill Label of all reconstructed
248     fLabelsRec[0][idxTrack] = label;
249
250     // -- Check if accepted by PID from TPC or TPC+TOF
251     Double_t pid[2];
252     if (!fHelper->IsTrackAcceptedPID(track, pid))
253       continue;
254
255     // -- Fill Label of all reconstructed && recPid_TPC+TOF    
256     fLabelsRec[1][idxTrack] = label;    
257     
258     // -- Check for contamination and fill contamination THnSparse
259     CheckContTrack(label, track->GetSign(), idxTrack);
260
261   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
262
263   return;
264 }
265
266 //________________________________________________________________________
267 void AliAnalysisNetParticleDCA::CheckDCATrack(Int_t label, Float_t sign, Float_t isDCArAccepted, Float_t dcar) {
268   // Check if DCA of contamination or correctly identified
269   // Fill contamination DCA THnSparse
270
271   TParticle* particle = fStack->Particle(label);
272   if (!particle)
273     return;
274
275   Bool_t isSecondaryFromWeakDecay = kFALSE;
276   Bool_t isSecondaryFromMaterial  = kFALSE;
277   
278   Int_t contPart = 0;
279
280   // -- Check if correctly identified 
281   if (particle->GetPdgCode() == (sign*fPdgCode)) {
282
283     // Check if is physical primary -> all ok 
284     if (fStack->IsPhysicalPrimary(label))
285       contPart = 1;    
286     // -- Check if secondaries from weak decay
287     else if(isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label))
288       contPart = 3;
289     // -- Check if secondaries from material decay
290     else if (isSecondaryFromMaterial  = fStack->IsSecondaryFromMaterial(label))
291       contPart = 4;
292   } 
293   // -- MissIdentification
294   else
295     contPart = 2;
296   
297   Double_t hnDCA[9] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, contPart, isDCArAccepted, dcar};
298   fHnDCA->Fill(hnContDCA);
299   
300   return;
301 }
302
303
304 #endif