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