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