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