coding violations corrected
[u/mrichter/AliRoot.git] / PWG4 / PartCorr / AliV0Reader.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt                        *\r
5  * Version 1.0                                                            *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 ////////////////////////////////////////////////\r
17 //--------------------------------------------- \r
18 // Class used to do analysis on conversion pairs\r
19 //---------------------------------------------\r
20 ////////////////////////////////////////////////\r
21 \r
22 // --- ROOT system ---\r
23 #include <TMath.h>\r
24 \r
25 //---- ANALYSIS system ----\r
26 #include "AliV0Reader.h"\r
27 #include "AliAnalysisManager.h"\r
28 #include "AliESDInputHandler.h"\r
29 #include "AliMCEvent.h"\r
30 #include "AliKFVertex.h"\r
31 \r
32 #include "AliStack.h"\r
33 #include "AliMCEventHandler.h"\r
34 \r
35 \r
36 class iostream;\r
37 class AliESDv0;\r
38 class TFormula;\r
39 \r
40 using namespace std;\r
41 \r
42 ClassImp(AliV0Reader)\r
43 \r
44 \r
45 \r
46   AliV0Reader::AliV0Reader() :\r
47     TObject(),\r
48     fMCStack(NULL),\r
49     fMCTruth(NULL),\r
50     fChain(NULL),\r
51     fESDHandler(NULL),\r
52     fESDEvent(NULL),\r
53     fHistograms(NULL),\r
54     fCurrentV0IndexNumber(0),\r
55     fCurrentV0(NULL),\r
56     fCurrentNegativeKFParticle(NULL),\r
57     fCurrentPositiveKFParticle(NULL),\r
58     fCurrentMotherKFCandidate(NULL),\r
59     fCurrentNegativeESDTrack(NULL),\r
60     fCurrentPositiveESDTrack(NULL),\r
61     fNegativeTrackLorentzVector(NULL),\r
62     fPositiveTrackLorentzVector(NULL),\r
63     fMotherCandidateLorentzVector(NULL),\r
64     fCurrentXValue(0),\r
65     fCurrentYValue(0),\r
66     fCurrentZValue(0),\r
67     fPositiveTrackPID(0),\r
68     fNegativeTrackPID(0),\r
69     fNegativeMCParticle(NULL),\r
70     fPositiveMCParticle(NULL),\r
71     fMotherMCParticle(NULL),\r
72     fMotherCandidateKFMass(0),\r
73     fMotherCandidateKFWidth(0),\r
74     fUseKFParticle(kTRUE),\r
75     fUseESDTrack(kFALSE),\r
76     fDoMC(kFALSE),\r
77     fMaxR(10000),// 100 meter(outside of ALICE)\r
78     fEtaCut(0.),\r
79     fPtCut(0.),\r
80     fChi2Cut(0.),\r
81     fPIDProbabilityCutNegativeParticle(0),\r
82     fPIDProbabilityCutPositiveParticle(0),\r
83     fXVertexCut(0.),\r
84     fYVertexCut(0.),\r
85     fZVertexCut(0.),\r
86     fNSigmaMass(0.),\r
87     fUseImprovedVertex(kFALSE),\r
88     fCurrentEventGoodV0s(),\r
89     fPreviousEventGoodV0s()\r
90 {\r
91 \r
92 }\r
93 \r
94 \r
95 AliV0Reader::AliV0Reader(const AliV0Reader & original) :\r
96   TObject(original),\r
97   fMCStack(original.fMCStack),\r
98   fMCTruth(original.fMCTruth),\r
99   fChain(original.fChain),\r
100   fESDHandler(original.fESDHandler),\r
101   fESDEvent(original.fESDEvent),\r
102   fHistograms(original.fHistograms),\r
103   fCurrentV0IndexNumber(original.fCurrentV0IndexNumber),\r
104   fCurrentV0(original.fCurrentV0),\r
105   fCurrentNegativeKFParticle(original.fCurrentNegativeKFParticle),\r
106   fCurrentPositiveKFParticle(original.fCurrentPositiveKFParticle),\r
107   fCurrentMotherKFCandidate(original.fCurrentMotherKFCandidate),\r
108   fCurrentNegativeESDTrack(original.fCurrentNegativeESDTrack),\r
109   fCurrentPositiveESDTrack(original.fCurrentPositiveESDTrack),\r
110   fNegativeTrackLorentzVector(original.fNegativeTrackLorentzVector),\r
111   fPositiveTrackLorentzVector(original.fPositiveTrackLorentzVector),\r
112   fMotherCandidateLorentzVector(original.fMotherCandidateLorentzVector),\r
113   fCurrentXValue(original.fCurrentXValue),\r
114   fCurrentYValue(original.fCurrentYValue),\r
115   fCurrentZValue(original.fCurrentZValue),\r
116   fPositiveTrackPID(original.fPositiveTrackPID),\r
117   fNegativeTrackPID(original.fNegativeTrackPID),\r
118   fNegativeMCParticle(original.fNegativeMCParticle),\r
119   fPositiveMCParticle(original.fPositiveMCParticle),\r
120   fMotherMCParticle(original.fMotherMCParticle),\r
121   fMotherCandidateKFMass(original.fMotherCandidateKFMass),\r
122   fMotherCandidateKFWidth(original.fMotherCandidateKFWidth),\r
123   fUseKFParticle(kTRUE),\r
124   fUseESDTrack(kFALSE),\r
125   fDoMC(kFALSE),\r
126   fMaxR(original.fMaxR),\r
127   fEtaCut(original.fEtaCut),\r
128   fPtCut(original.fPtCut),\r
129   fChi2Cut(original.fChi2Cut),\r
130   fPIDProbabilityCutNegativeParticle(original.fPIDProbabilityCutNegativeParticle),\r
131   fPIDProbabilityCutPositiveParticle(original.fPIDProbabilityCutPositiveParticle),\r
132   fXVertexCut(original.fXVertexCut),\r
133   fYVertexCut(original.fYVertexCut),\r
134   fZVertexCut(original.fZVertexCut),\r
135   fNSigmaMass(original.fNSigmaMass),\r
136   fUseImprovedVertex(original.fUseImprovedVertex),\r
137   fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),\r
138   fPreviousEventGoodV0s(original.fPreviousEventGoodV0s)\r
139 {\r
140 \r
141 }\r
142 \r
143 \r
144 AliV0Reader & AliV0Reader::operator = (const AliV0Reader & /*source*/)\r
145 {\r
146   // assignment operator\r
147   return *this;\r
148 }\r
149 \r
150 void AliV0Reader::Initialize(){\r
151   //see header file for documentation\r
152 \r
153   // Get the input handler from the manager\r
154   fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
155   if(fESDHandler == NULL){\r
156     //print warning here\r
157   }\r
158   \r
159   // Get pointer to esd event from input handler\r
160   fESDEvent = fESDHandler->GetEvent();\r
161   if(fESDEvent == NULL){\r
162     //print warning here\r
163   }\r
164 \r
165   //Get pointer to MCTruth\r
166   fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());\r
167   if(fMCTruth == NULL){\r
168     //print warning here\r
169   }\r
170 \r
171   //Get pointer to the mc stack\r
172   fMCStack = fMCTruth->MCEvent()->Stack();\r
173   if(fMCStack == NULL){\r
174     //print warning here\r
175   }\r
176 \r
177   AliKFParticle::SetField(fESDEvent->GetMagneticField());\r
178 \r
179 }\r
180 \r
181 AliESDv0* AliV0Reader::GetV0(Int_t index){\r
182   //see header file for documentation\r
183 \r
184   fCurrentV0 = fESDEvent->GetV0(index);\r
185   UpdateV0Information();\r
186   return fCurrentV0;\r
187 }\r
188 \r
189 \r
190 Bool_t AliV0Reader::NextV0(){\r
191   //see header file for documentation\r
192 \r
193   Bool_t iResult=kFALSE;\r
194   while(fCurrentV0IndexNumber<fESDEvent->GetNumberOfV0s()){\r
195     fCurrentV0 = fESDEvent->GetV0(fCurrentV0IndexNumber);\r
196     \r
197     //checks if on the fly mode is set\r
198     if ( !fCurrentV0->GetOnFlyStatus() ){\r
199       fCurrentV0IndexNumber++;\r
200       fHistograms->FillHistogram("V0MassDebugCut1",GetMotherCandidateMass());\r
201       continue;\r
202     }\r
203 \r
204     if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {//checks if we have a vertex\r
205       fCurrentV0IndexNumber++;\r
206       fHistograms->FillHistogram("V0MassDebugCut2",GetMotherCandidateMass());\r
207       continue;\r
208     }\r
209 \r
210     if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){\r
211       fCurrentV0IndexNumber++;\r
212       fHistograms->FillHistogram("V0MassDebugCut3",GetMotherCandidateMass());\r
213       continue;\r
214     }\r
215 \r
216 \r
217     fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);\r
218  \r
219     if(GetXYRadius()>fMaxR){ // cuts on distance from collision point\r
220       fCurrentV0IndexNumber++;\r
221       fHistograms->FillHistogram("V0MassDebugCut4",GetMotherCandidateMass());\r
222       continue;\r
223     }\r
224 \r
225     UpdateV0Information();\r
226         \r
227     if(fUseKFParticle){\r
228       if(fCurrentMotherKFCandidate->GetNDF()<=0){\r
229         fCurrentV0IndexNumber++;\r
230         fHistograms->FillHistogram("V0MassDebugCut5",GetMotherCandidateMass());\r
231         continue;\r
232       }\r
233       Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();\r
234       if(chi2V0 > fChi2Cut || chi2V0 <=0){\r
235         fCurrentV0IndexNumber++;\r
236         fHistograms->FillHistogram("V0MassDebugCut6",GetMotherCandidateMass());\r
237         continue;\r
238       }\r
239       \r
240       if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut){\r
241         fCurrentV0IndexNumber++;\r
242         fHistograms->FillHistogram("V0MassDebugCut7",GetMotherCandidateMass());\r
243         continue;\r
244       }\r
245       \r
246       if(fMotherCandidateLorentzVector->Pt()<fPtCut){\r
247         fCurrentV0IndexNumber++;\r
248         fHistograms->FillHistogram("V0MassDebugCut8",GetMotherCandidateMass());\r
249         continue;\r
250       }\r
251       \r
252     }\r
253     else if(fUseESDTrack){\r
254       //TODO\r
255     }\r
256 \r
257     iResult=kTRUE;//means we have a v0 who survived all the cuts applied\r
258 \r
259     fCurrentV0IndexNumber++;\r
260     \r
261     break;\r
262   }\r
263   return iResult; \r
264 }\r
265 \r
266 void AliV0Reader::UpdateV0Information(){\r
267   //see header file for documentation\r
268   \r
269   if(fCurrentNegativeKFParticle != NULL){\r
270     delete fCurrentNegativeKFParticle;\r
271   }\r
272   fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fNegativeTrackPID);\r
273   \r
274   if(fCurrentPositiveKFParticle != NULL){\r
275     delete fCurrentPositiveKFParticle;\r
276   }\r
277   fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fPositiveTrackPID);\r
278     \r
279   if(fCurrentMotherKFCandidate != NULL){\r
280     delete fCurrentMotherKFCandidate;\r
281   }\r
282   fCurrentMotherKFCandidate = new AliKFParticle(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);\r
283 \r
284   fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());\r
285 \r
286   fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());\r
287 \r
288   if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){\r
289     fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);\r
290   }\r
291 \r
292   if(fUseImprovedVertex == kTRUE){\r
293     AliKFVertex primaryVertexImproved(*GetPrimaryVertex());\r
294     primaryVertexImproved+=*fCurrentMotherKFCandidate;\r
295     fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);\r
296   }\r
297 \r
298   fCurrentMotherKFCandidate->GetMass(fMotherCandidateKFMass,fMotherCandidateKFWidth);\r
299 \r
300 \r
301   if(fNegativeTrackLorentzVector != NULL){\r
302     delete fNegativeTrackLorentzVector;\r
303   }\r
304   if(fUseKFParticle){\r
305     fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeKFParticle->Px(),fCurrentNegativeKFParticle->Py(),fCurrentNegativeKFParticle->Pz());\r
306   }\r
307   else if(fUseESDTrack){\r
308     fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeESDTrack->Px(),fCurrentNegativeESDTrack->Py(),fCurrentNegativeESDTrack->Pz());\r
309   }\r
310 \r
311   if(fPositiveTrackLorentzVector != NULL){\r
312     delete fPositiveTrackLorentzVector;\r
313   }\r
314   if(fUseKFParticle){\r
315     fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveKFParticle->Px(),fCurrentPositiveKFParticle->Py(),fCurrentPositiveKFParticle->Pz());\r
316   }\r
317   else if(fUseESDTrack){\r
318     fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveESDTrack->Px(),fCurrentPositiveESDTrack->Py(),fCurrentPositiveESDTrack->Pz());\r
319   }\r
320 \r
321   if(fMotherCandidateLorentzVector != NULL){\r
322     delete fMotherCandidateLorentzVector;\r
323   }\r
324   if(fUseKFParticle){\r
325     fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);\r
326   }\r
327   else if(fUseESDTrack){\r
328     fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);\r
329   }\r
330 \r
331   if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){\r
332     fMotherCandidateLorentzVector->SetXYZM(fMotherCandidateLorentzVector->Px() ,fMotherCandidateLorentzVector->Py(),fMotherCandidateLorentzVector->Pz(),0.); \r
333   }\r
334     \r
335   if(fDoMC){\r
336     fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));\r
337     fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));\r
338   }\r
339 }\r
340 \r
341 Bool_t AliV0Reader::HasSameMCMother(){\r
342   //see header file for documentation\r
343 \r
344   Bool_t iResult = kFALSE;\r
345   if(fDoMC == kTRUE){\r
346     if(fNegativeMCParticle != NULL && fPositiveMCParticle != NULL){\r
347       if(fNegativeMCParticle->GetMother(0) == fPositiveMCParticle->GetMother(0))\r
348         fMotherMCParticle = fMCStack->Particle(fPositiveMCParticle->GetMother(0));\r
349       iResult = kTRUE;\r
350     }\r
351   }\r
352   return iResult;\r
353 }\r
354 \r
355 Bool_t AliV0Reader::CheckPIDProbability(Double_t negProbCut, Double_t posProbCut){\r
356   //see header file for documentation\r
357 \r
358   Bool_t iResult=kFALSE;\r
359 \r
360   Double_t *posProbArray = new Double_t[10];\r
361   Double_t *negProbArray = new Double_t[10];\r
362   AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());\r
363   AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());\r
364   \r
365   negTrack->GetTPCpid(negProbArray);\r
366   posTrack->GetTPCpid(posProbArray);\r
367 \r
368   if(negProbArray!=NULL && posProbArray!=NULL){\r
369     if(negProbArray[GetSpeciesIndex(-1)]>=negProbCut && posProbArray[GetSpeciesIndex(1)]>=posProbCut){\r
370       iResult=kTRUE;\r
371     }\r
372   }\r
373   delete [] posProbArray;\r
374   delete [] negProbArray;\r
375   return iResult;\r
376 }\r
377 \r
378 void AliV0Reader::UpdateEventByEventData(){\r
379   //see header file for documentation\r
380 \r
381   if(fCurrentEventGoodV0s.size() >0 ){\r
382     fPreviousEventGoodV0s.clear();\r
383     fPreviousEventGoodV0s = fCurrentEventGoodV0s;\r
384   }\r
385   fCurrentEventGoodV0s.clear();\r
386   \r
387   fCurrentV0IndexNumber=0;\r
388 }\r
389 \r
390 Double_t AliV0Reader::GetNegativeTrackPhi() const{\r
391   //see header file for documentation\r
392 \r
393   Double_t offset=0;\r
394   if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){\r
395     offset = -2*TMath::Pi();\r
396   }\r
397   return fNegativeTrackLorentzVector->Phi()+offset;\r
398 }\r
399 \r
400 Double_t AliV0Reader::GetPositiveTrackPhi() const{\r
401   //see header file for documentation\r
402 \r
403   Double_t offset=0;\r
404   if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){\r
405     offset = -2*TMath::Pi();\r
406   }\r
407   return fPositiveTrackLorentzVector->Phi()+offset;\r
408 }\r
409 \r
410 Double_t AliV0Reader::GetMotherCandidatePhi() const{\r
411   //see header file for documentation\r
412 \r
413   Double_t offset=0;\r
414   if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){\r
415     offset = -2*TMath::Pi();\r
416   }\r
417   return fMotherCandidateLorentzVector->Phi()+offset;\r
418 }\r
419 \r
420 Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){\r
421   //see header file for documentation\r
422 \r
423   Int_t iResult = 10; // Unknown particle\r
424 \r
425   if(chargeOfTrack==-1){ //negative track\r
426     switch(abs(fNegativeTrackPID)){\r
427     case 11:       //electron\r
428       iResult = 0;\r
429       break;\r
430     case 13:       //muon\r
431       iResult = 1;\r
432       break;\r
433     case 211:      //pion\r
434       iResult = 2;\r
435       break;\r
436     case 321:      //kaon\r
437       iResult = 3;\r
438       break;\r
439     case 2212:     //proton\r
440       iResult = 4;\r
441       break;\r
442     case 22:       //photon\r
443       iResult = 5;\r
444       break;\r
445     case 111:      //pi0\r
446       iResult = 6;\r
447       break;\r
448     case 2112:     //neutron\r
449       iResult = 7;\r
450       break;\r
451     case 311:      //K0\r
452       iResult = 8;\r
453       break;\r
454       \r
455       //Put in here for kSPECIES::kEleCon  ????\r
456     }\r
457   }\r
458   else if(chargeOfTrack==1){ //positive track\r
459     switch(abs(fPositiveTrackPID)){\r
460     case 11:       //electron\r
461       iResult = 0;\r
462       break;\r
463     case 13:       //muon\r
464       iResult = 1;\r
465       break;\r
466     case 211:      //pion\r
467       iResult = 2;\r
468       break;\r
469     case 321:      //kaon\r
470       iResult = 3;\r
471       break;\r
472     case 2212:     //proton\r
473       iResult = 4;\r
474       break;\r
475     case 22:       //photon\r
476       iResult = 5;\r
477       break;\r
478     case 111:      //pi0\r
479       iResult = 6;\r
480       break;\r
481     case 2112:     //neutron\r
482       iResult = 7;\r
483       break;\r
484     case 311:      //K0\r
485       iResult = 8;\r
486       break;\r
487 \r
488       //Put in here for kSPECIES::kEleCon  ????\r
489     }\r
490   }\r
491   else{\r
492     //Wrong parameter.. Print warning\r
493   }\r
494   return iResult;\r
495 }\r