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