b5ecfbae8299f767323278e5e4749fd3470a312a
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / 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     fChi2CutConversion(0.),\r
81     fChi2CutMeson(0.),\r
82     fPIDProbabilityCutNegativeParticle(0),\r
83     fPIDProbabilityCutPositiveParticle(0),\r
84     fXVertexCut(0.),\r
85     fYVertexCut(0.),\r
86     fZVertexCut(0.),\r
87     fNSigmaMass(0.),\r
88     fUseImprovedVertex(kFALSE),\r
89     fCurrentEventGoodV0s(),\r
90     fPreviousEventGoodV0s()\r
91 {\r
92 \r
93 }\r
94 \r
95 \r
96 AliV0Reader::AliV0Reader(const AliV0Reader & original) :\r
97   TObject(original),\r
98   fMCStack(original.fMCStack),\r
99   fMCTruth(original.fMCTruth),\r
100   fChain(original.fChain),\r
101   fESDHandler(original.fESDHandler),\r
102   fESDEvent(original.fESDEvent),\r
103   fHistograms(original.fHistograms),\r
104   fCurrentV0IndexNumber(original.fCurrentV0IndexNumber),\r
105   fCurrentV0(original.fCurrentV0),\r
106   fCurrentNegativeKFParticle(original.fCurrentNegativeKFParticle),\r
107   fCurrentPositiveKFParticle(original.fCurrentPositiveKFParticle),\r
108   fCurrentMotherKFCandidate(original.fCurrentMotherKFCandidate),\r
109   fCurrentNegativeESDTrack(original.fCurrentNegativeESDTrack),\r
110   fCurrentPositiveESDTrack(original.fCurrentPositiveESDTrack),\r
111   fNegativeTrackLorentzVector(original.fNegativeTrackLorentzVector),\r
112   fPositiveTrackLorentzVector(original.fPositiveTrackLorentzVector),\r
113   fMotherCandidateLorentzVector(original.fMotherCandidateLorentzVector),\r
114   fCurrentXValue(original.fCurrentXValue),\r
115   fCurrentYValue(original.fCurrentYValue),\r
116   fCurrentZValue(original.fCurrentZValue),\r
117   fPositiveTrackPID(original.fPositiveTrackPID),\r
118   fNegativeTrackPID(original.fNegativeTrackPID),\r
119   fNegativeMCParticle(original.fNegativeMCParticle),\r
120   fPositiveMCParticle(original.fPositiveMCParticle),\r
121   fMotherMCParticle(original.fMotherMCParticle),\r
122   fMotherCandidateKFMass(original.fMotherCandidateKFMass),\r
123   fMotherCandidateKFWidth(original.fMotherCandidateKFWidth),\r
124   fUseKFParticle(kTRUE),\r
125   fUseESDTrack(kFALSE),\r
126   fDoMC(kFALSE),\r
127   fMaxR(original.fMaxR),\r
128   fEtaCut(original.fEtaCut),\r
129   fPtCut(original.fPtCut),\r
130   fChi2CutConversion(original.fChi2CutConversion),\r
131   fChi2CutMeson(original.fChi2CutMeson),\r
132   fPIDProbabilityCutNegativeParticle(original.fPIDProbabilityCutNegativeParticle),\r
133   fPIDProbabilityCutPositiveParticle(original.fPIDProbabilityCutPositiveParticle),\r
134   fXVertexCut(original.fXVertexCut),\r
135   fYVertexCut(original.fYVertexCut),\r
136   fZVertexCut(original.fZVertexCut),\r
137   fNSigmaMass(original.fNSigmaMass),\r
138   fUseImprovedVertex(original.fUseImprovedVertex),\r
139   fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),\r
140   fPreviousEventGoodV0s(original.fPreviousEventGoodV0s)\r
141 {\r
142 \r
143 }\r
144 \r
145 \r
146 AliV0Reader & AliV0Reader::operator = (const AliV0Reader & /*source*/)\r
147 {\r
148   // assignment operator\r
149   return *this;\r
150 }\r
151 \r
152 void AliV0Reader::Initialize(){\r
153   //see header file for documentation\r
154 \r
155   // Get the input handler from the manager\r
156   fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
157   if(fESDHandler == NULL){\r
158     //print warning here\r
159   }\r
160   \r
161   // Get pointer to esd event from input handler\r
162   fESDEvent = fESDHandler->GetEvent();\r
163   if(fESDEvent == NULL){\r
164     //print warning here\r
165   }\r
166 \r
167   //Get pointer to MCTruth\r
168   fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());\r
169   if(fMCTruth == NULL){\r
170     //print warning here\r
171   }\r
172 \r
173   //Get pointer to the mc stack\r
174   fMCStack = fMCTruth->MCEvent()->Stack();\r
175   if(fMCStack == NULL){\r
176     //print warning here\r
177   }\r
178 \r
179   AliKFParticle::SetField(fESDEvent->GetMagneticField());\r
180 \r
181 }\r
182 \r
183 AliESDv0* AliV0Reader::GetV0(Int_t index){\r
184   //see header file for documentation\r
185 \r
186   fCurrentV0 = fESDEvent->GetV0(index);\r
187   UpdateV0Information();\r
188   return fCurrentV0;\r
189 }\r
190 Bool_t AliV0Reader::CheckForPrimaryVertex(){\r
191   return fESDEvent->GetPrimaryVertex()->GetNContributors()>0;\r
192 }\r
193 \r
194 Bool_t AliV0Reader::NextV0(){\r
195   //see header file for documentation\r
196 \r
197   Bool_t iResult=kFALSE;\r
198   while(fCurrentV0IndexNumber<fESDEvent->GetNumberOfV0s()){\r
199     fCurrentV0 = fESDEvent->GetV0(fCurrentV0IndexNumber);\r
200     \r
201     //checks if on the fly mode is set\r
202     if ( !fCurrentV0->GetOnFlyStatus() ){\r
203       fCurrentV0IndexNumber++;\r
204       if(fHistograms != NULL){\r
205         fHistograms->FillHistogram("V0MassDebugCut1",GetMotherCandidateMass());\r
206       }\r
207       continue;\r
208     }\r
209 \r
210     if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {//checks if we have a vertex\r
211       fCurrentV0IndexNumber++;\r
212       if(fHistograms != NULL){\r
213         fHistograms->FillHistogram("V0MassDebugCut2",GetMotherCandidateMass());\r
214       }\r
215       continue;\r
216     }\r
217 \r
218     if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){\r
219       fCurrentV0IndexNumber++;\r
220       if(fHistograms != NULL){\r
221         fHistograms->FillHistogram("V0MassDebugCut3",GetMotherCandidateMass());\r
222       }\r
223       continue;\r
224     }\r
225 \r
226     fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);\r
227  \r
228     if(GetXYRadius()>fMaxR){ // cuts on distance from collision point\r
229       fCurrentV0IndexNumber++;\r
230       if(fHistograms != NULL){\r
231         fHistograms->FillHistogram("V0MassDebugCut4",GetMotherCandidateMass());\r
232       }\r
233       continue;\r
234     }\r
235 \r
236     UpdateV0Information();\r
237         \r
238     if(fUseKFParticle){\r
239       if(fCurrentMotherKFCandidate->GetNDF()<=0){\r
240         fCurrentV0IndexNumber++;\r
241         if(fHistograms != NULL){\r
242           fHistograms->FillHistogram("V0MassDebugCut5",GetMotherCandidateMass());\r
243         }\r
244         continue;\r
245       }\r
246       Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();\r
247       if(chi2V0 > fChi2CutConversion || chi2V0 <=0){\r
248         fCurrentV0IndexNumber++;\r
249         if(fHistograms != NULL){\r
250           fHistograms->FillHistogram("V0MassDebugCut6",GetMotherCandidateMass());\r
251         }\r
252         continue;\r
253       }\r
254       \r
255       if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut){\r
256         fCurrentV0IndexNumber++;\r
257         if(fHistograms != NULL){\r
258           fHistograms->FillHistogram("V0MassDebugCut7",GetMotherCandidateMass());\r
259         }\r
260         continue;\r
261       }\r
262       \r
263       if(fMotherCandidateLorentzVector->Pt()<fPtCut){\r
264         fCurrentV0IndexNumber++;\r
265         if(fHistograms != NULL){\r
266           fHistograms->FillHistogram("V0MassDebugCut8",GetMotherCandidateMass());\r
267         }\r
268         continue;\r
269       }\r
270     }\r
271     else if(fUseESDTrack){\r
272       //TODO\r
273     }\r
274 \r
275     iResult=kTRUE;//means we have a v0 who survived all the cuts applied\r
276 \r
277     fCurrentV0IndexNumber++;\r
278     \r
279     break;\r
280   }\r
281   return iResult; \r
282 }\r
283 \r
284 void AliV0Reader::UpdateV0Information(){\r
285   //see header file for documentation\r
286   \r
287   if(fCurrentNegativeKFParticle != NULL){\r
288     delete fCurrentNegativeKFParticle;\r
289   }\r
290   fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fNegativeTrackPID);\r
291   \r
292   if(fCurrentPositiveKFParticle != NULL){\r
293     delete fCurrentPositiveKFParticle;\r
294   }\r
295   fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fPositiveTrackPID);\r
296     \r
297   if(fCurrentMotherKFCandidate != NULL){\r
298     delete fCurrentMotherKFCandidate;\r
299   }\r
300   fCurrentMotherKFCandidate = new AliKFParticle(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);\r
301 \r
302   fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());\r
303 \r
304   fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());\r
305 \r
306   if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){\r
307     fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);\r
308   }\r
309 \r
310   if(fUseImprovedVertex == kTRUE){\r
311     AliKFVertex primaryVertexImproved(*GetPrimaryVertex());\r
312     primaryVertexImproved+=*fCurrentMotherKFCandidate;\r
313     fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);\r
314   }\r
315 \r
316   fCurrentMotherKFCandidate->GetMass(fMotherCandidateKFMass,fMotherCandidateKFWidth);\r
317 \r
318 \r
319   if(fNegativeTrackLorentzVector != NULL){\r
320     delete fNegativeTrackLorentzVector;\r
321   }\r
322   if(fUseKFParticle){\r
323     fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeKFParticle->Px(),fCurrentNegativeKFParticle->Py(),fCurrentNegativeKFParticle->Pz());\r
324   }\r
325   else if(fUseESDTrack){\r
326     fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeESDTrack->Px(),fCurrentNegativeESDTrack->Py(),fCurrentNegativeESDTrack->Pz());\r
327   }\r
328 \r
329   if(fPositiveTrackLorentzVector != NULL){\r
330     delete fPositiveTrackLorentzVector;\r
331   }\r
332   if(fUseKFParticle){\r
333     fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveKFParticle->Px(),fCurrentPositiveKFParticle->Py(),fCurrentPositiveKFParticle->Pz());\r
334   }\r
335   else if(fUseESDTrack){\r
336     fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveESDTrack->Px(),fCurrentPositiveESDTrack->Py(),fCurrentPositiveESDTrack->Pz());\r
337   }\r
338 \r
339   if(fMotherCandidateLorentzVector != NULL){\r
340     delete fMotherCandidateLorentzVector;\r
341   }\r
342   if(fUseKFParticle){\r
343     fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);\r
344   }\r
345   else if(fUseESDTrack){\r
346     fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);\r
347   }\r
348 \r
349   if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){\r
350     fMotherCandidateLorentzVector->SetXYZM(fMotherCandidateLorentzVector->Px() ,fMotherCandidateLorentzVector->Py(),fMotherCandidateLorentzVector->Pz(),0.); \r
351   }\r
352     \r
353   if(fDoMC == kTRUE){\r
354     fMotherMCParticle= NULL;\r
355     fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));\r
356     fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));\r
357     if(fPositiveMCParticle->GetMother(0)>-1){\r
358       fMotherMCParticle = fMCStack->Particle(fPositiveMCParticle->GetMother(0));\r
359     }\r
360   }\r
361   fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate);\r
362 }\r
363 \r
364 Bool_t AliV0Reader::HasSameMCMother(){\r
365   //see header file for documentation\r
366 \r
367   Bool_t iResult = kFALSE;\r
368   if(fDoMC == kTRUE){\r
369     if(fNegativeMCParticle != NULL && fPositiveMCParticle != NULL){\r
370       if(fNegativeMCParticle->GetMother(0) == fPositiveMCParticle->GetMother(0))\r
371         if(fMotherMCParticle){\r
372           iResult = kTRUE;\r
373         }\r
374     }\r
375   }\r
376   return iResult;\r
377 }\r
378 \r
379 Bool_t AliV0Reader::CheckPIDProbability(Double_t negProbCut, Double_t posProbCut){\r
380   //see header file for documentation\r
381 \r
382   Bool_t iResult=kFALSE;\r
383 \r
384   Double_t *posProbArray = new Double_t[10];\r
385   Double_t *negProbArray = new Double_t[10];\r
386   AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());\r
387   AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());\r
388   \r
389   negTrack->GetTPCpid(negProbArray);\r
390   posTrack->GetTPCpid(posProbArray);\r
391 \r
392   if(negProbArray!=NULL && posProbArray!=NULL){\r
393     if(negProbArray[GetSpeciesIndex(-1)]>=negProbCut && posProbArray[GetSpeciesIndex(1)]>=posProbCut){\r
394       iResult=kTRUE;\r
395     }\r
396   }\r
397   delete [] posProbArray;\r
398   delete [] negProbArray;\r
399   return iResult;\r
400 }\r
401 \r
402 void AliV0Reader::GetPIDProbability(Double_t &negPIDProb,Double_t & posPIDProb){\r
403 \r
404   Double_t *posProbArray = new Double_t[10];\r
405   Double_t *negProbArray = new Double_t[10];\r
406   AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());\r
407   AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());\r
408   \r
409   negTrack->GetTPCpid(negProbArray);\r
410   posTrack->GetTPCpid(posProbArray);\r
411 \r
412   if(negProbArray!=NULL && posProbArray!=NULL){\r
413     negPIDProb = negProbArray[GetSpeciesIndex(-1)];\r
414     posPIDProb = posProbArray[GetSpeciesIndex(1)];\r
415   }\r
416   delete [] posProbArray;\r
417   delete [] negProbArray;\r
418 }\r
419 \r
420 void AliV0Reader::UpdateEventByEventData(){\r
421   //see header file for documentation\r
422 \r
423   if(fCurrentEventGoodV0s.size() >0 ){\r
424     fPreviousEventGoodV0s.clear();\r
425     fPreviousEventGoodV0s = fCurrentEventGoodV0s;\r
426   }\r
427   fCurrentEventGoodV0s.clear();\r
428   \r
429   fCurrentV0IndexNumber=0;\r
430 }\r
431 \r
432 Double_t AliV0Reader::GetNegativeTrackPhi() const{\r
433   //see header file for documentation\r
434 \r
435   Double_t offset=0;\r
436   if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){\r
437     offset = -2*TMath::Pi();\r
438   }\r
439   return fNegativeTrackLorentzVector->Phi()+offset;\r
440 }\r
441 \r
442 Double_t AliV0Reader::GetPositiveTrackPhi() const{\r
443   //see header file for documentation\r
444 \r
445   Double_t offset=0;\r
446   if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){\r
447     offset = -2*TMath::Pi();\r
448   }\r
449   return fPositiveTrackLorentzVector->Phi()+offset;\r
450 }\r
451 \r
452 Double_t AliV0Reader::GetMotherCandidatePhi() const{\r
453   //see header file for documentation\r
454 \r
455   Double_t offset=0;\r
456   if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){\r
457     offset = -2*TMath::Pi();\r
458   }\r
459   return fMotherCandidateLorentzVector->Phi()+offset;\r
460 }\r
461 \r
462 Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){\r
463   //see header file for documentation\r
464 \r
465   Int_t iResult = 10; // Unknown particle\r
466 \r
467   if(chargeOfTrack==-1){ //negative track\r
468     switch(abs(fNegativeTrackPID)){\r
469     case 11:       //electron\r
470       iResult = 0;\r
471       break;\r
472     case 13:       //muon\r
473       iResult = 1;\r
474       break;\r
475     case 211:      //pion\r
476       iResult = 2;\r
477       break;\r
478     case 321:      //kaon\r
479       iResult = 3;\r
480       break;\r
481     case 2212:     //proton\r
482       iResult = 4;\r
483       break;\r
484     case 22:       //photon\r
485       iResult = 5;\r
486       break;\r
487     case 111:      //pi0\r
488       iResult = 6;\r
489       break;\r
490     case 2112:     //neutron\r
491       iResult = 7;\r
492       break;\r
493     case 311:      //K0\r
494       iResult = 8;\r
495       break;\r
496       \r
497       //Put in here for kSPECIES::kEleCon  ????\r
498     }\r
499   }\r
500   else if(chargeOfTrack==1){ //positive track\r
501     switch(abs(fPositiveTrackPID)){\r
502     case 11:       //electron\r
503       iResult = 0;\r
504       break;\r
505     case 13:       //muon\r
506       iResult = 1;\r
507       break;\r
508     case 211:      //pion\r
509       iResult = 2;\r
510       break;\r
511     case 321:      //kaon\r
512       iResult = 3;\r
513       break;\r
514     case 2212:     //proton\r
515       iResult = 4;\r
516       break;\r
517     case 22:       //photon\r
518       iResult = 5;\r
519       break;\r
520     case 111:      //pi0\r
521       iResult = 6;\r
522       break;\r
523     case 2112:     //neutron\r
524       iResult = 7;\r
525       break;\r
526     case 311:      //K0\r
527       iResult = 8;\r
528       break;\r
529 \r
530       //Put in here for kSPECIES::kEleCon  ????\r
531     }\r
532   }\r
533   else{\r
534     //Wrong parameter.. Print warning\r
535   }\r
536   return iResult;\r
537 }\r