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