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