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