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 "AliPID.h"
  30 #include "AliESDtrack.h"
  31 #include "AliMCEvent.h"
  32 #include "AliKFVertex.h"
  33
  34 #include "AliStack.h"
  35 #include "AliMCEventHandler.h"
  36 #include "AliESDpid.h"
  37 #include "AliGammaConversionBGHandler.h"
  38 #include "AliESDtrackCuts.h"
  39
  40
  41 class iostream;
  42 class AliESDv0;
  43 class TFormula;
  44
  45 using namespace std;
  46
  47 ClassImp(AliV0Reader)
  48
  49
  50 AliESDpid* AliV0Reader::fgESDpid = 0x0;
  51
  52 AliV0Reader::AliV0Reader() :
  53   TObject(),
  54   fMCStack(NULL),
  55   // fMCTruth(NULL),
  56   fMCEvent(NULL),    // for CF
  57   fChain(NULL),
  58   // fESDHandler(NULL),
  59   fESDEvent(NULL),
  60   fCFManager(NULL),
  61   //fESDpid(NULL),
  62   fHistograms(NULL),
  63   fCurrentV0IndexNumber(0),
  64   fCurrentV0(NULL),
  65   fCurrentNegativeKFParticle(NULL),
  66   fCurrentPositiveKFParticle(NULL),
  67   fCurrentMotherKFCandidate(NULL),
  68   fCurrentNegativeESDTrack(NULL),
  69   fCurrentPositiveESDTrack(NULL),
  70   fNegativeTrackLorentzVector(NULL),
  71   fPositiveTrackLorentzVector(NULL),
  72   fMotherCandidateLorentzVector(NULL),
  73   fCurrentXValue(0),
  74   fCurrentYValue(0),
  75   fCurrentZValue(0),
  76   fPositiveTrackPID(0),
  77   fNegativeTrackPID(0),
  78   fNegativeMCParticle(NULL),
  79   fPositiveMCParticle(NULL),
  80   fMotherMCParticle(NULL),
  81   fMotherCandidateKFMass(0),
  82   fMotherCandidateKFWidth(0),
  83   fUseKFParticle(kTRUE),
  84   fUseESDTrack(kFALSE),
  85   fDoMC(kFALSE),
  86   fMaxVertexZ(100.),// 100 cm(from the 0)
  87   fMaxR(10000),// 100 meter(outside of ALICE)
  88   fMinR(0),// 100 meter(outside of ALICE)
  89   fEtaCut(0.),
  90   fRapidityMesonCut(0.),
  91   fPtCut(0.),
  92   fSinglePtCut(0.),
  93   fMaxZ(0.),
  94   fMinClsTPC(0.),
  95   fMinClsTPCToF(0.),
  96   fLineCutZRSlope(0.),
  97   fLineCutZValue(0.),
  98   fChi2CutConversion(0.),
  99   fChi2CutMeson(0.),
 100   fAlphaCutMeson(1.),
 101   fAlphaMinCutMeson(0.),
 102   fPIDProbabilityCutNegativeParticle(0),
 103   fPIDProbabilityCutPositiveParticle(0),
 104   fDodEdxSigmaCut(kFALSE),
 105   fPIDnSigmaAboveElectronLine(100),
 106   fPIDnSigmaBelowElectronLine(-100),
 107   fPIDnSigmaAbovePionLine(-100),
 108   fPIDMinPnSigmaAbovePionLine(100),
 109   fPIDMaxPnSigmaAbovePionLine(100),
 110   fDoKaonRejectionLowP(kFALSE),
 111   fDoProtonRejectionLowP(kFALSE),
 112   fDoPionRejectionLowP(kFALSE),
 113   fPIDnSigmaAtLowPAroundKaonLine(0),
 114   fPIDnSigmaAtLowPAroundProtonLine(0),
 115   fPIDnSigmaAtLowPAroundPionLine(0),
 116   fPIDMinPKaonRejectionLowP(0),
 117   fPIDMinPProtonRejectionLowP(0),
 118   fPIDMinPPionRejectionLowP(0),
 119   fDoQtGammaSelection(kFALSE),
 120   fDoHighPtQtGammaSelection(kFALSE), // RRnew
 121   fQtMax(100.),
 122   fHighPtQtMax(100.), // RRnew
 123   fPtBorderForQt(100.), // RRnew
 124   fXVertexCut(0.),
 125   fYVertexCut(0.),
 126   fZVertexCut(0.),
 127   fNSigmaMass(0.),
 128   fUseImprovedVertex(kFALSE),
 129   fUseOwnXYZCalculation(kFALSE),
 130   fUseConstructGamma(kFALSE),
 131   fDoCF(kFALSE),
 132   fUseOnFlyV0Finder(kTRUE),
 133   fUpdateV0AlreadyCalled(kFALSE),
 134   fCurrentEventGoodV0s(NULL),
 135   fV0Pindex(),
 136   fV0Nindex(),
 137 //  fPreviousEventGoodV0s(),
 138   fCalculateBackground(kFALSE),
 139   fBGEventHandler(NULL),
 140   fBGEventInitialized(kFALSE),
 141   fEsdTrackCuts(NULL),
 142   fNumberOfESDTracks(0),
 143   fNEventsForBGCalculation(20),
 144   fUseChargedTrackMultiplicityForBG(kTRUE),
 145   fNumberOfGoodV0s(0),
 146   fIsHeavyIon(0),
 147   fUseCorrectedTPCClsInfo(kFALSE)
 148 {
 149   //fESDpid = new AliESDpid;
 150 }
 151
 152
 153 AliV0Reader::AliV0Reader(const AliV0Reader & original) :
 154   TObject(original),
 155   fMCStack(original.fMCStack),
 156   // fMCTruth(original.fMCTruth),
 157   fMCEvent(original.fMCEvent),  // for CF
 158   fChain(original.fChain),
 159   //  fESDHandler(original.fESDHandler),
 160   fESDEvent(original.fESDEvent),
 161   fCFManager(original.fCFManager),
 162   // fESDpid(original.fESDpid),
 163   fHistograms(original.fHistograms),
 164   fCurrentV0IndexNumber(original.fCurrentV0IndexNumber),
 165   fCurrentV0(original.fCurrentV0),
 166   fCurrentNegativeKFParticle(original.fCurrentNegativeKFParticle),
 167   fCurrentPositiveKFParticle(original.fCurrentPositiveKFParticle),
 168   fCurrentMotherKFCandidate(original.fCurrentMotherKFCandidate),
 169   fCurrentNegativeESDTrack(original.fCurrentNegativeESDTrack),
 170   fCurrentPositiveESDTrack(original.fCurrentPositiveESDTrack),
 171   fNegativeTrackLorentzVector(original.fNegativeTrackLorentzVector),
 172   fPositiveTrackLorentzVector(original.fPositiveTrackLorentzVector),
 173   fMotherCandidateLorentzVector(original.fMotherCandidateLorentzVector),
 174   fCurrentXValue(original.fCurrentXValue),
 175   fCurrentYValue(original.fCurrentYValue),
 176   fCurrentZValue(original.fCurrentZValue),
 177   fPositiveTrackPID(original.fPositiveTrackPID),
 178   fNegativeTrackPID(original.fNegativeTrackPID),
 179   fNegativeMCParticle(original.fNegativeMCParticle),
 180   fPositiveMCParticle(original.fPositiveMCParticle),
 181   fMotherMCParticle(original.fMotherMCParticle),
 182   fMotherCandidateKFMass(original.fMotherCandidateKFMass),
 183   fMotherCandidateKFWidth(original.fMotherCandidateKFWidth),
 184   fUseKFParticle(kTRUE),
 185   fUseESDTrack(kFALSE),
 186   fDoMC(kFALSE),
 187   fMaxVertexZ(original.fMaxVertexZ),
 188   fMaxR(original.fMaxR),
 189   fMinR(original.fMinR),
 190   fEtaCut(original.fEtaCut),
 191   fRapidityMesonCut(original.fRapidityMesonCut),
 192   fPtCut(original.fPtCut),
 193   fSinglePtCut(original.fSinglePtCut),
 194   fMaxZ(original.fMaxZ),
 195   fMinClsTPC(original.fMinClsTPC),
 196   fMinClsTPCToF(original.fMinClsTPCToF),
 197   fLineCutZRSlope(original.fLineCutZRSlope),
 198   fLineCutZValue(original.fLineCutZValue),
 199   fChi2CutConversion(original.fChi2CutConversion),
 200   fChi2CutMeson(original.fChi2CutMeson),
 201   fAlphaCutMeson(original.fAlphaCutMeson),
 202   fAlphaMinCutMeson(original.fAlphaMinCutMeson),
 203   fPIDProbabilityCutNegativeParticle(original.fPIDProbabilityCutNegativeParticle),
 204   fPIDProbabilityCutPositiveParticle(original.fPIDProbabilityCutPositiveParticle),
 205   fDodEdxSigmaCut(original.fDodEdxSigmaCut),
 206   fPIDnSigmaAboveElectronLine(original.fPIDnSigmaAboveElectronLine),
 207   fPIDnSigmaBelowElectronLine(original.fPIDnSigmaBelowElectronLine),
 208   fPIDnSigmaAbovePionLine(original.fPIDnSigmaAbovePionLine),
 209   fPIDMinPnSigmaAbovePionLine(original.fPIDMinPnSigmaAbovePionLine),
 210   fPIDMaxPnSigmaAbovePionLine(original.fPIDMaxPnSigmaAbovePionLine),
 211   fDoKaonRejectionLowP(original.fDoKaonRejectionLowP),
 212   fDoProtonRejectionLowP(original.fDoProtonRejectionLowP),
 213   fDoPionRejectionLowP(original.fDoPionRejectionLowP),
 214   fPIDnSigmaAtLowPAroundKaonLine(original.fPIDnSigmaAtLowPAroundKaonLine),
 215   fPIDnSigmaAtLowPAroundProtonLine(original.fPIDnSigmaAtLowPAroundProtonLine),
 216   fPIDnSigmaAtLowPAroundPionLine(original.fPIDnSigmaAtLowPAroundPionLine),
 217   fPIDMinPKaonRejectionLowP(original.fPIDMinPKaonRejectionLowP),
 218   fPIDMinPProtonRejectionLowP(original.fPIDMinPProtonRejectionLowP),
 219   fPIDMinPPionRejectionLowP(original.fPIDMinPPionRejectionLowP),
 220   fDoQtGammaSelection(original.fDoQtGammaSelection),
 221   fDoHighPtQtGammaSelection(original.fDoHighPtQtGammaSelection), // RRnew
 222   fQtMax(original.fQtMax),
 223   fHighPtQtMax(original.fHighPtQtMax), // RRnew
 224   fPtBorderForQt(original.fPtBorderForQt), // RRnew
 225   fXVertexCut(original.fXVertexCut),
 226   fYVertexCut(original.fYVertexCut),
 227   fZVertexCut(original.fZVertexCut),
 228   fNSigmaMass(original.fNSigmaMass),
 229   fUseImprovedVertex(original.fUseImprovedVertex),
 230   fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
 231   fUseConstructGamma(original.fUseConstructGamma),
 232   fDoCF(original.fDoCF),
 233   fUseOnFlyV0Finder(original.fUseOnFlyV0Finder),
 234   fUpdateV0AlreadyCalled(original.fUpdateV0AlreadyCalled),
 235   fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),
 236   fV0Pindex(original.fV0Pindex),
 237   fV0Nindex(original.fV0Nindex),
 238   //  fPreviousEventGoodV0s(original.fPreviousEventGoodV0s),
 239   fCalculateBackground(original.fCalculateBackground),
 240   fBGEventHandler(original.fBGEventHandler),
 241   fBGEventInitialized(original.fBGEventInitialized),
 242   fEsdTrackCuts(original.fEsdTrackCuts),
 243   fNumberOfESDTracks(original.fNumberOfESDTracks),
 244   fNEventsForBGCalculation(original.fNEventsForBGCalculation),
 245   fUseChargedTrackMultiplicityForBG(original.fUseChargedTrackMultiplicityForBG),
 246   fNumberOfGoodV0s(original.fNumberOfGoodV0s),
 247   fIsHeavyIon(original.fIsHeavyIon),
 248   fUseCorrectedTPCClsInfo(original.fUseCorrectedTPCClsInfo)
 249 {
 250
 251 }
 252
 253
 254 AliV0Reader & AliV0Reader::operator = (const AliV0Reader & /*source*/)
 255 {
 256   // assignment operator
 257   return *this;
 258 }
 259 AliV0Reader::~AliV0Reader()
 260 {
 261   //  if(fESDpid){
 262   // delete fESDpid;
 263   //}
 264 }
 265
 266 //____________________________________________________________________________
 267 void AliV0Reader::SetInputAndMCEvent(AliVEvent* esd, AliMCEvent* mc) {
 268   // Connect the data pointers
 269
 270   SetInputEvent(esd);
 271   SetMC(mc);
 272
 273 }
 274
 275
 276 void AliV0Reader::Initialize(){
 277   //see header file for documentation
 278
 279   fUpdateV0AlreadyCalled = kFALSE;
 280
 281   /*
 282   // Get the input handler from the manager
 283   fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 284   if(fESDHandler == NULL){
 285     //print warning here
 286   }
 287
 288   // Get pointer to esd event from input handler
 289   fESDEvent = fESDHandler->GetEvent();
 290   if(fESDEvent == NULL){
 291     //print warning here
 292   }
 293
 294   //Get pointer to MCTruth
 295   fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 296   */
 297
 298
 299
 300   //  fMCTruth = mcH->MCEvent();
 301   //  fMC = mcH->MCEvent();
 302   // stack = fMC->Stack();
 303
 304
 305   //if(fMCTruth == NULL){
 306     //print warning here
 307   // fDoMC = kFALSE;
 308   //}
 309
 310   if(fMCEvent == NULL){
 311    fDoMC = kFALSE;
 312   }
 313
 314   //Get pointer to the mc stack
 315   //  if(fMCTruth){
 316   if(fMCEvent){
 317     fMCStack = fMCEvent->Stack();
 318     if(fMCStack == NULL){
 319       //print warning here
 320     }
 321     // Better parameters for MonteCarlo from A. Kalweit 2010/01/8
 322 //     fESDpid->GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
 323 //                                    1.75295e+01,
 324 //                                    3.40030e-09,
 325 //                                    1.96178e+00,
 326 //                                    3.91720e+00);
 327   }
 328   else{
 329     // Better parameters for data from A. Kalweit 2010/01/8
 330  //    fESDpid->GetTPCResponse().SetBetheBlochParameters(0.0283086,
 331 //                                   2.63394e+01,
 332 //                                   5.04114e-11,
 333 //                                   2.12543e+00,
 334 //                                   4.88663e+00);
 335   }
 336
 337   // for CF
 338   //Get pointer to the mc event
 339   if(fDoCF && fDoMC){
 340     //fMCEvent = fMCTruth->MCEvent();
 341     if(fMCEvent == NULL){
 342       //print warning here
 343       fDoCF = kFALSE;
 344     }
 345   }
 346
 347   AliKFParticle::SetField(fESDEvent->GetMagneticField());
 348
 349   //  fCurrentEventGoodV0s = new TClonesArray("TClonesArray", 0);
 350   if(fCurrentEventGoodV0s == NULL){
 351     fCurrentEventGoodV0s = new TClonesArray("AliKFParticle", 0);
 352   }
 353
 354   fV0Pindex.clear();
 355   fV0Nindex.clear();
 356
 357   if(fCalculateBackground == kTRUE){
 358     if(fBGEventInitialized == kFALSE){
 359
 360
 361       Double_t *zBinLimitsArray = new Double_t[9];
 362       zBinLimitsArray[0] = -50.00;
 363       zBinLimitsArray[1] = -3.375;
 364       zBinLimitsArray[2] = -1.605;
 365       zBinLimitsArray[3] = -0.225;
 366       zBinLimitsArray[4] = 1.065;
 367       zBinLimitsArray[5] = 2.445;
 368       zBinLimitsArray[6] = 4.245;
 369       zBinLimitsArray[7] = 50.00;
 370       zBinLimitsArray[8] = 1000.00;
 371
 372       Double_t *multiplicityBinLimitsArray= new Double_t[6];
 373       if(fUseChargedTrackMultiplicityForBG == kTRUE){
 374         multiplicityBinLimitsArray[0] = 0;
 375         multiplicityBinLimitsArray[1] = 8.5;
 376         multiplicityBinLimitsArray[2] = 16.5;
 377         multiplicityBinLimitsArray[3] = 27.5;
 378         multiplicityBinLimitsArray[4] = 41.5;
 379         multiplicityBinLimitsArray[5] = 100.;
 380         if(fIsHeavyIon){
 381           multiplicityBinLimitsArray[0] = 0;
 382           multiplicityBinLimitsArray[1] = 200.;
 383           multiplicityBinLimitsArray[2] = 500.;
 384           multiplicityBinLimitsArray[3] = 1000.;
 385           multiplicityBinLimitsArray[4] = 1500.;
 386           multiplicityBinLimitsArray[5] = 3000.;
 387         }
 388         fBGEventHandler = new AliGammaConversionBGHandler(9,6,fNEventsForBGCalculation);
 389       }
 390       else{
 391         multiplicityBinLimitsArray[0] = 2;
 392         multiplicityBinLimitsArray[1] = 3;
 393         multiplicityBinLimitsArray[2] = 4;
 394         multiplicityBinLimitsArray[3] = 5;
 395         multiplicityBinLimitsArray[4] = 9999;
 396         if(fIsHeavyIon){
 397           multiplicityBinLimitsArray[0] = 2;
 398           multiplicityBinLimitsArray[1] = 10;
 399           multiplicityBinLimitsArray[2] = 30;
 400           multiplicityBinLimitsArray[3] = 50;
 401           multiplicityBinLimitsArray[4] = 9999;
 402         }
 403
 404         fBGEventHandler = new AliGammaConversionBGHandler(9,5,fNEventsForBGCalculation);
 405       }
 406
 407
 408
 409       /*
 410       // ---------------------------------
 411       Double_t *zBinLimitsArray = new Double_t[1];
 412       zBinLimitsArray[0] = 999999.00;
 413
 414       Double_t *multiplicityBinLimitsArray= new Double_t[1];
 415       multiplicityBinLimitsArray[0] = 99999999.00;
 416       fBGEventHandler = new AliGammaConversionBGHandler(1,1,10);
 417       // ---------------------------------
 418       */
 419       fBGEventHandler->Initialize(zBinLimitsArray, multiplicityBinLimitsArray);
 420       fBGEventInitialized = kTRUE;
 421     }
 422   }
 423 }
 424
 425 AliESDv0* AliV0Reader::GetV0(Int_t index){
 426   //see header file for documentation
 427   fCurrentV0 = fESDEvent->GetV0(index);
 428   UpdateV0Information();
 429   return fCurrentV0;
 430 }
 431
 432 Int_t AliV0Reader::GetNumberOfContributorsVtx(){
 433   // returns number of contributors to the vertex
 434   if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
 435     return fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
 436   }
 437
 438   if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
 439     //    return 0;
 440     //-AM test pi0s without SPD only vertex
 441     if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
 442       return fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
 443
 444     }
 445     if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
 446       cout<<"number of contributors from bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 447       return 0;
 448     }
 449   }
 450   return 0;
 451 }
 452 Bool_t AliV0Reader::CheckForPrimaryVertex(){
 453   //see headerfile for documentation
 454
 455   if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
 456     return 1;
 457   }
 458
 459   if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
 460   // SPD vertex
 461     if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
 462       // return 0;
 463       //-AM test pi0s without SPD only vertex
 464       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 465       return 1;
 466
 467     }
 468     if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
 469       //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 470       return 0;
 471     }
 472   }
 473   return 0;
 474   //  return fESDEvent->GetPrimaryVertex()->GetNContributors()>0;
 475 }
 476
 477 Bool_t AliV0Reader::CheckForPrimaryVertexZ(){
 478   //see headerfile for documentation
 479
 480   if(TMath::Abs(fESDEvent->GetPrimaryVertex()->GetZ())<GetMaxVertexZ()){
 481     return kTRUE;
 482   }else{
 483     return kFALSE;
 484   }
 485   return kTRUE;
 486 }
 487
 488 Bool_t AliV0Reader::CheckV0FinderStatus(Int_t index){
 489   // see headerfile for documentation
 490   if(fUseOnFlyV0Finder){
 491     if(!GetV0(index)->GetOnFlyStatus()){
 492       return kFALSE;
 493     }
 494   }
 495   if(!fUseOnFlyV0Finder){
 496     if(GetV0(index)->GetOnFlyStatus()){
 497       return kFALSE;
 498     }
 499   }
 500   return kTRUE;
 501 }
 502
 503
 504
 505 Bool_t AliV0Reader::NextV0(){
 506   //see header file for documentation
 507
 508   Bool_t iResult=kFALSE;
 509   while(fCurrentV0IndexNumber<fESDEvent->GetNumberOfV0s()){
 510     fCurrentV0 = fESDEvent->GetV0(fCurrentV0IndexNumber);
 511
 512     fUpdateV0AlreadyCalled=kFALSE;
 513
 514     if(fHistograms != NULL){
 515       fHistograms->FillHistogram("ESD_AllV0s_InvMass",GetMotherCandidateMass());
 516     }
 517
 518     // moved it up here so that the correction framework can access pt and eta information
 519     if(UpdateV0Information() == kFALSE){
 520       fCurrentV0IndexNumber++;
 521       continue;
 522     }
 523
 524     Double_t containerInput[3];
 525     if(fDoCF){
 526       containerInput[0] = GetMotherCandidatePt();
 527       containerInput[1] = GetMotherCandidateEta();
 528       containerInput[2] = GetMotherCandidateMass();
 529     }
 530     /*
 531     if(fDoCF){
 532       containerInput[0] = GetMotherCandidatePt();
 533       containerInput[1] = GetMotherCandidateEta();
 534       containerInput[2] = GetMotherCandidateMass();
 535
 536       fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign);           // for CF
 537       fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit);           // for CF
 538       fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks);              // for CF
 539     }
 540     */
 541
 542     //checks if on the fly mode is set
 543     if ( !CheckV0FinderStatus(fCurrentV0IndexNumber) ){
 544       if(fHistograms != NULL){
 545         fHistograms->FillHistogram("ESD_CutGetOnFly_InvMass",GetMotherCandidateMass());
 546       }
 547       fCurrentV0IndexNumber++;
 548       continue;
 549     }
 550     if(fDoCF){
 551       fCFManager->GetParticleContainer()->Fill(containerInput,kStepGetOnFly);           // for CF
 552     }
 553
 554     if(fHistograms != NULL){
 555       fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_InvMass",GetMotherCandidateMass());
 556     }
 557
 558     Double_t armenterosQtAlfa[2];
 559     GetArmenterosQtAlfa(GetNegativeKFParticle(),
 560                         GetPositiveKFParticle(),
 561                         GetMotherCandidateKFCombination(),
 562                         armenterosQtAlfa);
 563
 564     fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
 565
 566
 567     if(fCurrentNegativeESDTrack->Charge() == fCurrentPositiveESDTrack->Charge()){             // avoid like sign
 568       //  iResult=kFALSE;
 569       if(fHistograms != NULL ){
 570         fHistograms->FillHistogram("ESD_CutLikeSign_InvMass",GetMotherCandidateMass());
 571         // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 572         // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 573         //      fUpdateV0AlreadyCalled = kTRUE;
 574       }
 575       fCurrentV0IndexNumber++;
 576       continue;
 577     }
 578     if(fDoCF){
 579       fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign);           // for CF
 580     }
 581
 582
 583     if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ||
 584         !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ){
 585       //  if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kITSrefit) ||
 586       //      !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kITSrefit) ){
 587       //  iResult=kFALSE;
 588       if(fHistograms != NULL){
 589         fHistograms->FillHistogram("ESD_CutRefit_InvMass",GetMotherCandidateMass());
 590         // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 591         // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 592         //fUpdateV0AlreadyCalled = kTRUE;
 593       }
 594       fCurrentV0IndexNumber++;
 595       continue;
 596     }
 597     if(fDoCF){
 598       fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit);           // for CF
 599     }
 600
 601
 602
 603     if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 ||
 604         fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {
 605       //iResult=kFALSE;
 606       if(fHistograms != NULL ){
 607         fHistograms->FillHistogram("ESD_CutKink_InvMass",GetMotherCandidateMass());
 608         // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 609         // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 610         //fUpdateV0AlreadyCalled = kTRUE;
 611       }
 612       fCurrentV0IndexNumber++;
 613       continue;
 614     }
 615
 616     if(fDoCF){
 617       fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks);              // for CF
 618     }
 619
 620     fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_goodtracks_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
 621
 622     if(fDodEdxSigmaCut == kTRUE){
 623       if( fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
 624           fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ||
 625           fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
 626           fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ){
 627         //iResult=kFALSE;
 628         if(fHistograms != NULL ){
 629           fHistograms->FillHistogram("ESD_CutdEdxSigmaElectronLine_InvMass",GetMotherCandidateMass());
 630           // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 631           // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 632           //fUpdateV0AlreadyCalled = kTRUE;
 633         }
 634         fCurrentV0IndexNumber++;
 635         continue;
 636       }
 637       if(fDoCF){
 638         fCFManager->GetParticleContainer()->Fill(containerInput,kStepdEdxElectronselection);               // for CF
 639       }
 640
 641       if( fCurrentPositiveESDTrack->P()>fPIDMinPnSigmaAbovePionLine && fCurrentPositiveESDTrack->P()<fPIDMaxPnSigmaAbovePionLine ){
 642         if(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
 643            fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
 644            fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
 645           //      iResult=kFALSE;
 646           if(fHistograms != NULL){
 647             fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
 648             // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 649             // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 650             //fUpdateV0AlreadyCalled = kTRUE;
 651           }
 652           fCurrentV0IndexNumber++;
 653           continue;
 654         }
 655       }
 656
 657       if( fCurrentNegativeESDTrack->P()>fPIDMinPnSigmaAbovePionLine && fCurrentNegativeESDTrack->P()<fPIDMaxPnSigmaAbovePionLine){
 658         if(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
 659            fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
 660            fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
 661           //      iResult=kFALSE;
 662           if(fHistograms != NULL){
 663             fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
 664             // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 665             // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 666             //fUpdateV0AlreadyCalled = kTRUE;
 667           }
 668           fCurrentV0IndexNumber++;
 669           continue;
 670         }
 671       }
 672       if(fDoCF){
 673         fCFManager->GetParticleContainer()->Fill(containerInput,kStepdEdxPionrejection);               // for CF
 674       }
 675
 676     }
 677
 678     if(fDoKaonRejectionLowP == kTRUE){
 679       if( fCurrentNegativeESDTrack->P()<fPIDMinPKaonRejectionLowP ){
 680         if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kKaon))<fPIDnSigmaAtLowPAroundKaonLine){
 681           if(fHistograms != NULL){
 682             fHistograms->FillHistogram("ESD_CutKaonRejectionLowP_InvMass",GetMotherCandidateMass());
 683             // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 684             // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 685             //fUpdateV0AlreadyCalled = kTRUE;
 686           }
 687           fCurrentV0IndexNumber++;
 688           continue;
 689         }
 690       }
 691       if( fCurrentPositiveESDTrack->P()<fPIDMinPKaonRejectionLowP ){
 692         if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kKaon))<fPIDnSigmaAtLowPAroundKaonLine){
 693           if(fHistograms != NULL){
 694             fHistograms->FillHistogram("ESD_CutKaonRejectionLowP_InvMass",GetMotherCandidateMass());
 695             // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 696             // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 697             //fUpdateV0AlreadyCalled = kTRUE;
 698           }
 699           fCurrentV0IndexNumber++;
 700           continue;
 701         }
 702       }
 703     }
 704
 705     if(fDoProtonRejectionLowP == kTRUE){
 706       if( fCurrentNegativeESDTrack->P()<fPIDMinPProtonRejectionLowP){
 707         if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kProton))<fPIDnSigmaAtLowPAroundProtonLine){
 708           if(fHistograms != NULL){
 709             fHistograms->FillHistogram("ESD_CutProtonRejectionLowP_InvMass",GetMotherCandidateMass());
 710             // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 711             // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 712             //fUpdateV0AlreadyCalled = kTRUE;
 713           }
 714           fCurrentV0IndexNumber++;
 715           continue;
 716         }
 717       }
 718       if( fCurrentPositiveESDTrack->P()<fPIDMinPProtonRejectionLowP ){
 719         if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kProton))<fPIDnSigmaAtLowPAroundProtonLine){
 720           if(fHistograms != NULL){
 721             fHistograms->FillHistogram("ESD_CutProtonRejectionLowP_InvMass",GetMotherCandidateMass());
 722             // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 723             // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 724             //fUpdateV0AlreadyCalled = kTRUE;
 725           }
 726           fCurrentV0IndexNumber++;
 727           continue;
 728         }
 729       }
 730
 731     }
 732
 733     if(fDoPionRejectionLowP == kTRUE){
 734       if( fCurrentNegativeESDTrack->P()<fPIDMinPPionRejectionLowP ){
 735         if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion))<fPIDnSigmaAtLowPAroundPionLine){
 736           if(fHistograms != NULL){
 737             fHistograms->FillHistogram("ESD_CutPionRejectionLowP_InvMass",GetMotherCandidateMass());
 738             // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 739             // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 740             //fUpdateV0AlreadyCalled = kTRUE;
 741           }
 742           fCurrentV0IndexNumber++;
 743           continue;
 744         }
 745       }
 746       if( fCurrentPositiveESDTrack->P()<fPIDMinPPionRejectionLowP ){
 747         if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion))<fPIDnSigmaAtLowPAroundPionLine){
 748           if(fHistograms != NULL){
 749             fHistograms->FillHistogram("ESD_CutPionRejectionLowP_InvMass",GetMotherCandidateMass());
 750             // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 751             // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 752             //fUpdateV0AlreadyCalled = kTRUE;
 753           }
 754           fCurrentV0IndexNumber++;
 755           continue;
 756         }
 757       }
 758     }
 759
 760
 761     // Gamma selection based on QT from Armenteros
 762     if(fDoQtGammaSelection == kTRUE){ // RRnew start : apply different qT-cut above/below
 763         if(fDoHighPtQtGammaSelection){
 764                 if(GetMotherCandidatePt() < fPtBorderForQt){
 765                         if(armenterosQtAlfa[0]>fQtMax){
 766                                 if(fHistograms != NULL){
 767                                 fHistograms->FillHistogram("ESD_CutQt_InvMass",GetMotherCandidateMass());
 768                                 }
 769                                 fCurrentV0IndexNumber++;
 770                                 continue;
 771                         }
 772                 }
 773                 else{
 774                         if(armenterosQtAlfa[0]>fHighPtQtMax)  {
 775                                 if(fHistograms != NULL){
 776                                 fHistograms->FillHistogram("ESD_CutQt_InvMass",GetMotherCandidateMass());
 777                                 }
 778                                 fCurrentV0IndexNumber++;
 779                                 continue;
 780                         }
 781                 }
 782
 783         }
 784         else{
 785                 if(armenterosQtAlfa[0]>fQtMax){
 786                         if(fHistograms != NULL){
 787                         fHistograms->FillHistogram("ESD_CutQt_InvMass",GetMotherCandidateMass());
 788                         }
 789                         fCurrentV0IndexNumber++;
 790                         continue;
 791                 }
 792         }
 793     }   // RRnew end
 794
 795     //checks if we have a prim vertex
 796     //if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {
 797     if(GetNumberOfContributorsVtx()<=0) {
 798       if(fHistograms != NULL){
 799         fHistograms->FillHistogram("ESD_CutNContributors_InvMass",GetMotherCandidateMass());
 800       }
 801       fCurrentV0IndexNumber++;
 802       continue;
 803     }
 804     if(fDoCF){
 805       fCFManager->GetParticleContainer()->Fill(containerInput,kStepNContributors);              // for CF
 806     }
 807
 808     //Check the pid probability
 809     if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){
 810       if(fHistograms != NULL){
 811         fHistograms->FillHistogram("ESD_CutPIDProb_InvMass",GetMotherCandidateMass());
 812       }
 813       fCurrentV0IndexNumber++;
 814       continue;
 815     }
 816     if(fDoCF){
 817       fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCPID);                     // for CF
 818     }
 819
 820     if(GetXYRadius()>fMaxR){ // cuts on distance from collision point
 821       if(fHistograms != NULL){
 822         fHistograms->FillHistogram("ESD_CutR_InvMass",GetMotherCandidateMass());
 823       }
 824       fCurrentV0IndexNumber++;
 825       continue;
 826     }
 827     if(fDoCF){
 828       fCFManager->GetParticleContainer()->Fill(containerInput,kStepR);                  // for CF
 829     }
 830     if(GetXYRadius()<fMinR){ // cuts on distance from collision point
 831       if(fHistograms != NULL){
 832         fHistograms->FillHistogram("ESD_CutMinR_InvMass",GetMotherCandidateMass());
 833       }
 834       fCurrentV0IndexNumber++;
 835       continue;
 836     }
 837
 838
 839
 840     if((TMath::Abs(fCurrentZValue)*fLineCutZRSlope)-fLineCutZValue > GetXYRadius() ){ // cuts out regions where we do not reconstruct
 841       if(fHistograms != NULL){
 842         fHistograms->FillHistogram("ESD_CutLine_InvMass",GetMotherCandidateMass());
 843       }
 844       fCurrentV0IndexNumber++;
 845       continue;
 846     }
 847     if(fDoCF){
 848       fCFManager->GetParticleContainer()->Fill(containerInput,kStepLine);                       // for CF
 849     }
 850
 851     if(TMath::Abs(fCurrentZValue) > fMaxZ ){ // cuts out regions where we do not reconstruct
 852       if(fHistograms != NULL){
 853         fHistograms->FillHistogram("ESD_CutZ_InvMass",GetMotherCandidateMass());
 854       }
 855       fCurrentV0IndexNumber++;
 856       continue;
 857     }
 858     if(fDoCF){
 859       fCFManager->GetParticleContainer()->Fill(containerInput,kStepZ);          // for CF
 860     }
 861
 862     /* Moved further up so corr framework can work
 863        if(UpdateV0Information() == kFALSE){
 864        fCurrentV0IndexNumber++;
 865        continue;
 866        }
 867     */
 868     if(fCurrentNegativeESDTrack->GetNcls(1) < fMinClsTPC ||  fCurrentPositiveESDTrack->GetNcls(1) < fMinClsTPC ){
 869       if(fHistograms != NULL){
 870         fHistograms->FillHistogram("ESD_CutMinNClsTPC_InvMass",GetMotherCandidateMass());
 871       }
 872       fCurrentV0IndexNumber++;
 873       continue;
 874     }
 875     if(fDoCF){
 876       fCFManager->GetParticleContainer()->Fill(containerInput,kStepMinClsTPC);          // for CF
 877     }
 878     Double_t negclsToF = 0.;
 879     if (!fUseCorrectedTPCClsInfo ){
 880       if(fCurrentNegativeESDTrack->GetTPCNclsF()!=0  ){
 881         negclsToF = (Double_t)fCurrentNegativeESDTrack->GetNcls(1)/(Double_t)fCurrentNegativeESDTrack->GetTPCNclsF();
 882       }
 883     }else{
 884       negclsToF = fCurrentNegativeESDTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
 885     }
 886
 887     Double_t posclsToF = 0.;
 888     if (!fUseCorrectedTPCClsInfo ){
 889       if(fCurrentPositiveESDTrack->GetTPCNclsF()!=0  ){
 890         posclsToF = (Double_t)fCurrentPositiveESDTrack->GetNcls(1)/(Double_t)fCurrentPositiveESDTrack->GetTPCNclsF();
 891       }
 892     }else{
 893       posclsToF = fCurrentPositiveESDTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
 894     }
 895
 896     if( negclsToF < fMinClsTPCToF ||  posclsToF < fMinClsTPCToF ){
 897       if(fHistograms != NULL){
 898         fHistograms->FillHistogram("ESD_CutMinNClsTPCToF_InvMass",GetMotherCandidateMass());
 899       }
 900       fCurrentV0IndexNumber++;
 901       continue;
 902     }
 903
 904
 905
 906
 907     if(fUseKFParticle){
 908
 909
 910       if( fCurrentNegativeKFParticle->GetPt()< fSinglePtCut ||  fCurrentPositiveKFParticle->GetPt()< fSinglePtCut){
 911         if(fHistograms != NULL){
 912           fHistograms->FillHistogram("ESD_CutSinglePt_InvMass",GetMotherCandidateMass());
 913         }
 914         fCurrentV0IndexNumber++;
 915         continue;
 916       }
 917       if(fDoCF){
 918         fCFManager->GetParticleContainer()->Fill(containerInput,kStepSinglePt);         // for CF
 919       }
 920
 921
 922       if(fCurrentMotherKFCandidate->GetNDF()<=0){
 923         if(fHistograms != NULL){
 924           fHistograms->FillHistogram("ESD_CutNDF_InvMass",GetMotherCandidateMass());
 925         }
 926         fCurrentV0IndexNumber++;
 927         continue;
 928       }
 929       if(fDoCF){
 930         fCFManager->GetParticleContainer()->Fill(containerInput,kStepNDF);              // for CF
 931       }
 932
 933       Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();
 934       if(chi2V0 > fChi2CutConversion || chi2V0 <=0){
 935         if(fHistograms != NULL){
 936           fHistograms->FillHistogram("ESD_CutChi2_InvMass",GetMotherCandidateMass());
 937         }
 938         fCurrentV0IndexNumber++;
 939         continue;
 940       }
 941       if(fDoCF){
 942         fCFManager->GetParticleContainer()->Fill(containerInput,kStepChi2);                     // for CF
 943       }
 944
 945       if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut){
 946         if(fHistograms != NULL){
 947           fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
 948         }
 949         fCurrentV0IndexNumber++;
 950         continue;
 951       }
 952
 953       if(TMath::Abs(fCurrentNegativeKFParticle->GetEta())> fEtaCut){
 954         if(fHistograms != NULL){
 955           fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
 956         }
 957         fCurrentV0IndexNumber++;
 958         continue;
 959       }
 960
 961       if(TMath::Abs(fCurrentPositiveKFParticle->GetEta())> fEtaCut){
 962         if(fHistograms != NULL){
 963           fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
 964         }
 965         fCurrentV0IndexNumber++;
 966         continue;
 967       }
 968
 969       if(fDoCF){
 970         fCFManager->GetParticleContainer()->Fill(containerInput,kStepEta);                      // for CF
 971       }
 972
 973       if(fMotherCandidateLorentzVector->Pt()<fPtCut){
 974         if(fHistograms != NULL){
 975           fHistograms->FillHistogram("ESD_CutPt_InvMass",GetMotherCandidateMass());
 976         }
 977         fCurrentV0IndexNumber++;
 978         continue;
 979       }
 980       if(fDoCF){
 981         fCFManager->GetParticleContainer()->Fill(containerInput,kStepPt);                       // for CF
 982       }
 983
 984     }
 985     else if(fUseESDTrack){
 986       //TODO
 987     }
 988
 989     if(fHistograms != NULL){
 990       fHistograms->FillHistogram("ESD_GoodV0s_InvMass",GetMotherCandidateMass());
 991     }
 992
 993     //    fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate);
 994
 995     new((*fCurrentEventGoodV0s)[fCurrentEventGoodV0s->GetEntriesFast()])  AliKFParticle(*fCurrentMotherKFCandidate);
 996     fV0Pindex.push_back(fCurrentV0->GetPindex());
 997     fV0Nindex.push_back(fCurrentV0->GetNindex());
 998
 999     iResult=kTRUE;//means we have a v0 who survived all the cuts applied
1000
1001     fNumberOfGoodV0s++;
1002
1003     fCurrentV0IndexNumber++;
1004
1005     break;
1006   }
1007   return iResult;
1008 }
1009
1010 Bool_t AliV0Reader::UpdateV0Information(){
1011   //see header file for documentation
1012
1013   Bool_t iResult=kTRUE;                                         // for taking out not refitted, kinks and like sign tracks
1014
1015   Bool_t switchTracks = kFALSE;
1016
1017   fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1018   fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1019
1020   if(fCurrentPositiveESDTrack->GetSign() == -1 && fCurrentNegativeESDTrack->GetSign() == 1){  // switch wrong signed tracks
1021     fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1022     fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1023     switchTracks = kTRUE;
1024   }
1025
1026   if(fCurrentNegativeKFParticle != NULL){
1027     delete fCurrentNegativeKFParticle;
1028   }
1029   if(switchTracks == kFALSE){
1030     fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fNegativeTrackPID);
1031   }
1032   else{
1033     fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fNegativeTrackPID);
1034   }
1035
1036   if(fCurrentPositiveKFParticle != NULL){
1037     delete fCurrentPositiveKFParticle;
1038   }
1039   if(switchTracks == kFALSE){
1040     fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fPositiveTrackPID);
1041   }
1042   else{
1043     fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fPositiveTrackPID);
1044   }
1045
1046   if(fCurrentMotherKFCandidate != NULL){
1047     delete fCurrentMotherKFCandidate;
1048   }
1049
1050   if(fUseConstructGamma==kTRUE){
1051     fCurrentMotherKFCandidate = new AliKFParticle;//(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
1052     fCurrentMotherKFCandidate->ConstructGamma(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
1053   }else{
1054     fCurrentMotherKFCandidate = new AliKFParticle(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
1055     if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
1056       fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);
1057     }
1058   }
1059   if(fUseImprovedVertex == kTRUE){
1060     AliKFVertex primaryVertexImproved(*GetPrimaryVertex());
1061     primaryVertexImproved+=*fCurrentMotherKFCandidate;
1062     fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);
1063   }
1064
1065   fCurrentMotherKFCandidate->GetMass(fMotherCandidateKFMass,fMotherCandidateKFWidth);
1066
1067   if(fNegativeTrackLorentzVector != NULL){
1068     delete fNegativeTrackLorentzVector;
1069   }
1070   if(fUseKFParticle){
1071     fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeKFParticle->Px(),fCurrentNegativeKFParticle->Py(),fCurrentNegativeKFParticle->Pz());
1072   }
1073   else if(fUseESDTrack){
1074     fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeESDTrack->Px(),fCurrentNegativeESDTrack->Py(),fCurrentNegativeESDTrack->Pz());
1075   }
1076
1077   if(fPositiveTrackLorentzVector != NULL){
1078     delete fPositiveTrackLorentzVector;
1079   }
1080   if(fUseKFParticle){
1081     fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveKFParticle->Px(),fCurrentPositiveKFParticle->Py(),fCurrentPositiveKFParticle->Pz());
1082   }
1083   else if(fUseESDTrack){
1084     fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveESDTrack->Px(),fCurrentPositiveESDTrack->Py(),fCurrentPositiveESDTrack->Pz());
1085   }
1086
1087   if(fMotherCandidateLorentzVector != NULL){
1088     delete fMotherCandidateLorentzVector;
1089   }
1090   if(fUseKFParticle){
1091     fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
1092   }
1093   else if(fUseESDTrack){
1094     fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
1095   }
1096
1097   if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
1098     fMotherCandidateLorentzVector->SetXYZM(fMotherCandidateLorentzVector->Px() ,fMotherCandidateLorentzVector->Py(),fMotherCandidateLorentzVector->Pz(),0.);
1099   }
1100
1101
1102   if(fDoMC == kTRUE){
1103     fMotherMCParticle= NULL;
1104     fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
1105     fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
1106     if(fPositiveMCParticle->GetMother(0)>-1){
1107       fMotherMCParticle = fMCStack->Particle(fPositiveMCParticle->GetMother(0));
1108     }
1109   }
1110
1111
1112   // for CF
1113 //   Double_t containerInput[3];
1114 //   if(fDoCF){
1115 //     containerInput[0] = GetMotherCandidatePt();
1116 //     containerInput[1] = GetMotherCandidateEta();
1117 //     containerInput[2] = GetMotherCandidateMass();
1118
1119 //     fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign);          // for CF
1120 //     fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit);          // for CF
1121 //     fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks);             // for CF
1122 //   }
1123
1124
1125   if(fUseOwnXYZCalculation == kFALSE){
1126     if(fUseConstructGamma == kFALSE){
1127       fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);
1128     }else{
1129       fCurrentXValue=GetMotherCandidateKFCombination()->GetX();
1130       fCurrentYValue=GetMotherCandidateKFCombination()->GetY();
1131       fCurrentZValue=GetMotherCandidateKFCombination()->GetZ();
1132     }
1133   }
1134   else{
1135     Double_t convpos[2];
1136     convpos[0]=0;
1137     convpos[1]=0;
1138
1139     GetConvPosXY(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField(),convpos);
1140     fCurrentXValue = convpos[0];
1141     fCurrentYValue = convpos[1];
1142     fCurrentZValue = GetConvPosZ(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField());
1143   }
1144   /*
1145   if(fCurrentNegativeESDTrack->GetSign() == fCurrentPositiveESDTrack->GetSign()){             // avoid like sign
1146     iResult=kFALSE;
1147     if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE && doFillHistos == kTRUE){
1148       fHistograms->FillHistogram("ESD_CutLikeSign_InvMass",GetMotherCandidateMass());
1149       // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1150       // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1151       fUpdateV0AlreadyCalled = kTRUE;
1152     }
1153   }
1154
1155   if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ||
1156       !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ){
1157     //  if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kITSrefit) ||
1158     //      !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kITSrefit) ){
1159     iResult=kFALSE;
1160     if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE  && doFillHistos == kTRUE){
1161       fHistograms->FillHistogram("ESD_CutRefit_InvMass",GetMotherCandidateMass());
1162       // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1163       // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1164       fUpdateV0AlreadyCalled = kTRUE;
1165     }
1166   }
1167
1168   if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 ||
1169       fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {
1170
1171     iResult=kFALSE;
1172     if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE && doFillHistos == kTRUE ){
1173       fHistograms->FillHistogram("ESD_CutKink_InvMass",GetMotherCandidateMass());
1174       // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1175       // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1176       fUpdateV0AlreadyCalled = kTRUE;
1177     }
1178   }
1179
1180   if(fDodEdxSigmaCut == kTRUE){
1181
1182     if( fESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
1183         fESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ||
1184         fESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
1185         fESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ){
1186       iResult=kFALSE;
1187       if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE  && doFillHistos == kTRUE){
1188         fHistograms->FillHistogram("ESD_CutdEdxSigmaElectronLine_InvMass",GetMotherCandidateMass());
1189         // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1190         // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1191         fUpdateV0AlreadyCalled = kTRUE;
1192       }
1193     }
1194     if( fCurrentPositiveESDTrack->P()>fPIDMinPnSigmaAbovePionLine){
1195       if(fESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
1196          fESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
1197          fESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
1198         iResult=kFALSE;
1199         if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE  && doFillHistos == kTRUE){
1200           fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
1201           // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1202           // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1203           fUpdateV0AlreadyCalled = kTRUE;
1204         }
1205       }
1206     }
1207
1208     if( fCurrentNegativeESDTrack->P()>fPIDMinPnSigmaAbovePionLine){
1209       if(fESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
1210          fESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
1211          fESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
1212         iResult=kFALSE;
1213         if(fHistograms != NULL && fUpdateV0AlreadyCalled == kFALSE && doFillHistos == kTRUE ){
1214           fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
1215           // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1216           // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1217           fUpdateV0AlreadyCalled = kTRUE;
1218         }
1219       }
1220     }
1221   }
1222   */
1223   fUpdateV0AlreadyCalled = kTRUE;
1224
1225   return iResult;
1226 }
1227
1228
1229
1230 Bool_t AliV0Reader::HasSameMCMother(){
1231   //see header file for documentation
1232
1233   Bool_t iResult = kFALSE;
1234   if(fDoMC == kTRUE){
1235     if(fNegativeMCParticle != NULL && fPositiveMCParticle != NULL){
1236       if(fNegativeMCParticle->GetMother(0) == fPositiveMCParticle->GetMother(0))
1237         if(fMotherMCParticle){
1238           iResult = kTRUE;
1239         }
1240     }
1241   }
1242   return iResult;
1243 }
1244
1245 Bool_t AliV0Reader::CheckPIDProbability(Double_t negProbCut, Double_t posProbCut){
1246   //see header file for documentation
1247
1248   Bool_t iResult=kFALSE;
1249
1250   //  Double_t *posProbArray = new Double_t[10];
1251   //  Double_t *negProbArray = new Double_t[10];
1252   //-AM The TPCpid method expects an array of length kSPECIES that is 5 not 10
1253
1254   Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1255   Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1256
1257   AliESDtrack* negTrack  = GetNegativeESDTrack();
1258   AliESDtrack* posTrack  = GetPositiveESDTrack();
1259   //fESDEvent->GetTrack(fCurrentV0->GetNindex());
1260     //fESDEvent->GetTrack(fCurrentV0->GetPindex());
1261   //-AM for switchtracks==true the above is a bug
1262
1263   negTrack->GetTPCpid(negProbArray);
1264   posTrack->GetTPCpid(posProbArray);
1265
1266   //  if(negProbArray != NULL && posProbArray != NULL){ // this is not allowed anymore for some reason(RC19)
1267   if(negProbArray && posProbArray){
1268     if(negProbArray[GetSpeciesIndex(-1)]>=negProbCut && posProbArray[GetSpeciesIndex(1)]>=posProbCut){
1269       iResult=kTRUE;
1270     }
1271   }
1272   delete [] posProbArray;
1273   delete [] negProbArray;
1274   return iResult;
1275 }
1276
1277 void AliV0Reader::GetPIDProbability(Double_t &negPIDProb,Double_t & posPIDProb){
1278   // see header file for documentation
1279
1280   //Double_t *posProbArray = new Double_t[10];
1281   // Double_t *negProbArray = new Double_t[10];
1282   //-AM The TPCpid method expects an array of length kSPECIES that is 5 not 10
1283   Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1284   Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1285
1286 //   AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1287 //   AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1288   //-AM for switchtracks the above is a bug
1289   AliESDtrack* negTrack  = GetNegativeESDTrack();
1290   AliESDtrack* posTrack  = GetPositiveESDTrack();
1291
1292
1293   negTrack->GetTPCpid(negProbArray);
1294   posTrack->GetTPCpid(posProbArray);
1295
1296   //  if(negProbArray!=NULL && posProbArray!=NULL){ // this is not allowed anymore for some reason(RC19)
1297   if(negProbArray && posProbArray){
1298     negPIDProb = negProbArray[GetSpeciesIndex(-1)];
1299     posPIDProb = posProbArray[GetSpeciesIndex(1)];
1300   }
1301   delete [] posProbArray;
1302   delete [] negProbArray;
1303 }
1304 void AliV0Reader::GetPIDProbabilityMuonPion(Double_t &negPIDProb,Double_t & posPIDProb){
1305   // see header file for documentation
1306
1307
1308   Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1309   Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1310
1311   // AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1312   // AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1313   //-AM for switchtracks the above is a bug
1314
1315   AliESDtrack* negTrack  = GetNegativeESDTrack();
1316   AliESDtrack* posTrack  = GetPositiveESDTrack();
1317
1318   negTrack->GetTPCpid(negProbArray);
1319   posTrack->GetTPCpid(posProbArray);
1320
1321   //  if(negProbArray!=NULL && posProbArray!=NULL){ // this is not allowed anymore for some reason(RC19)
1322   if(negProbArray && posProbArray){
1323     negPIDProb = negProbArray[1]+negProbArray[2];
1324     posPIDProb = posProbArray[1]+posProbArray[2];
1325   }
1326   delete [] posProbArray;
1327   delete [] negProbArray;
1328 }
1329
1330 void AliV0Reader::UpdateEventByEventData(){
1331   //see header file for documentation
1332   if(fCurrentEventGoodV0s->GetEntriesFast() >0 ){
1333     if(fCalculateBackground){
1334       if(fUseChargedTrackMultiplicityForBG == kTRUE){
1335         fBGEventHandler->AddEvent(fCurrentEventGoodV0s,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
1336         //filling z and multiplicity histograms
1337         fHistograms->FillHistogram("ESD_Z_distribution",fESDEvent->GetPrimaryVertex()->GetZ());
1338         fHistograms->FillHistogram("ESD_multiplicity_distribution",CountESDTracks());
1339         fHistograms->FillHistogram("ESD_ZvsMultiplicity",fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
1340       }
1341       else{ // means we use #V0s for multiplicity
1342         fBGEventHandler->AddEvent(fCurrentEventGoodV0s,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
1343         //filling z and multiplicity histograms
1344         fHistograms->FillHistogram("ESD_Z_distribution",fESDEvent->GetPrimaryVertex()->GetZ());
1345         fHistograms->FillHistogram("ESD_multiplicity_distribution",fNumberOfGoodV0s);
1346         fHistograms->FillHistogram("ESD_ZvsMultiplicity",fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
1347       }
1348     }
1349   }
1350   fCurrentEventGoodV0s->Delete();
1351   fCurrentV0IndexNumber=0;
1352   fNumberOfESDTracks=0;
1353
1354   fV0Pindex.clear();
1355   fV0Nindex.clear();
1356
1357
1358
1359   //  fBGEventHandler->PrintBGArray(); // for debugging
1360 }
1361
1362
1363 Double_t AliV0Reader::GetNegativeTrackPhi() const{
1364   //see header file for documentation
1365
1366   Double_t offset=0;
1367   if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){
1368     offset = -2*TMath::Pi();
1369   }
1370   return fNegativeTrackLorentzVector->Phi()+offset;
1371 }
1372
1373 Double_t AliV0Reader::GetPositiveTrackPhi() const{
1374   //see header file for documentation
1375
1376   Double_t offset=0;
1377   if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){
1378     offset = -2*TMath::Pi();
1379   }
1380   return fPositiveTrackLorentzVector->Phi()+offset;
1381 }
1382
1383 Double_t AliV0Reader::GetMotherCandidatePhi() const{
1384   //see header file for documentation
1385
1386   Double_t offset=0;
1387   if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){
1388     offset = -2*TMath::Pi();
1389   }
1390   return fMotherCandidateLorentzVector->Phi()+offset;
1391 }
1392
1393
1394 Double_t AliV0Reader::GetMotherCandidateRapidity() const{
1395   //see header file for documentation
1396
1397   Double_t rapidity=0;
1398   if(fMotherCandidateLorentzVector->Energy() - fMotherCandidateLorentzVector->Pz() == 0 || fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz() == 0) rapidity=0;
1399   else rapidity = 0.5*(TMath::Log((fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz()) / (fMotherCandidateLorentzVector->Energy()-fMotherCandidateLorentzVector->Pz())));
1400   return rapidity;
1401
1402 }
1403
1404
1405
1406
1407
1408 Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){
1409   //see header file for documentation
1410
1411   Int_t iResult = 10; // Unknown particle
1412
1413   if(chargeOfTrack==-1){ //negative track
1414     switch(abs(fNegativeTrackPID)){
1415     case 11:       //electron
1416       iResult = 0;
1417       break;
1418     case 13:       //muon
1419       iResult = 1;
1420       break;
1421     case 211:      //pion
1422       iResult = 2;
1423       break;
1424     case 321:      //kaon
1425       iResult = 3;
1426       break;
1427     case 2212:     //proton
1428       iResult = 4;
1429       break;
1430     case 22:       //photon
1431       iResult = 5;
1432       break;
1433     case 111:      //pi0
1434       iResult = 6;
1435       break;
1436     case 2112:     //neutron
1437       iResult = 7;
1438       break;
1439     case 311:      //K0
1440       iResult = 8;
1441       break;
1442
1443       //Put in here for kSPECIES::kEleCon  ????
1444     }
1445   }
1446   else if(chargeOfTrack==1){ //positive track
1447     switch(abs(fPositiveTrackPID)){
1448     case 11:       //electron
1449       iResult = 0;
1450       break;
1451     case 13:       //muon
1452       iResult = 1;
1453       break;
1454     case 211:      //pion
1455       iResult = 2;
1456       break;
1457     case 321:      //kaon
1458       iResult = 3;
1459       break;
1460     case 2212:     //proton
1461       iResult = 4;
1462       break;
1463     case 22:       //photon
1464       iResult = 5;
1465       break;
1466     case 111:      //pi0
1467       iResult = 6;
1468       break;
1469     case 2112:     //neutron
1470       iResult = 7;
1471       break;
1472     case 311:      //K0
1473       iResult = 8;
1474       break;
1475
1476       //Put in here for kSPECIES::kEleCon  ????
1477     }
1478   }
1479   else{
1480     //Wrong parameter.. Print warning
1481   }
1482   return iResult;
1483 }
1484
1485 Bool_t AliV0Reader::GetHelixCenter(AliESDtrack* track, Double_t b,Int_t charge, Double_t center[2]){
1486   // see header file for documentation
1487
1488   Double_t pi = 3.14159265358979323846;
1489
1490   Double_t  helix[6];
1491   track->GetHelixParameters(helix,b);
1492
1493   Double_t xpos =  helix[5];
1494   Double_t ypos =  helix[0];
1495   Double_t radius = TMath::Abs(1./helix[4]);
1496   Double_t phi = helix[2];
1497
1498   if(phi < 0){
1499     phi = phi + 2*pi;
1500   }
1501
1502   phi -= pi/2.;
1503   Double_t xpoint =  radius * TMath::Cos(phi);
1504   Double_t ypoint =  radius * TMath::Sin(phi);
1505
1506   if(b<0){
1507     if(charge > 0){
1508       xpoint = - xpoint;
1509       ypoint = - ypoint;
1510     }
1511
1512     if(charge < 0){
1513       xpoint =  xpoint;
1514       ypoint =  ypoint;
1515     }
1516   }
1517   if(b>0){
1518     if(charge > 0){
1519       xpoint =  xpoint;
1520       ypoint =  ypoint;
1521     }
1522
1523     if(charge < 0){
1524       xpoint = - xpoint;
1525       ypoint = - ypoint;
1526     }
1527   }
1528   center[0] =  xpos + xpoint;
1529   center[1] =  ypos + ypoint;
1530
1531   return 1;
1532 }
1533
1534 Bool_t AliV0Reader::GetConvPosXY(AliESDtrack* ptrack, AliESDtrack* ntrack, Double_t b, Double_t convpos[2]){
1535   //see header file for documentation
1536
1537   Double_t helixcenterpos[2];
1538   GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
1539
1540   Double_t helixcenterneg[2];
1541   GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
1542
1543   Double_t  poshelix[6];
1544   ptrack->GetHelixParameters(poshelix,b);
1545   Double_t posradius = TMath::Abs(1./poshelix[4]);
1546
1547   Double_t  neghelix[6];
1548   ntrack->GetHelixParameters(neghelix,b);
1549   Double_t negradius = TMath::Abs(1./neghelix[4]);
1550
1551   Double_t xpos = helixcenterpos[0];
1552   Double_t ypos = helixcenterpos[1];
1553   Double_t xneg = helixcenterneg[0];
1554   Double_t yneg = helixcenterneg[1];
1555
1556   convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1557   convpos[1] = (ypos*negradius+  yneg*posradius)/(negradius+posradius);
1558
1559   return 1;
1560 }
1561
1562
1563
1564 Double_t AliV0Reader::GetConvPosZ(AliESDtrack* ptrack,AliESDtrack* ntrack, Double_t b){
1565   //see header file for documentation
1566
1567   Double_t  helixpos[6];
1568   ptrack->GetHelixParameters(helixpos,b);
1569
1570   Double_t  helixneg[6];
1571   ntrack->GetHelixParameters(helixneg,b);
1572
1573   Double_t negtrackradius =  TMath::Abs(1./helixneg[4]);
1574   Double_t postrackradius =  TMath::Abs(1./helixpos[4]);
1575
1576   Double_t pi = 3.14159265358979323846;
1577
1578   Double_t convpos[2];
1579   GetConvPosXY(ptrack,ntrack,b,convpos);
1580
1581    Double_t convposx = convpos[0];
1582    Double_t convposy = convpos[1];
1583
1584    Double_t helixcenterpos[2];
1585    GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
1586
1587    Double_t helixcenterneg[2];
1588    GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
1589
1590    Double_t xpos = helixcenterpos[0];
1591    Double_t ypos = helixcenterpos[1];
1592    Double_t xneg = helixcenterneg[0];
1593    Double_t yneg = helixcenterneg[1];
1594
1595    Double_t deltaXPos = convposx -  xpos;
1596    Double_t deltaYPos = convposy -  ypos;
1597
1598    Double_t deltaXNeg = convposx -  xneg;
1599    Double_t deltaYNeg = convposy -  yneg;
1600
1601    Double_t alphaPos =  pi + TMath::ATan2(-deltaYPos,-deltaXPos);
1602    Double_t alphaNeg =  pi + TMath::ATan2(-deltaYNeg,-deltaXNeg);
1603
1604    Double_t vertexXNeg =  xneg +  TMath::Abs(negtrackradius)*
1605    TMath::Cos(alphaNeg);
1606    Double_t vertexYNeg =  yneg +  TMath::Abs(negtrackradius)*
1607    TMath::Sin(alphaNeg);
1608
1609    Double_t vertexXPos =  xpos +  TMath::Abs(postrackradius)*
1610    TMath::Cos(alphaPos);
1611    Double_t vertexYPos =  ypos +  TMath::Abs(postrackradius)*
1612    TMath::Sin(alphaPos);
1613
1614    Double_t x0neg =   helixneg[5];
1615    Double_t y0neg =   helixneg[0];
1616
1617    Double_t x0pos =   helixpos[5];
1618    Double_t y0pos =   helixpos[0];
1619
1620    Double_t dNeg = TMath::Sqrt((vertexXNeg -  x0neg)*(vertexXNeg - x0neg)
1621                                +(vertexYNeg -  y0neg)*(vertexYNeg - y0neg));
1622
1623    Double_t dPos = TMath::Sqrt((vertexXPos -  x0pos)*(vertexXPos - x0pos)
1624                                +(vertexYPos -  y0pos)*(vertexYPos - y0pos));
1625
1626    Double_t rNeg =  TMath::Sqrt(negtrackradius*negtrackradius -
1627    dNeg*dNeg/4.);
1628
1629    Double_t rPos = TMath::Sqrt(postrackradius*postrackradius -
1630    dPos*dPos/4.);
1631
1632    Double_t deltabetaNeg =  2*(pi +   TMath::ATan2(-dNeg/2.,-rNeg));
1633    Double_t deltabetaPos = 2*(pi + TMath::ATan2(-dPos/2.,-rPos));
1634
1635    Double_t deltaUNeg = negtrackradius*deltabetaNeg;
1636    Double_t deltaUPos = postrackradius*deltabetaPos;
1637
1638    Double_t zphaseNeg = ntrack->GetZ() +  deltaUNeg * ntrack->GetTgl();
1639    Double_t zphasePos = ptrack->GetZ() +  deltaUPos * ptrack->GetTgl();
1640
1641    Double_t convposz = (zphasePos*negtrackradius+zphaseNeg*postrackradius)/(negtrackradius+postrackradius);
1642
1643    return convposz;
1644 }
1645
1646 AliGammaConversionKFVector* AliV0Reader::GetBGGoodV0s(Int_t /*event*/) const{
1647   /*
1648   if(fUseChargedTrackMultiplicityForBG == kTRUE){
1649     return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
1650   }
1651   else{ // means we use #v0s as multiplicity
1652     return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
1653   }
1654   */
1655   return NULL;
1656 }
1657
1658 Int_t AliV0Reader::CountESDTracks(){
1659   // see header file for documentation
1660   if(fNumberOfESDTracks == 0){ // count the good esd tracks
1661     for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1662       AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1663       if(!curTrack){
1664         continue;
1665       }
1666       if(fEsdTrackCuts->AcceptTrack(curTrack) ){
1667         fNumberOfESDTracks++;
1668       }
1669     }
1670   }
1671
1672   return fNumberOfESDTracks;
1673 }
1674
1675 Bool_t AliV0Reader::CheckIfPi0IsMother(Int_t label){
1676   // see headerfile for documentation
1677   Bool_t iResult=kFALSE;
1678   //  cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
1679   if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 111){
1680     iResult=kTRUE;
1681   }
1682   return iResult;
1683 }
1684
1685 Bool_t AliV0Reader::CheckIfEtaIsMother(Int_t label){
1686   // see headerfile for documentation
1687   Bool_t iResult=kFALSE;
1688   //  cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
1689   if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 221){
1690     iResult=kTRUE;
1691   }
1692   return iResult;
1693 }
1694
1695
1696 Bool_t AliV0Reader::GetArmenterosQtAlfa(const AliKFParticle* positiveKFParticle, const AliKFParticle * negativeKFParticle, const AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlfa[2] ){
1697   //see header file for documentation
1698
1699   TVector3 momentumVectorPositiveKF(positiveKFParticle->GetPx(),positiveKFParticle->GetPy(),positiveKFParticle->GetPz());
1700   TVector3 momentumVectorNegativeKF(negativeKFParticle->GetPx(),negativeKFParticle->GetPy(),negativeKFParticle->GetPz());
1701   TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
1702
1703   Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
1704   Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
1705
1706   Float_t alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
1707     ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
1708
1709
1710   Float_t qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
1711
1712   armenterosQtAlfa[0]=qt;
1713   armenterosQtAlfa[1]=alfa;
1714
1715   return 1;
1716
1717 }
1718
1719 Int_t AliV0Reader::GetFirstTPCRow(Double_t radius){
1720
1721
1722   Int_t firstTPCRow=0;
1723   Double_t radiusI  =  84.8;
1724   Double_t radiusO  = 134.6;
1725   Double_t radiusOB = 198.;
1726   Double_t rSizeI   = 0.75;
1727   Double_t rSizeO   = 1.;
1728   Double_t rSizeOB  = 1.5;
1729   Int_t nClsI=63;
1730   Int_t nClsIO=127;
1731
1732   if(radius <= radiusI){
1733     return firstTPCRow;
1734   }
1735   if(radius>radiusI && radius<=radiusO){
1736     firstTPCRow = (Int_t)((radius-radiusI)/rSizeI);
1737   }
1738   if(radius>radiusO && radius<=radiusOB){
1739     firstTPCRow = (Int_t)(nClsI+(radius-radiusO)/rSizeO);
1740   }
1741
1742   if(radius>radiusOB){
1743     firstTPCRow =(Int_t)(nClsIO+(radius-radiusOB)/rSizeOB);
1744   }
1745
1746
1747   return firstTPCRow;
1748 }