Added new background scheeme, did some cleanup. added new bethe block parameters...
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / 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 ////////////////////////////////////////////////
17 //--------------------------------------------- 
18 // Class used to do analysis on conversion pairs
19 //---------------------------------------------
20 ////////////////////////////////////////////////
21
22 // --- ROOT system ---
23 #include <TMath.h>
24
25 //---- ANALYSIS system ----
26 #include "AliV0Reader.h"
27 #include "AliAnalysisManager.h"
28 #include "AliESDInputHandler.h"
29 #include "AliESDtrack.h"
30 #include "AliMCEvent.h"
31 #include "AliKFVertex.h"
32
33 #include "AliStack.h"
34 #include "AliMCEventHandler.h"
35 #include "AliTPCpidESD.h"
36 #include "AliGammaConversionBGHandler.h"
37
38 class iostream;
39 class AliESDv0;
40 class TFormula;
41
42 using namespace std;
43
44 ClassImp(AliV0Reader)
45
46
47
48 AliV0Reader::AliV0Reader() :
49   TObject(),
50   fMCStack(NULL),
51   fMCTruth(NULL),
52   fMCEvent(NULL),    // for CF
53   fChain(NULL),
54   fESDHandler(NULL),
55   fESDEvent(NULL),
56   fCFManager(NULL),
57   fTPCpid(NULL),
58   fHistograms(NULL),
59   fCurrentV0IndexNumber(0),
60   fCurrentV0(NULL),
61   fCurrentNegativeKFParticle(NULL),
62   fCurrentPositiveKFParticle(NULL),
63   fCurrentMotherKFCandidate(NULL),
64   fCurrentNegativeESDTrack(NULL),
65   fCurrentPositiveESDTrack(NULL),
66   fNegativeTrackLorentzVector(NULL),
67   fPositiveTrackLorentzVector(NULL),
68   fMotherCandidateLorentzVector(NULL),
69   fCurrentXValue(0),
70   fCurrentYValue(0),
71   fCurrentZValue(0),
72   fPositiveTrackPID(0),
73   fNegativeTrackPID(0),
74   fNegativeMCParticle(NULL),
75   fPositiveMCParticle(NULL),
76   fMotherMCParticle(NULL),
77   fMotherCandidateKFMass(0),
78   fMotherCandidateKFWidth(0),
79   fUseKFParticle(kTRUE),
80   fUseESDTrack(kFALSE),
81   fDoMC(kFALSE),
82   fMaxR(10000),// 100 meter(outside of ALICE)
83   fEtaCut(0.),
84   fPtCut(0.),
85   fMaxZ(0.),
86   fLineCutZRSlope(0.),
87   fLineCutZValue(0.),
88   fChi2CutConversion(0.),
89   fChi2CutMeson(0.),
90   fPIDProbabilityCutNegativeParticle(0),
91   fPIDProbabilityCutPositiveParticle(0),
92   fDodEdxSigmaCut(kFALSE),
93   fPIDnSigmaAboveElectronLine(100),
94   fPIDnSigmaBelowElectronLine(-100),
95   fPIDnSigmaAbovePionLine(-100), 
96   fPIDMinPnSigmaAbovePionLine(100), 
97   fXVertexCut(0.),
98   fYVertexCut(0.),
99   fZVertexCut(0.),
100   fNSigmaMass(0.),
101   fUseImprovedVertex(kFALSE),
102   fUseOwnXYZCalculation(kFALSE),
103   fDoCF(kFALSE),
104   fUseOnFlyV0Finder(kTRUE),
105   fUpdateV0AlreadyCalled(kFALSE),
106   fCurrentEventGoodV0s(NULL),
107 //  fPreviousEventGoodV0s(),
108   fBGEventHandler(NULL),
109   fBGEventInitialized(kFALSE)
110 {
111   fTPCpid = new AliTPCpidESD;   
112 }
113
114
115 AliV0Reader::AliV0Reader(const AliV0Reader & original) :
116   TObject(original),
117   fMCStack(original.fMCStack),
118   fMCTruth(original.fMCTruth),
119   fMCEvent(original.fMCEvent),  // for CF
120   fChain(original.fChain),
121   fESDHandler(original.fESDHandler),
122   fESDEvent(original.fESDEvent),
123   fCFManager(original.fCFManager),
124   fTPCpid(original.fTPCpid),
125   fHistograms(original.fHistograms),
126   fCurrentV0IndexNumber(original.fCurrentV0IndexNumber),
127   fCurrentV0(original.fCurrentV0),
128   fCurrentNegativeKFParticle(original.fCurrentNegativeKFParticle),
129   fCurrentPositiveKFParticle(original.fCurrentPositiveKFParticle),
130   fCurrentMotherKFCandidate(original.fCurrentMotherKFCandidate),
131   fCurrentNegativeESDTrack(original.fCurrentNegativeESDTrack),
132   fCurrentPositiveESDTrack(original.fCurrentPositiveESDTrack),
133   fNegativeTrackLorentzVector(original.fNegativeTrackLorentzVector),
134   fPositiveTrackLorentzVector(original.fPositiveTrackLorentzVector),
135   fMotherCandidateLorentzVector(original.fMotherCandidateLorentzVector),
136   fCurrentXValue(original.fCurrentXValue),
137   fCurrentYValue(original.fCurrentYValue),
138   fCurrentZValue(original.fCurrentZValue),
139   fPositiveTrackPID(original.fPositiveTrackPID),
140   fNegativeTrackPID(original.fNegativeTrackPID),
141   fNegativeMCParticle(original.fNegativeMCParticle),
142   fPositiveMCParticle(original.fPositiveMCParticle),
143   fMotherMCParticle(original.fMotherMCParticle),
144   fMotherCandidateKFMass(original.fMotherCandidateKFMass),
145   fMotherCandidateKFWidth(original.fMotherCandidateKFWidth),
146   fUseKFParticle(kTRUE),
147   fUseESDTrack(kFALSE),
148   fDoMC(kFALSE),
149   fMaxR(original.fMaxR),
150   fEtaCut(original.fEtaCut),
151   fPtCut(original.fPtCut),
152   fMaxZ(original.fMaxZ),
153   fLineCutZRSlope(original.fLineCutZRSlope),
154   fLineCutZValue(original.fLineCutZValue),
155   fChi2CutConversion(original.fChi2CutConversion),
156   fChi2CutMeson(original.fChi2CutMeson),
157   fPIDProbabilityCutNegativeParticle(original.fPIDProbabilityCutNegativeParticle),
158   fPIDProbabilityCutPositiveParticle(original.fPIDProbabilityCutPositiveParticle),
159   fDodEdxSigmaCut(original.fDodEdxSigmaCut),
160   fPIDnSigmaAboveElectronLine(original.fPIDnSigmaAboveElectronLine),
161   fPIDnSigmaBelowElectronLine(original.fPIDnSigmaBelowElectronLine),
162   fPIDnSigmaAbovePionLine(original.fPIDnSigmaAbovePionLine), 
163   fPIDMinPnSigmaAbovePionLine(original.fPIDMinPnSigmaAbovePionLine), 
164   fXVertexCut(original.fXVertexCut),
165   fYVertexCut(original.fYVertexCut),
166   fZVertexCut(original.fZVertexCut),
167   fNSigmaMass(original.fNSigmaMass),
168   fUseImprovedVertex(original.fUseImprovedVertex),
169   fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
170   fDoCF(original.fDoCF),
171   fUseOnFlyV0Finder(original.fUseOnFlyV0Finder),
172   fUpdateV0AlreadyCalled(original.fUpdateV0AlreadyCalled),
173   fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),
174   //  fPreviousEventGoodV0s(original.fPreviousEventGoodV0s),
175   fBGEventHandler(original.fBGEventHandler),
176   fBGEventInitialized(original.fBGEventInitialized)
177 {
178         
179 }
180
181
182 AliV0Reader & AliV0Reader::operator = (const AliV0Reader & /*source*/)
183 {
184   // assignment operator
185   return *this;
186 }
187 AliV0Reader::~AliV0Reader()
188 {
189   if(fTPCpid){
190     delete fTPCpid;
191   }
192 }
193
194 void AliV0Reader::Initialize(){
195   //see header file for documentation
196
197   fUpdateV0AlreadyCalled = kFALSE;      
198   // Get the input handler from the manager
199   fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
200   if(fESDHandler == NULL){
201     //print warning here
202   }
203         
204   // Get pointer to esd event from input handler
205   fESDEvent = fESDHandler->GetEvent();
206   if(fESDEvent == NULL){
207     //print warning here
208   }
209         
210   //Get pointer to MCTruth
211   fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
212   if(fMCTruth == NULL){
213     //print warning here
214     fDoMC = kFALSE;
215   }
216         
217   //Get pointer to the mc stack
218   if(fMCTruth){
219     fMCStack = fMCTruth->MCEvent()->Stack();
220     if(fMCStack == NULL){
221       //print warning here
222     }
223     // Better parameters for MonteCarlo from A. Kalweit 2010/01/8
224     fTPCpid->SetBetheBlochParameters( 2.15898e+00/50.,
225                                       1.75295e+01,
226                                       3.40030e-09,
227                                       1.96178e+00,
228                                       3.91720e+00);
229   }
230   else{
231     // Better parameters for data from A. Kalweit 2010/01/8
232     fTPCpid->SetBetheBlochParameters(0.0283086,
233                                      2.63394e+01,
234                                      5.04114e-11,
235                                      2.12543e+00,
236                                      4.88663e+00);
237   }
238         
239   // for CF
240   //Get pointer to the mc event
241   if(fDoCF && fDoMC){
242     fMCEvent = fMCTruth->MCEvent();
243     if(fMCEvent == NULL){
244       //print warning here
245       fDoCF = kFALSE;
246     }   
247   }
248
249
250
251         
252   AliKFParticle::SetField(fESDEvent->GetMagneticField());
253
254   //  fCurrentEventGoodV0s = new TClonesArray("TClonesArray", 0);
255   fCurrentEventGoodV0s = new TClonesArray("AliKFParticle", 0);
256
257   if(fBGEventInitialized == kFALSE){
258     Double_t *zBinLimitsArray = new Double_t[8];//{-7,-5,-3,-1,1,3,5,7};
259     zBinLimitsArray[0] = -7;
260     zBinLimitsArray[1] = -5;
261     zBinLimitsArray[2] = -3;
262     zBinLimitsArray[3] = -1;
263     zBinLimitsArray[4] = 1;
264     zBinLimitsArray[5] = 3;
265     zBinLimitsArray[6] = 5;
266     zBinLimitsArray[7] = 7;
267     
268     Double_t *multiplicityBinLimitsArray= new Double_t[4];//={0,10,20,500};
269     multiplicityBinLimitsArray[0] = 0;
270     multiplicityBinLimitsArray[1] = 10;
271     multiplicityBinLimitsArray[2] = 20;
272     multiplicityBinLimitsArray[3] = 500;
273     
274     
275     fBGEventHandler = new AliGammaConversionBGHandler(8,4,10);
276     
277     fBGEventHandler->Initialize(zBinLimitsArray, multiplicityBinLimitsArray);
278     fBGEventInitialized = kTRUE;
279   }
280 }
281
282 AliESDv0* AliV0Reader::GetV0(Int_t index){
283   //see header file for documentation
284   fCurrentV0 = fESDEvent->GetV0(index);
285   UpdateV0Information();
286   return fCurrentV0;
287 }
288
289 Bool_t AliV0Reader::CheckForPrimaryVertex(){
290   return fESDEvent->GetPrimaryVertex()->GetNContributors()>0;
291 }
292
293
294 Bool_t AliV0Reader::CheckV0FinderStatus(Int_t index){
295
296   if(fUseOnFlyV0Finder){
297     if(!GetV0(index)->GetOnFlyStatus()){
298       return kFALSE;
299     }
300   }
301   if(!fUseOnFlyV0Finder){
302     if(!GetV0(index)->GetOnFlyStatus()){
303       return kFALSE;
304     }
305   }
306   return kTRUE;
307 }
308
309
310
311 Bool_t AliV0Reader::NextV0(){
312   //see header file for documentation
313         
314   Bool_t iResult=kFALSE;
315   while(fCurrentV0IndexNumber<fESDEvent->GetNumberOfV0s()){
316     fCurrentV0 = fESDEvent->GetV0(fCurrentV0IndexNumber);
317                 
318     // moved it up here so that the correction framework can access pt and eta information
319     if(UpdateV0Information() == kFALSE){
320       fCurrentV0IndexNumber++;
321       continue;
322     }
323  
324     Double_t containerInput[3];
325     if(fDoCF){
326       containerInput[0] = GetMotherCandidatePt();
327       containerInput[1] = GetMotherCandidateEta();
328       containerInput[2] = GetMotherCandidateMass();
329     }
330     
331    
332     //checks if on the fly mode is set
333     if ( !fCurrentV0->GetOnFlyStatus() ){
334       if(fHistograms != NULL){
335         fHistograms->FillHistogram("ESD_CutGetOnFly_InvMass",GetMotherCandidateMass());
336       }
337       fCurrentV0IndexNumber++;
338       continue;
339     }
340     if(fDoCF){
341       fCFManager->GetParticleContainer()->Fill(containerInput,kStepGetOnFly);           // for CF       
342     }
343     
344     //checks if we have a prim vertex
345     if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) { 
346       if(fHistograms != NULL){
347         fHistograms->FillHistogram("ESD_CutNContributors_InvMass",GetMotherCandidateMass());
348       }
349       fCurrentV0IndexNumber++;
350       continue;
351     }
352     if(fDoCF){
353       fCFManager->GetParticleContainer()->Fill(containerInput,kStepNContributors);              // for CF       
354     }
355                 
356     //Check the pid probability
357     if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){
358       if(fHistograms != NULL){
359         fHistograms->FillHistogram("ESD_CutPIDProb_InvMass",GetMotherCandidateMass());
360       }
361       fCurrentV0IndexNumber++;
362       continue;
363     }
364     if(fDoCF){
365       fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCPID);                     // for CF
366     }
367                 
368     if(GetXYRadius()>fMaxR){ // cuts on distance from collision point
369       if(fHistograms != NULL){
370         fHistograms->FillHistogram("ESD_CutR_InvMass",GetMotherCandidateMass());
371       }
372       fCurrentV0IndexNumber++;
373       continue;
374     }   
375     if(fDoCF){
376       fCFManager->GetParticleContainer()->Fill(containerInput,kStepR);                  // for CF
377     }
378                 
379                 
380     if((TMath::Abs(fCurrentZValue)*fLineCutZRSlope)-fLineCutZValue > GetXYRadius() ){ // cuts out regions where we do not reconstruct
381       if(fHistograms != NULL){
382         fHistograms->FillHistogram("ESD_CutLine_InvMass",GetMotherCandidateMass());
383       }
384       fCurrentV0IndexNumber++;
385       continue;
386     }
387     if(fDoCF){
388       fCFManager->GetParticleContainer()->Fill(containerInput,kStepLine);                       // for CF
389     }
390                 
391     if(TMath::Abs(fCurrentZValue) > fMaxZ ){ // cuts out regions where we do not reconstruct
392       if(fHistograms != NULL){
393         fHistograms->FillHistogram("ESD_CutZ_InvMass",GetMotherCandidateMass());
394       }
395       fCurrentV0IndexNumber++;
396       continue;
397     }
398     if(fDoCF){
399       fCFManager->GetParticleContainer()->Fill(containerInput,kStepZ);          // for CF       
400     }
401                 
402     /* Moved further up so corr framework can work
403        if(UpdateV0Information() == kFALSE){
404        fCurrentV0IndexNumber++;
405        continue;
406        }
407     */
408
409                 
410     if(fUseKFParticle){
411       if(fCurrentMotherKFCandidate->GetNDF()<=0){
412         if(fHistograms != NULL){
413           fHistograms->FillHistogram("ESD_CutNDF_InvMass",GetMotherCandidateMass());
414         }
415         fCurrentV0IndexNumber++;
416         continue;
417       }
418       if(fDoCF){
419         fCFManager->GetParticleContainer()->Fill(containerInput,kStepNDF);              // for CF       
420       }
421                         
422       Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();
423       if(chi2V0 > fChi2CutConversion || chi2V0 <=0){
424         if(fHistograms != NULL){
425           fHistograms->FillHistogram("ESD_CutChi2_InvMass",GetMotherCandidateMass());
426         }
427         fCurrentV0IndexNumber++;
428         continue;
429       }
430       if(fDoCF){
431         fCFManager->GetParticleContainer()->Fill(containerInput,kStepChi2);                     // for CF
432       }
433                         
434       if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut){
435         if(fHistograms != NULL){
436           fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
437         }
438         fCurrentV0IndexNumber++;
439         continue;
440       }
441       if(fDoCF){
442         fCFManager->GetParticleContainer()->Fill(containerInput,kStepEta);                      // for CF
443       }
444                         
445       if(fMotherCandidateLorentzVector->Pt()<fPtCut){
446         if(fHistograms != NULL){
447           fHistograms->FillHistogram("ESD_CutPt_InvMass",GetMotherCandidateMass());
448         }
449         fCurrentV0IndexNumber++;
450         continue;
451       }
452       if(fDoCF){
453         fCFManager->GetParticleContainer()->Fill(containerInput,kStepPt);                       // for CF
454       }
455                         
456     }
457     else if(fUseESDTrack){
458       //TODO
459     }
460
461     if(fHistograms != NULL){
462       fHistograms->FillHistogram("ESD_GoodV0s_InvMass",GetMotherCandidateMass());
463     }
464
465     //    fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate);
466
467     new((*fCurrentEventGoodV0s)[fCurrentEventGoodV0s->GetEntriesFast()])  AliKFParticle(*fCurrentMotherKFCandidate);
468
469     iResult=kTRUE;//means we have a v0 who survived all the cuts applied
470                 
471     fCurrentV0IndexNumber++;
472                 
473     break;
474   }
475   return iResult; 
476 }
477
478 Bool_t AliV0Reader::UpdateV0Information(){
479   //see header file for documentation
480         
481   Bool_t iResult=kTRUE;                                         // for taking out not refitted, kinks and like sign tracks 
482         
483   Bool_t switchTracks = kFALSE;
484         
485   fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
486   fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
487         
488   if(fCurrentNegativeESDTrack->GetSign() == fCurrentPositiveESDTrack->GetSign()){             // avoid like sign
489     iResult=kFALSE;
490     if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
491       fHistograms->FillHistogram("ESD_CutLikeSign_InvMass",GetMotherCandidateMass());
492     }
493   }
494         
495   if(fCurrentPositiveESDTrack->GetSign() == -1 && fCurrentNegativeESDTrack->GetSign() == 1){  // switch wrong signed tracks
496     fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
497     fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
498     switchTracks = kTRUE;
499   }
500         
501   if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTPCrefit) || 
502       !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ){
503     //  if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kITSrefit) || 
504     //      !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kITSrefit) ){
505     iResult=kFALSE;
506     if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
507       fHistograms->FillHistogram("ESD_CutRefit_InvMass",GetMotherCandidateMass());
508     }
509   }
510         
511   if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 || 
512       fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {                  
513                 
514     iResult=kFALSE;
515     if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
516       fHistograms->FillHistogram("ESD_CutKink_InvMass",GetMotherCandidateMass());
517     }
518   }
519
520   if(fDodEdxSigmaCut == kTRUE){
521
522     if( fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
523         fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ||
524         fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
525         fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ){
526       iResult=kFALSE;
527       if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
528         fHistograms->FillHistogram("ESD_CutdEdxSigmaElectronLine_InvMass",GetMotherCandidateMass());
529       }
530     }
531     if( fCurrentPositiveESDTrack->P()>fPIDMinPnSigmaAbovePionLine){
532       if(fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
533          fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
534          fTPCpid->GetNumberOfSigmas(fCurrentPositiveESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
535         iResult=kFALSE;
536         if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
537           fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
538         }
539       }
540     }
541
542     if( fCurrentNegativeESDTrack->P()>fPIDMinPnSigmaAbovePionLine){
543       if(fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
544          fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
545          fTPCpid->GetNumberOfSigmas(fCurrentNegativeESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
546         iResult=kFALSE;
547         if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE){
548           fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
549         }
550       }
551     }
552   }
553         
554   if(fCurrentNegativeKFParticle != NULL){
555     delete fCurrentNegativeKFParticle;
556   }
557   if(switchTracks == kFALSE){
558     fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fNegativeTrackPID);
559   }
560   else{
561     fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fNegativeTrackPID);
562   }
563         
564   if(fCurrentPositiveKFParticle != NULL){
565     delete fCurrentPositiveKFParticle;
566   }
567   if(switchTracks == kFALSE){
568     fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fPositiveTrackPID);
569   }
570   else{
571     fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fPositiveTrackPID);
572   }
573     
574   if(fCurrentMotherKFCandidate != NULL){
575     delete fCurrentMotherKFCandidate;
576   }
577   fCurrentMotherKFCandidate = new AliKFParticle(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
578         
579         
580   if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
581     fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);
582   }
583         
584   if(fUseImprovedVertex == kTRUE){
585     AliKFVertex primaryVertexImproved(*GetPrimaryVertex());
586     primaryVertexImproved+=*fCurrentMotherKFCandidate;
587     fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);
588   }
589         
590   fCurrentMotherKFCandidate->GetMass(fMotherCandidateKFMass,fMotherCandidateKFWidth);
591         
592         
593   if(fNegativeTrackLorentzVector != NULL){
594     delete fNegativeTrackLorentzVector;
595   }
596   if(fUseKFParticle){
597     fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeKFParticle->Px(),fCurrentNegativeKFParticle->Py(),fCurrentNegativeKFParticle->Pz());
598   }
599   else if(fUseESDTrack){
600     fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeESDTrack->Px(),fCurrentNegativeESDTrack->Py(),fCurrentNegativeESDTrack->Pz());
601   }
602         
603   if(fPositiveTrackLorentzVector != NULL){
604     delete fPositiveTrackLorentzVector;
605   }
606   if(fUseKFParticle){
607     fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveKFParticle->Px(),fCurrentPositiveKFParticle->Py(),fCurrentPositiveKFParticle->Pz());
608   }
609   else if(fUseESDTrack){
610     fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveESDTrack->Px(),fCurrentPositiveESDTrack->Py(),fCurrentPositiveESDTrack->Pz());
611   }
612         
613   if(fMotherCandidateLorentzVector != NULL){
614     delete fMotherCandidateLorentzVector;
615   }
616   if(fUseKFParticle){
617     fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
618   }
619   else if(fUseESDTrack){
620     fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
621   }
622         
623   if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
624     fMotherCandidateLorentzVector->SetXYZM(fMotherCandidateLorentzVector->Px() ,fMotherCandidateLorentzVector->Py(),fMotherCandidateLorentzVector->Pz(),0.); 
625   }
626     
627         
628   if(fDoMC == kTRUE){
629     fMotherMCParticle= NULL;
630     fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
631     fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
632     if(fPositiveMCParticle->GetMother(0)>-1){
633       fMotherMCParticle = fMCStack->Particle(fPositiveMCParticle->GetMother(0));
634     }
635   }
636         
637
638   // for CF
639   Double_t containerInput[3];
640   if(fDoCF){
641     containerInput[0] = GetMotherCandidatePt();
642     containerInput[1] = GetMotherCandidateEta();
643     containerInput[2] = GetMotherCandidateMass();
644     
645     fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign);             // for CF       
646     fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit);             // for CF       
647     fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks);                // for CF       
648   }
649   
650
651   if(fUseOwnXYZCalculation == kFALSE){
652     fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);
653   }
654   else{
655     Double_t convpos[2];
656     convpos[0]=0;
657     convpos[1]=0;
658     GetConvPosXY(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField(),convpos);
659     fCurrentXValue = convpos[0];
660     fCurrentYValue = convpos[1];
661     fCurrentZValue = GetConvPosZ(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField());
662   }
663
664
665   fUpdateV0AlreadyCalled = kTRUE;
666
667   return iResult;
668 }
669
670
671
672 Bool_t AliV0Reader::HasSameMCMother(){
673   //see header file for documentation
674         
675   Bool_t iResult = kFALSE;
676   if(fDoMC == kTRUE){
677     if(fNegativeMCParticle != NULL && fPositiveMCParticle != NULL){
678       if(fNegativeMCParticle->GetMother(0) == fPositiveMCParticle->GetMother(0))
679         if(fMotherMCParticle){
680           iResult = kTRUE;
681         }
682     }
683   }
684   return iResult;
685 }
686
687 Bool_t AliV0Reader::CheckPIDProbability(Double_t negProbCut, Double_t posProbCut){
688   //see header file for documentation
689         
690   Bool_t iResult=kFALSE;
691         
692   Double_t *posProbArray = new Double_t[10];
693   Double_t *negProbArray = new Double_t[10];
694   AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());
695   AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());
696         
697   negTrack->GetTPCpid(negProbArray);
698   posTrack->GetTPCpid(posProbArray);
699         
700   //  if(negProbArray != NULL && posProbArray != NULL){ // this is not allowed anymore for some reason(RC19)
701   if(negProbArray && posProbArray){
702     if(negProbArray[GetSpeciesIndex(-1)]>=negProbCut && posProbArray[GetSpeciesIndex(1)]>=posProbCut){
703       iResult=kTRUE;
704     }
705   }
706   delete [] posProbArray;
707   delete [] negProbArray;
708   return iResult;
709 }
710
711 void AliV0Reader::GetPIDProbability(Double_t &negPIDProb,Double_t & posPIDProb){
712   // see header file for documentation
713
714   Double_t *posProbArray = new Double_t[10];
715   Double_t *negProbArray = new Double_t[10];
716   AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());
717   AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());
718         
719   negTrack->GetTPCpid(negProbArray);
720   posTrack->GetTPCpid(posProbArray);
721         
722   //  if(negProbArray!=NULL && posProbArray!=NULL){ // this is not allowed anymore for some reason(RC19)
723   if(negProbArray && posProbArray){
724     negPIDProb = negProbArray[GetSpeciesIndex(-1)];
725     posPIDProb = posProbArray[GetSpeciesIndex(1)];
726   }
727   delete [] posProbArray;
728   delete [] negProbArray;
729 }
730
731 void AliV0Reader::UpdateEventByEventData(){
732   //see header file for documentation
733   if(fCurrentEventGoodV0s->GetEntriesFast() >0 ){
734     fBGEventHandler->AddEvent(fCurrentEventGoodV0s,fESDEvent->GetPrimaryVertex()->GetZ(),fESDEvent->GetNumberOfTracks());
735   }
736   fCurrentEventGoodV0s->Delete();
737   fCurrentV0IndexNumber=0;
738 }
739
740
741 Double_t AliV0Reader::GetNegativeTrackPhi() const{
742   //see header file for documentation
743         
744   Double_t offset=0;
745   if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){
746     offset = -2*TMath::Pi();
747   }
748   return fNegativeTrackLorentzVector->Phi()+offset;
749 }
750
751 Double_t AliV0Reader::GetPositiveTrackPhi() const{
752   //see header file for documentation
753         
754   Double_t offset=0;
755   if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){
756     offset = -2*TMath::Pi();
757   }
758   return fPositiveTrackLorentzVector->Phi()+offset;
759 }
760
761 Double_t AliV0Reader::GetMotherCandidatePhi() const{
762   //see header file for documentation
763         
764   Double_t offset=0;
765   if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){
766     offset = -2*TMath::Pi();
767   }
768   return fMotherCandidateLorentzVector->Phi()+offset;
769 }
770
771
772 Double_t AliV0Reader::GetMotherCandidateRapidity() const{
773   //see header file for documentation
774         
775   Double_t rapidity=0;
776   if(fMotherCandidateLorentzVector->Energy() - fMotherCandidateLorentzVector->Pz() == 0 || fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz() == 0) rapidity=0;
777   else rapidity = 0.5*(TMath::Log((fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz()) / (fMotherCandidateLorentzVector->Energy()-fMotherCandidateLorentzVector->Pz())));
778   return rapidity;
779         
780 }
781
782
783
784
785
786 Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){
787   //see header file for documentation
788         
789   Int_t iResult = 10; // Unknown particle
790         
791   if(chargeOfTrack==-1){ //negative track
792     switch(abs(fNegativeTrackPID)){
793     case 11:       //electron
794       iResult = 0;
795       break;
796     case 13:       //muon
797       iResult = 1;
798       break;
799     case 211:      //pion
800       iResult = 2;
801       break;
802     case 321:      //kaon
803       iResult = 3;
804       break;
805     case 2212:     //proton
806       iResult = 4;
807       break;
808     case 22:       //photon
809       iResult = 5;
810       break;
811     case 111:      //pi0
812       iResult = 6;
813       break;
814     case 2112:     //neutron
815       iResult = 7;
816       break;
817     case 311:      //K0
818       iResult = 8;
819       break;
820                                 
821       //Put in here for kSPECIES::kEleCon  ????
822     }
823   }
824   else if(chargeOfTrack==1){ //positive track
825     switch(abs(fPositiveTrackPID)){
826     case 11:       //electron
827       iResult = 0;
828       break;
829     case 13:       //muon
830       iResult = 1;
831       break;
832     case 211:      //pion
833       iResult = 2;
834       break;
835     case 321:      //kaon
836       iResult = 3;
837       break;
838     case 2212:     //proton
839       iResult = 4;
840       break;
841     case 22:       //photon
842       iResult = 5;
843       break;
844     case 111:      //pi0
845       iResult = 6;
846       break;
847     case 2112:     //neutron
848       iResult = 7;
849       break;
850     case 311:      //K0
851       iResult = 8;
852       break;
853                                 
854       //Put in here for kSPECIES::kEleCon  ????
855     }
856   }
857   else{
858     //Wrong parameter.. Print warning
859   }
860   return iResult;
861 }
862
863 Bool_t AliV0Reader::GetHelixCenter(AliESDtrack* track, Double_t b,Int_t charge, Double_t center[2]){
864   // see header file for documentation
865   
866   Double_t pi = 3.14159265358979323846;
867   
868   Double_t  helix[6];
869   track->GetHelixParameters(helix,b);
870   
871   Double_t xpos =  helix[5];
872   Double_t ypos =  helix[0];
873   Double_t radius = TMath::Abs(1./helix[4]);
874   Double_t phi = helix[2];
875
876   if(phi < 0){
877     phi = phi + 2*pi;
878   }
879
880   phi -= pi/2.;
881   Double_t xpoint =  radius * TMath::Cos(phi);
882   Double_t ypoint =  radius * TMath::Sin(phi);
883
884   if(charge > 0){
885     xpoint = - xpoint;
886     ypoint = - ypoint;
887   }
888
889   if(charge < 0){
890     xpoint =  xpoint;
891     ypoint =  ypoint;
892   }
893   center[0] =  xpos + xpoint;
894   center[1] =  ypos + ypoint;
895
896   return 1;
897 }
898
899 Bool_t AliV0Reader::GetConvPosXY(AliESDtrack* ptrack, AliESDtrack* ntrack, Double_t b, Double_t convpos[2]){
900   //see header file for documentation
901
902   Double_t helixcenterpos[2];
903   GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
904
905   Double_t helixcenterneg[2];
906   GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
907
908   Double_t  poshelix[6];
909   ptrack->GetHelixParameters(poshelix,b);
910   Double_t posradius = TMath::Abs(1./poshelix[4]);
911
912   Double_t  neghelix[6];
913   ntrack->GetHelixParameters(neghelix,b);
914   Double_t negradius = TMath::Abs(1./neghelix[4]);
915
916   Double_t xpos = helixcenterpos[0];
917   Double_t ypos = helixcenterpos[1];
918   Double_t xneg = helixcenterneg[0];
919   Double_t yneg = helixcenterneg[1];
920
921   convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
922   convpos[1] = (ypos*negradius+  yneg*posradius)/(negradius+posradius);
923
924   return 1;
925 }
926
927
928
929 Double_t AliV0Reader::GetConvPosZ(AliESDtrack* ptrack,AliESDtrack* ntrack, Double_t b){
930   //see header file for documentation
931
932   Double_t  helixpos[6];
933   ptrack->GetHelixParameters(helixpos,b);
934
935   Double_t  helixneg[6];
936   ntrack->GetHelixParameters(helixneg,b);
937
938   Double_t negtrackradius =  TMath::Abs(1./helixneg[4]);
939   Double_t postrackradius =  TMath::Abs(1./helixpos[4]);
940
941   Double_t pi = 3.14159265358979323846;
942
943   Double_t convpos[2];
944   GetConvPosXY(ptrack,ntrack,b,convpos);
945
946    Double_t convposx = convpos[0];
947    Double_t convposy = convpos[1];
948
949    Double_t helixcenterpos[2];
950    GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
951
952    Double_t helixcenterneg[2];
953    GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
954
955    Double_t xpos = helixcenterpos[0];
956    Double_t ypos = helixcenterpos[1];
957    Double_t xneg = helixcenterneg[0];
958    Double_t yneg = helixcenterneg[1];
959
960    Double_t deltaXPos = convposx -  xpos;
961    Double_t deltaYPos = convposy -  ypos;
962
963    Double_t deltaXNeg = convposx -  xneg;
964    Double_t deltaYNeg = convposy -  yneg;
965
966    Double_t alphaPos =  pi + TMath::ATan2(-deltaYPos,-deltaXPos);
967    Double_t alphaNeg =  pi + TMath::ATan2(-deltaYNeg,-deltaXNeg);
968
969    Double_t vertexXNeg =  xneg +  TMath::Abs(negtrackradius)*
970    TMath::Cos(alphaNeg);
971    Double_t vertexYNeg =  yneg +  TMath::Abs(negtrackradius)*
972    TMath::Sin(alphaNeg);
973
974    Double_t vertexXPos =  xpos +  TMath::Abs(postrackradius)*
975    TMath::Cos(alphaPos);
976    Double_t vertexYPos =  ypos +  TMath::Abs(postrackradius)*
977    TMath::Sin(alphaPos);
978
979    Double_t x0neg =   helixneg[5];
980    Double_t y0neg =   helixneg[0];
981
982    Double_t x0pos =   helixpos[5];
983    Double_t y0pos =   helixpos[0];
984
985    Double_t dNeg = TMath::Sqrt((vertexXNeg -  x0neg)*(vertexXNeg - x0neg)
986                                +(vertexYNeg -  y0neg)*(vertexYNeg - y0neg));
987
988    Double_t dPos = TMath::Sqrt((vertexXPos -  x0pos)*(vertexXPos - x0pos)
989                                +(vertexYPos -  y0pos)*(vertexYPos - y0pos));
990
991    Double_t rNeg =  TMath::Sqrt(negtrackradius*negtrackradius -
992    dNeg*dNeg/4.);
993
994    Double_t rPos = TMath::Sqrt(postrackradius*postrackradius -
995    dPos*dPos/4.);
996
997    Double_t deltabetaNeg =  2*(pi +   TMath::ATan2(-dNeg/2.,-rNeg));
998    Double_t deltabetaPos = 2*(pi + TMath::ATan2(-dPos/2.,-rPos));
999
1000    Double_t deltaUNeg = negtrackradius*deltabetaNeg;
1001    Double_t deltaUPos = postrackradius*deltabetaPos;
1002
1003    Double_t zphaseNeg = ntrack->GetZ() +  deltaUNeg * ntrack->GetTgl();
1004    Double_t zphasePos = ptrack->GetZ() +  deltaUPos * ptrack->GetTgl();
1005
1006    Double_t convposz =
1007    (zphasePos*negtrackradius+zphaseNeg*postrackradius)/(negtrackradius+postrackradius);
1008
1009    return convposz;
1010 }
1011
1012 AliGammaConversionKFVector* AliV0Reader::GetBGGoodV0s(Int_t event){
1013
1014   return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),fESDEvent->GetNumberOfTracks());
1015 }