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