PWGGA/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 #include "AliKFParticle.h"
  34 #include "AliStack.h"
  35 #include "AliMCEventHandler.h"
  36 #include "AliESDpid.h"
  37 #include "AliGammaConversionBGHandler.h"
  38 #include "AliESDtrackCuts.h"
  39 #include "TRandom3.h"
  40 #include "AliGenCocktailEventHeader.h"
  41 #include "TList.h"
  42
  43 class iostream;
  44 class AliESDv0;
  45 class TFormula;
  46 class TRandom3;
  47
  48 using namespace std;
  49
  50 ClassImp(AliV0Reader)
  51
  52
  53 AliESDpid* AliV0Reader::fgESDpid = 0x0;
  54
  55 AliV0Reader::AliV0Reader() :
  56         TObject(),
  57         fMCStack(NULL),
  58         // fMCTruth(NULL),
  59         fMCEvent(NULL),         // for CF
  60         fChain(NULL),
  61         // fESDHandler(NULL),
  62         fESDEvent(NULL),
  63         fCFManager(NULL),
  64         //fESDpid(NULL),
  65         fHistograms(NULL),
  66         fCurrentV0IndexNumber(0),
  67         fCurrentV0(NULL),
  68         fCurrentNegativeKFParticle(NULL),
  69         fCurrentPositiveKFParticle(NULL),
  70         fCurrentMotherKFCandidate(NULL),
  71         fCurrentNegativeESDTrack(NULL),
  72         fCurrentPositiveESDTrack(NULL),
  73         fNegativeTrackLorentzVector(NULL),
  74         fPositiveTrackLorentzVector(NULL),
  75         fMotherCandidateLorentzVector(NULL),
  76         fCurrentXValue(0),
  77         fCurrentYValue(0),
  78         fCurrentZValue(0),
  79         fPositiveTrackPID(0),
  80         fNegativeTrackPID(0),
  81         fNegativeMCParticle(NULL),
  82         fPositiveMCParticle(NULL),
  83         fMotherMCParticle(NULL),
  84         fMotherCandidateKFMass(0),
  85         fMotherCandidateKFWidth(0),
  86         fUseKFParticle(kTRUE),
  87         fUseESDTrack(kFALSE),
  88         fDoMC(kFALSE),
  89         fMaxVertexZ(100.),// 100 cm(from the 0)
  90         fMaxR(10000),// 100 meter(outside of ALICE)
  91         fMinR(0),// 100 meter(outside of ALICE)
  92         fEtaCut(0.),
  93         fEtaCutMin(-0.1),
  94         fRapidityMesonCut(0.),
  95         fPtCut(0.),
  96         fSinglePtCut(0.),
  97         fMaxZ(0.),
  98         fMinClsTPC(0.),
  99         fMinClsTPCToF(0.),
 100         fLineCutZRSlope(0.),
 101         fLineCutZValue(0.),
 102         fLineCutZRSlopeMin(0.),
 103         fLineCutZValueMin(0.),
 104         fChi2CutConversion(0.),
 105         fChi2CutMeson(0.),
 106         fAlphaCutMeson(1.),
 107         fAlphaMinCutMeson(0.),
 108         fPIDProbabilityCutNegativeParticle(0),
 109         fPIDProbabilityCutPositiveParticle(0),
 110         fDodEdxSigmaCut(kFALSE),
 111         fDoTOFsigmaCut(kFALSE), // RRnewTOF
 112         fPIDnSigmaAboveElectronLine(100),
 113         fPIDnSigmaBelowElectronLine(-100),
 114         fTofPIDnSigmaAboveElectronLine(100), // RRnewTOF
 115         fTofPIDnSigmaBelowElectronLine(-100), // RRnewTOF
 116         fPIDnSigmaAbovePionLine(-100),
 117         fPIDnSigmaAbovePionLineHighPt(-100),
 118         fPIDMinPnSigmaAbovePionLine(100),
 119         fPIDMaxPnSigmaAbovePionLine(100),
 120         fDoKaonRejectionLowP(kFALSE),
 121         fDoProtonRejectionLowP(kFALSE),
 122         fDoPionRejectionLowP(kFALSE),
 123         fPIDnSigmaAtLowPAroundKaonLine(0),
 124         fPIDnSigmaAtLowPAroundProtonLine(0),
 125         fPIDnSigmaAtLowPAroundPionLine(0),
 126         fPIDMinPKaonRejectionLowP(0),
 127         fPIDMinPProtonRejectionLowP(0),
 128         fPIDMinPPionRejectionLowP(0),
 129         fDoQtGammaSelection(kFALSE),
 130         fDoHighPtQtGammaSelection(kFALSE), // RRnew
 131         fQtMax(100.),
 132         fHighPtQtMax(100.), // RRnew
 133         fPtBorderForQt(100.), // RRnew
 134         fXVertexCut(0.),
 135         fYVertexCut(0.),
 136         fZVertexCut(0.),
 137         fPsiPairCut(0.),
 138         fCosinePointCut(0.),
 139         fNSigmaMass(0.),
 140         fUseImprovedVertex(kFALSE),
 141         fUseOwnXYZCalculation(kFALSE),
 142         fUseConstructGamma(kFALSE),
 143         fDoCF(kFALSE),
 144         fUseEtaMinCut(kFALSE),
 145         fUseOnFlyV0Finder(kTRUE),
 146         fUpdateV0AlreadyCalled(kFALSE),
 147         fCurrentEventGoodV0s(NULL),
 148         fV0Pindex(),
 149         fV0Nindex(),
 150 //      fPreviousEventGoodV0s(),
 151         fCalculateBackground(kFALSE),
 152         fBGEventHandler(NULL),
 153         fBGEventInitialized(kFALSE),
 154         fEsdTrackCuts(NULL),
 155         fNumberOfESDTracks(0),
 156         fNEventsForBGCalculation(20),
 157         fUseChargedTrackMultiplicityForBG(kTRUE),
 158         fNumberOfGoodV0s(0),
 159         fIsHeavyIon(0),
 160         fUseCorrectedTPCClsInfo(kFALSE),
 161         fUseMCPSmearing(kTRUE),
 162         fPBremSmearing(1.),
 163         fPSigSmearing(0.),
 164         fPSigSmearingCte(0.),
 165         fRandom(0),
 166         fBrem(NULL),
 167         fDoPhotonAsymmetryCut(0),
 168         fdoESDQtCut(0),
 169         fMinPPhotonAsymmetryCut(100.),
 170         fMinPhotonAsymmetry(0.),
 171         fExcludeBackgroundEventForGammaCorrection(0),
 172         fNumberOfPrimerisFromHijingAndPythia(0)
 173 {
 174         //fESDpid = new AliESDpid;
 175 }
 176
 177
 178 AliV0Reader::AliV0Reader(const AliV0Reader & original) :
 179         TObject(original),
 180         fMCStack(original.fMCStack),
 181         // fMCTruth(original.fMCTruth),
 182         fMCEvent(original.fMCEvent),    // for CF
 183         fChain(original.fChain),
 184         //      fESDHandler(original.fESDHandler),
 185         fESDEvent(original.fESDEvent),
 186         fCFManager(original.fCFManager),
 187         // fESDpid(original.fESDpid),
 188         fHistograms(original.fHistograms),
 189         fCurrentV0IndexNumber(original.fCurrentV0IndexNumber),
 190         fCurrentV0(original.fCurrentV0),
 191         fCurrentNegativeKFParticle(original.fCurrentNegativeKFParticle),
 192         fCurrentPositiveKFParticle(original.fCurrentPositiveKFParticle),
 193         fCurrentMotherKFCandidate(original.fCurrentMotherKFCandidate),
 194         fCurrentNegativeESDTrack(original.fCurrentNegativeESDTrack),
 195         fCurrentPositiveESDTrack(original.fCurrentPositiveESDTrack),
 196         fNegativeTrackLorentzVector(original.fNegativeTrackLorentzVector),
 197         fPositiveTrackLorentzVector(original.fPositiveTrackLorentzVector),
 198         fMotherCandidateLorentzVector(original.fMotherCandidateLorentzVector),
 199         fCurrentXValue(original.fCurrentXValue),
 200         fCurrentYValue(original.fCurrentYValue),
 201         fCurrentZValue(original.fCurrentZValue),
 202         fPositiveTrackPID(original.fPositiveTrackPID),
 203         fNegativeTrackPID(original.fNegativeTrackPID),
 204         fNegativeMCParticle(original.fNegativeMCParticle),
 205         fPositiveMCParticle(original.fPositiveMCParticle),
 206         fMotherMCParticle(original.fMotherMCParticle),
 207         fMotherCandidateKFMass(original.fMotherCandidateKFMass),
 208         fMotherCandidateKFWidth(original.fMotherCandidateKFWidth),
 209         fUseKFParticle(kTRUE),
 210         fUseESDTrack(kFALSE),
 211         fDoMC(kFALSE),
 212         fMaxVertexZ(original.fMaxVertexZ),
 213         fMaxR(original.fMaxR),
 214         fMinR(original.fMinR),
 215         fEtaCut(original.fEtaCut),
 216         fEtaCutMin(original.fEtaCutMin),
 217         fRapidityMesonCut(original.fRapidityMesonCut),
 218         fPtCut(original.fPtCut),
 219         fSinglePtCut(original.fSinglePtCut),
 220         fMaxZ(original.fMaxZ),
 221         fMinClsTPC(original.fMinClsTPC),
 222         fMinClsTPCToF(original.fMinClsTPCToF),
 223         fLineCutZRSlope(original.fLineCutZRSlope),
 224         fLineCutZValue(original.fLineCutZValue),
 225         fLineCutZRSlopeMin(original.fLineCutZRSlopeMin),
 226         fLineCutZValueMin(original.fLineCutZValueMin),
 227         fChi2CutConversion(original.fChi2CutConversion),
 228         fChi2CutMeson(original.fChi2CutMeson),
 229         fAlphaCutMeson(original.fAlphaCutMeson),
 230         fAlphaMinCutMeson(original.fAlphaMinCutMeson),
 231         fPIDProbabilityCutNegativeParticle(original.fPIDProbabilityCutNegativeParticle),
 232         fPIDProbabilityCutPositiveParticle(original.fPIDProbabilityCutPositiveParticle),
 233         fDodEdxSigmaCut(original.fDodEdxSigmaCut),
 234         fDoTOFsigmaCut(original.fDoTOFsigmaCut), // RRnewTOF
 235         fPIDnSigmaAboveElectronLine(original.fPIDnSigmaAboveElectronLine),
 236         fPIDnSigmaBelowElectronLine(original.fPIDnSigmaBelowElectronLine),
 237         fTofPIDnSigmaAboveElectronLine(original.fTofPIDnSigmaAboveElectronLine), // RRnewTOF
 238         fTofPIDnSigmaBelowElectronLine(original.fTofPIDnSigmaBelowElectronLine), // RRnewTOF
 239         fPIDnSigmaAbovePionLine(original.fPIDnSigmaAbovePionLine),
 240         fPIDnSigmaAbovePionLineHighPt(original.fPIDnSigmaAbovePionLineHighPt),
 241         fPIDMinPnSigmaAbovePionLine(original.fPIDMinPnSigmaAbovePionLine),
 242         fPIDMaxPnSigmaAbovePionLine(original.fPIDMaxPnSigmaAbovePionLine),
 243         fDoKaonRejectionLowP(original.fDoKaonRejectionLowP),
 244         fDoProtonRejectionLowP(original.fDoProtonRejectionLowP),
 245         fDoPionRejectionLowP(original.fDoPionRejectionLowP),
 246         fPIDnSigmaAtLowPAroundKaonLine(original.fPIDnSigmaAtLowPAroundKaonLine),
 247         fPIDnSigmaAtLowPAroundProtonLine(original.fPIDnSigmaAtLowPAroundProtonLine),
 248         fPIDnSigmaAtLowPAroundPionLine(original.fPIDnSigmaAtLowPAroundPionLine),
 249         fPIDMinPKaonRejectionLowP(original.fPIDMinPKaonRejectionLowP),
 250         fPIDMinPProtonRejectionLowP(original.fPIDMinPProtonRejectionLowP),
 251         fPIDMinPPionRejectionLowP(original.fPIDMinPPionRejectionLowP),
 252         fDoQtGammaSelection(original.fDoQtGammaSelection),
 253         fDoHighPtQtGammaSelection(original.fDoHighPtQtGammaSelection), // RRnew
 254         fQtMax(original.fQtMax),
 255         fHighPtQtMax(original.fHighPtQtMax), // RRnew
 256         fPtBorderForQt(original.fPtBorderForQt), // RRnew
 257         fXVertexCut(original.fXVertexCut),
 258         fYVertexCut(original.fYVertexCut),
 259         fZVertexCut(original.fZVertexCut),
 260         fPsiPairCut(original.fPsiPairCut),
 261         fCosinePointCut(original.fCosinePointCut),
 262         fNSigmaMass(original.fNSigmaMass),
 263         fUseImprovedVertex(original.fUseImprovedVertex),
 264         fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
 265         fUseConstructGamma(original.fUseConstructGamma),
 266         fDoCF(original.fDoCF),
 267         fUseEtaMinCut(original.fUseEtaMinCut),
 268         fUseOnFlyV0Finder(original.fUseOnFlyV0Finder),
 269         fUpdateV0AlreadyCalled(original.fUpdateV0AlreadyCalled),
 270         fCurrentEventGoodV0s(original.fCurrentEventGoodV0s),
 271         fV0Pindex(original.fV0Pindex),
 272         fV0Nindex(original.fV0Nindex),
 273         //      fPreviousEventGoodV0s(original.fPreviousEventGoodV0s),
 274         fCalculateBackground(original.fCalculateBackground),
 275         fBGEventHandler(original.fBGEventHandler),
 276         fBGEventInitialized(original.fBGEventInitialized),
 277         fEsdTrackCuts(original.fEsdTrackCuts),
 278         fNumberOfESDTracks(original.fNumberOfESDTracks),
 279         fNEventsForBGCalculation(original.fNEventsForBGCalculation),
 280         fUseChargedTrackMultiplicityForBG(original.fUseChargedTrackMultiplicityForBG),
 281         fNumberOfGoodV0s(original.fNumberOfGoodV0s),
 282         fIsHeavyIon(original.fIsHeavyIon),
 283         fUseCorrectedTPCClsInfo(original.fUseCorrectedTPCClsInfo),
 284         fUseMCPSmearing(original.fUseMCPSmearing),
 285         fPBremSmearing(original.fPBremSmearing),
 286         fPSigSmearing(original.fPSigSmearing),
 287         fPSigSmearingCte(original.fPSigSmearingCte),
 288         fRandom(original.fRandom),
 289         fBrem(original.fBrem),
 290         fDoPhotonAsymmetryCut(original.fDoPhotonAsymmetryCut),
 291         fdoESDQtCut(original.fdoESDQtCut),
 292         fMinPPhotonAsymmetryCut(original.fMinPPhotonAsymmetryCut),
 293         fMinPhotonAsymmetry(original.fMinPhotonAsymmetry),
 294         fExcludeBackgroundEventForGammaCorrection(original.fExcludeBackgroundEventForGammaCorrection),
 295         fNumberOfPrimerisFromHijingAndPythia(original.fNumberOfPrimerisFromHijingAndPythia)
 296 {
 297
 298 }
 299
 300
 301 AliV0Reader & AliV0Reader::operator = (const AliV0Reader & /*source*/)
 302 {
 303         // assignment operator
 304         return *this;
 305 }
 306 AliV0Reader::~AliV0Reader()
 307 {
 308         //      if(fESDpid){
 309         // delete fESDpid;
 310         //}
 311 }
 312
 313 //____________________________________________________________________________
 314 void AliV0Reader::SetInputAndMCEvent(AliVEvent* esd, AliMCEvent* mc) {
 315         // Connect the data pointers
 316
 317         SetInputEvent(esd);
 318         SetMC(mc);
 319
 320 }
 321
 322
 323 void AliV0Reader::Initialize(){
 324         //see header file for documentation
 325
 326         fUpdateV0AlreadyCalled = kFALSE;
 327
 328         /*
 329         // Get the input handler from the manager
 330         fESDHandler = (AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 331         if(fESDHandler == NULL){
 332                 //print warning here
 333         }
 334
 335         // Get pointer to esd event from input handler
 336         fESDEvent = fESDHandler->GetEvent();
 337         if(fESDEvent == NULL){
 338                 //print warning here
 339         }
 340
 341         //Get pointer to MCTruth
 342         fMCTruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 343         */
 344
 345
 346
 347         //      fMCTruth = mcH->MCEvent();
 348         //      fMC = mcH->MCEvent();
 349         // stack = fMC->Stack();
 350
 351
 352         //if(fMCTruth == NULL){
 353                 //print warning here
 354         // fDoMC = kFALSE;
 355         //}
 356
 357         if(fMCEvent == NULL){
 358                 fDoMC = kFALSE;
 359         }
 360
 361         //Get pointer to the mc stack
 362         //      if(fMCTruth){
 363         if(fMCEvent){
 364                 fMCStack = fMCEvent->Stack();
 365                 // if(fMCStack == NULL){
 366                 //      //print warning here
 367                 // }
 368                 // Better parameters for MonteCarlo from A. Kalweit 2010/01/8
 369 //               fESDpid->GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
 370 //                                                      1.75295e+01,
 371 //                                                      3.40030e-09,
 372 //                                                      1.96178e+00,
 373 //                                                      3.91720e+00);
 374         }
 375         else{
 376                 // Better parameters for data from A. Kalweit 2010/01/8
 377  //             fESDpid->GetTPCResponse().SetBetheBlochParameters(0.0283086,
 378 //                                               2.63394e+01,
 379 //                                               5.04114e-11,
 380 //                                               2.12543e+00,
 381 //                                               4.88663e+00);
 382         }
 383
 384         // for CF
 385         //Get pointer to the mc event
 386         if(fDoCF && fDoMC){
 387                 //fMCEvent = fMCTruth->MCEvent();
 388                 if(fMCEvent == NULL){
 389                         //print warning here
 390                         fDoCF = kFALSE;
 391                 }
 392         }
 393
 394         fUseEtaMinCut = kFALSE;
 395         if ( fEtaCutMin != -0.1) {
 396                 fUseEtaMinCut = kTRUE;
 397         }
 398
 399
 400         AliKFParticle::SetField(fESDEvent->GetMagneticField());
 401
 402         //      fCurrentEventGoodV0s = new TClonesArray("TClonesArray", 0);
 403         if(fCurrentEventGoodV0s == NULL){
 404                 fCurrentEventGoodV0s = new TClonesArray("AliKFParticle", 0);
 405         }
 406
 407         fV0Pindex.clear();
 408         fV0Nindex.clear();
 409
 410 //      if(gRandom != NULL){
 411 //              delete gRandom;
 412 //              gRandom= new TRandom3(0);
 413 //      }
 414
 415
 416         if (fBrem == NULL){
 417                 fBrem = new TF1("fBrem","pow(-log(x),[0]/log(2.0)-1.0)/TMath::Gamma([0]/log(2.0))",0.00001,0.999999999);
 418                 // tests done with 1.0e-14
 419                 fBrem->SetParameter(0,fPBremSmearing);
 420                 fBrem->SetNpx(100000);
 421         }
 422
 423         if(fCalculateBackground == kTRUE){
 424                 if(fBGEventInitialized == kFALSE){
 425
 426                         Double_t *zBinLimitsArray = new Double_t[9];
 427                         zBinLimitsArray[0] = -50.00;
 428                         zBinLimitsArray[1] = -3.375;
 429                         zBinLimitsArray[2] = -1.605;
 430                         zBinLimitsArray[3] = -0.225;
 431                         zBinLimitsArray[4] = 1.065;
 432                         zBinLimitsArray[5] = 2.445;
 433                         zBinLimitsArray[6] = 4.245;
 434                         zBinLimitsArray[7] = 50.00;
 435                         zBinLimitsArray[8] = 1000.00;
 436
 437                    Double_t *multiplicityBinLimitsArray = new Double_t[6];
 438                         if(fUseChargedTrackMultiplicityForBG == kTRUE){
 439                                 multiplicityBinLimitsArray[0] = 0;
 440                                 multiplicityBinLimitsArray[1] = 8.5;
 441                                 multiplicityBinLimitsArray[2] = 16.5;
 442                                 multiplicityBinLimitsArray[3] = 27.5;
 443                                 multiplicityBinLimitsArray[4] = 41.5;
 444                                 multiplicityBinLimitsArray[5] = 100.;
 445                                 if(fIsHeavyIon){
 446                                         multiplicityBinLimitsArray[0] = 0;
 447                                         multiplicityBinLimitsArray[1] = 200.;
 448                                         multiplicityBinLimitsArray[2] = 500.;
 449                                         multiplicityBinLimitsArray[3] = 1000.;
 450                                         multiplicityBinLimitsArray[4] = 1500.;
 451                                         multiplicityBinLimitsArray[5] = 3000.;
 452                                 }
 453                                 fBGEventHandler = new AliGammaConversionBGHandler(9,6,fNEventsForBGCalculation);
 454                         } else {
 455                                 multiplicityBinLimitsArray[0] = 2;
 456                                 multiplicityBinLimitsArray[1] = 3;
 457                                 multiplicityBinLimitsArray[2] = 4;
 458                                 multiplicityBinLimitsArray[3] = 5;
 459                                 multiplicityBinLimitsArray[4] = 9999;
 460                                 if(fIsHeavyIon){
 461                                         multiplicityBinLimitsArray[0] = 2;
 462                                         multiplicityBinLimitsArray[1] = 10;
 463                                         multiplicityBinLimitsArray[2] = 30;
 464                                         multiplicityBinLimitsArray[3] = 50;
 465                                         multiplicityBinLimitsArray[4] = 9999;
 466                                 }
 467
 468                                 fBGEventHandler = new AliGammaConversionBGHandler(9,5,fNEventsForBGCalculation);
 469                         }
 470
 471
 472
 473                         /*
 474                         // ---------------------------------
 475                         Double_t *zBinLimitsArray = new Double_t[1];
 476                         zBinLimitsArray[0] = 999999.00;
 477
 478                         Double_t *multiplicityBinLimitsArray= new Double_t[1];
 479                         multiplicityBinLimitsArray[0] = 99999999.00;
 480                         fBGEventHandler = new AliGammaConversionBGHandler(1,1,10);
 481                         // ---------------------------------
 482                         */
 483                         fBGEventHandler->Initialize(zBinLimitsArray, multiplicityBinLimitsArray);
 484                         fBGEventInitialized = kTRUE;
 485                 }
 486         }
 487
 488         if(fDoMC && fExcludeBackgroundEventForGammaCorrection){
 489            fNumberOfPrimerisFromHijingAndPythia = GetNumberOfHijingPlusPythiaPrimeries(fExcludeBackgroundEventForGammaCorrection);
 490         }
 491 }
 492
 493 AliESDv0* AliV0Reader::GetV0(Int_t index){
 494         //see header file for documentation
 495         fCurrentV0 = fESDEvent->GetV0(index);
 496         UpdateV0Information();
 497         return fCurrentV0;
 498 }
 499
 500 Int_t AliV0Reader::GetNumberOfContributorsVtx(){
 501         // returns number of contributors to the vertex
 502         if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
 503                 return fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
 504         }
 505
 506         if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
 507                 //              return 0;
 508                 //-AM test pi0s without SPD only vertex
 509                 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
 510                         return fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
 511
 512                 }
 513                 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
 514                         cout<<"number of contributors from bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 515                         return 0;
 516                 }
 517         }
 518         return 0;
 519 }
 520
 521 Bool_t AliV0Reader::CheckForPrimaryVertex(){
 522         //see headerfile for documentation
 523         if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
 524                 return 1;
 525         }
 526         if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
 527         // SPD vertex
 528                 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
 529                         // return 0;
 530                         //-AM test pi0s without SPD only vertex
 531                         //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 532                         return 1;
 533                 }
 534                 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
 535                         //                      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 536                         return 0;
 537                 }
 538         }
 539         return 0;
 540         //      return fESDEvent->GetPrimaryVertex()->GetNContributors()>0;
 541 }
 542
 543 Bool_t AliV0Reader::CheckForPrimaryVertexZ(){
 544         //see headerfile for documentation
 545         if(TMath::Abs(fESDEvent->GetPrimaryVertex()->GetZ())<GetMaxVertexZ()){
 546                 return kTRUE;
 547         }else{
 548                 return kFALSE;
 549         }
 550         return kTRUE;
 551 }
 552
 553
 554 Bool_t AliV0Reader::CheckV0FinderStatus(AliESDv0 *v0){
 555         // see headerfile for documentation
 556         if(fUseOnFlyV0Finder){
 557                 if(!v0->GetOnFlyStatus()){
 558                         return kFALSE;
 559                 }
 560         }
 561         if(!fUseOnFlyV0Finder){
 562                 if(v0->GetOnFlyStatus()){
 563                         return kFALSE;
 564                 }
 565         }
 566         return kTRUE;
 567 }
 568
 569
 570
 571
 572 Bool_t AliV0Reader::NextV0(){
 573         //see header file for documentation
 574         Bool_t iResult=kFALSE;
 575         while(fCurrentV0IndexNumber<fESDEvent->GetNumberOfV0s()){
 576                 fCurrentV0 = fESDEvent->GetV0(fCurrentV0IndexNumber);
 577
 578                 fUpdateV0AlreadyCalled=kFALSE;
 579
 580                 if(fHistograms != NULL){
 581                         fHistograms->FillHistogram("ESD_AllV0s_InvMass",GetMotherCandidateMass());
 582                         fHistograms->FillHistogram("ESD_AllV0s_Pt",GetMotherCandidatePt());
 583                 }
 584
 585                 // moved it up here so that the correction framework can access pt and eta information
 586                 if(UpdateV0Information() == kFALSE){
 587                         fCurrentV0IndexNumber++;
 588                         continue;
 589                 }
 590
 591
 592                 if(fDoMC && fExcludeBackgroundEventForGammaCorrection){ // Remove all V0s from BGEvent
 593                    Bool_t isFromBGEvent = kFALSE;
 594                    isFromBGEvent = IsParticleFromBGEvent(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
 595                    if(isFromBGEvent){
 596                       fHistograms->FillHistogram("ESD_CutMCBgEvent_InvMass",GetMotherCandidateMass());
 597                       fHistograms->FillHistogram("ESD_CutMCBgEvent_Pt",GetMotherCandidatePt());
 598                       fCurrentV0IndexNumber++;
 599                       continue;
 600                    }
 601                    isFromBGEvent = IsParticleFromBGEvent(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
 602                    if(isFromBGEvent){
 603                       fHistograms->FillHistogram("ESD_CutMCBgEvent_InvMass",GetMotherCandidateMass());
 604                       fHistograms->FillHistogram("ESD_CutMCBgEvent_Pt",GetMotherCandidatePt());
 605                       fCurrentV0IndexNumber++;
 606                       continue;
 607                    }
 608                 }
 609
 610
 611
 612                 Double_t containerInput[3];
 613                 if(fDoCF){
 614                         containerInput[0] = GetMotherCandidatePt();
 615                         containerInput[1] = GetMotherCandidateEta();
 616                         containerInput[2] = GetMotherCandidateMass();
 617                 }
 618                 /*
 619                 if(fDoCF){
 620                         containerInput[0] = GetMotherCandidatePt();
 621                         containerInput[1] = GetMotherCandidateEta();
 622                         containerInput[2] = GetMotherCandidateMass();
 623
 624                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign);         // for CF
 625                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit);         // for CF
 626                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks);            // for CF
 627                 }
 628                 */
 629
 630                 //checks if on the fly mode is set
 631                 if ( !CheckV0FinderStatus(fCurrentV0) ){
 632
 633                         if(fHistograms != NULL){
 634                                 fHistograms->FillHistogram("ESD_CutGetOnFly_InvMass",GetMotherCandidateMass());
 635                                 fHistograms->FillHistogram("ESD_CutGetOnFly_Pt",GetMotherCandidatePt());
 636                         }
 637                         fCurrentV0IndexNumber++;
 638                         continue;
 639                 }
 640                 if(fDoCF){
 641                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepGetOnFly);         // for CF
 642                 }
 643
 644                 if(fHistograms != NULL){
 645                         fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_InvMass",GetMotherCandidateMass());
 646                         fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_Pt",GetMotherCandidatePt());
 647                 }
 648
 649                 Double_t armenterosQtAlpha[2] = {0,0};
 650                 Double_t armenterosQtAlphaKF[2] = {0,0};
 651                 Double_t armenterosQtAlphaESD[2] = {0,0};
 652                 Double_t armenterosQtAlphaKFNew[2] = {0,0};
 653                 Double_t armenterosQtAlphaESDMC[2] = {0,0};
 654                 Double_t armenterosQtAlphaMC[2] = {0,0};
 655
 656                 GetArmenterosQtAlpha(GetNegativeKFParticle(), // old KF way calculating Qt Alpha
 657                                     GetPositiveKFParticle(),
 658                                     GetMotherCandidateKFCombination(),
 659                                     armenterosQtAlphaKF);
 660                 GetArmenterosQtAlpha(fCurrentV0,armenterosQtAlphaESD); // ESD way calculating Qt Alpha
 661                 GetArmenterosQtAlpha(GetNegativeKFParticle(), // new KF way calculating Qt Alpha
 662                                     GetPositiveKFParticle(),
 663                                     armenterosQtAlphaKFNew,fdoESDQtCut);
 664
 665                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinderKF_alfa_qt",armenterosQtAlphaKF[1],armenterosQtAlphaKF[0]);
 666                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinderESD_alfa_qt",armenterosQtAlphaESD[1],armenterosQtAlphaESD[0]);
 667                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinderKFNew_alfa_qt",armenterosQtAlphaKFNew[1],armenterosQtAlphaKFNew[0]);
 668
 669                 if(fdoESDQtCut == 0){
 670                    armenterosQtAlpha[0] = armenterosQtAlphaKF[0];
 671                    armenterosQtAlpha[1] = armenterosQtAlphaKF[1];
 672                 }
 673                 else if(fdoESDQtCut == 1){
 674                    armenterosQtAlpha[0] = armenterosQtAlphaESD[0];
 675                    armenterosQtAlpha[1] = armenterosQtAlphaESD[1];
 676                 }
 677
 678                 else if(fdoESDQtCut == 2 || fdoESDQtCut == 3){
 679                    armenterosQtAlpha[0] = armenterosQtAlphaKFNew[0];
 680                    armenterosQtAlpha[1] = armenterosQtAlphaKFNew[1];
 681                 }
 682
 683                 if(fCurrentNegativeESDTrack->Charge() == fCurrentPositiveESDTrack->Charge()){                                            // avoid like sign
 684                         //      iResult=kFALSE;
 685                         if(fHistograms != NULL ){
 686                                 fHistograms->FillHistogram("ESD_CutLikeSign_InvMass",GetMotherCandidateMass());
 687                                 fHistograms->FillHistogram("ESD_CutLikeSign_Pt",GetMotherCandidatePt());
 688                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 689                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 690                                 //      fUpdateV0AlreadyCalled = kTRUE;
 691                         }
 692                         fCurrentV0IndexNumber++;
 693                         continue;
 694                 }
 695                 if(fDoCF){
 696                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign);         // for CF
 697                 }
 698
 699                 if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ||
 700                         !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTPCrefit) ){
 701                         //      if( !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kITSrefit) ||
 702                         //                      !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kITSrefit) ){
 703                         //      iResult=kFALSE;
 704                         if(fHistograms != NULL){
 705                                 fHistograms->FillHistogram("ESD_CutRefit_InvMass",GetMotherCandidateMass());
 706                                 fHistograms->FillHistogram("ESD_CutRefit_Pt",GetMotherCandidatePt());
 707                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 708                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 709                                 //fUpdateV0AlreadyCalled = kTRUE;
 710                         }
 711                         fCurrentV0IndexNumber++;
 712                         continue;
 713                 }
 714                 if(fDoCF){
 715                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit);         // for CF
 716                 }
 717
 718
 719
 720                 if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 ||
 721                         fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {
 722                         //iResult=kFALSE;
 723                         if(fHistograms != NULL ){
 724                                 fHistograms->FillHistogram("ESD_CutKink_InvMass",GetMotherCandidateMass());
 725                                 fHistograms->FillHistogram("ESD_CutKink_Pt",GetMotherCandidatePt());
 726                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 727                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 728                                 //fUpdateV0AlreadyCalled = kTRUE;
 729                         }
 730                         fCurrentV0IndexNumber++;
 731                         continue;
 732                 }
 733
 734                 if(fDoCF){
 735                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks);            // for CF
 736                 }
 737
 738                 if(GetXYRadius()>fMaxR){ // cuts on distance from collision point
 739                         if(fHistograms != NULL){
 740                                 fHistograms->FillHistogram("ESD_CutR_InvMass",GetMotherCandidateMass());
 741                                 fHistograms->FillHistogram("ESD_CutR_Pt",GetMotherCandidatePt());
 742                         }
 743                         fCurrentV0IndexNumber++;
 744                         continue;
 745                 }
 746                 if(fDoCF){
 747                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepR);                        // for CF
 748                 }
 749                 if(GetXYRadius()<fMinR){ // cuts on distance from collision point
 750                         if(fHistograms != NULL){
 751                                 fHistograms->FillHistogram("ESD_CutMinR_InvMass",GetMotherCandidateMass());
 752                                 fHistograms->FillHistogram("ESD_CutMinR_Pt",GetMotherCandidatePt());
 753                         }
 754                         fCurrentV0IndexNumber++;
 755                         continue;
 756                 }
 757
 758                 //if((TMath::Abs(fCurrentZValue)*fLineCutZRSlope)-fLineCutZValue > GetXYRadius() ) { // cuts out regions where we do not reconstruct
 759                 if( GetXYRadius() <= ((TMath::Abs(fCurrentZValue)*fLineCutZRSlope)-fLineCutZValue)){
 760                         if(fHistograms != NULL){
 761                                 fHistograms->FillHistogram("ESD_CutLine_InvMass",GetMotherCandidateMass());
 762                                 fHistograms->FillHistogram("ESD_CutLine_Pt",GetMotherCandidatePt());
 763                         }
 764                         fCurrentV0IndexNumber++;
 765                         continue;
 766                 } else if (fUseEtaMinCut &&  GetXYRadius() >= ((TMath::Abs(fCurrentZValue)*fLineCutZRSlopeMin)-fLineCutZValueMin )){
 767                         if(fHistograms != NULL){
 768                                 fHistograms->FillHistogram("ESD_CutLine_InvMass",GetMotherCandidateMass());
 769                                 fHistograms->FillHistogram("ESD_CutLine_Pt",GetMotherCandidatePt());
 770                         }
 771                         fCurrentV0IndexNumber++;
 772                         continue;
 773                 }
 774
 775                 if(fDoCF){
 776                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepLine);                     // for CF
 777                 }
 778
 779                 if(TMath::Abs(fCurrentZValue) > fMaxZ ){ // cuts out regions where we do not reconstruct
 780                         if(fHistograms != NULL){
 781                                 fHistograms->FillHistogram("ESD_CutZ_InvMass",GetMotherCandidateMass());
 782                                 fHistograms->FillHistogram("ESD_CutZ_Pt",GetMotherCandidatePt());
 783                         }
 784                         fCurrentV0IndexNumber++;
 785                         continue;
 786                 }
 787                 if(fDoCF){
 788                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepZ);                // for CF
 789                 }
 790
 791                 if(fUseKFParticle){
 792                         if(TMath::Abs(fMotherCandidateLorentzVector->Eta())> fEtaCut || TMath::Abs(fMotherCandidateLorentzVector->Eta())< fEtaCutMin){
 793                                 if(fHistograms != NULL){
 794                                         fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
 795                                         fHistograms->FillHistogram("ESD_CutEta_Pt",GetMotherCandidatePt());
 796                                 }
 797                                 fCurrentV0IndexNumber++;
 798                                 continue;
 799                         }
 800
 801                         if(TMath::Abs(fCurrentNegativeKFParticle->GetEta())> fEtaCut || TMath::Abs(fCurrentNegativeKFParticle->GetEta())< fEtaCutMin){
 802                                 if(fHistograms != NULL){
 803                                         fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
 804                                         fHistograms->FillHistogram("ESD_CutEta_Pt",GetMotherCandidatePt());
 805                                 }
 806                                 fCurrentV0IndexNumber++;
 807                                 continue;
 808                         }
 809
 810                         if(TMath::Abs(fCurrentPositiveKFParticle->GetEta())> fEtaCut || TMath::Abs(fCurrentPositiveKFParticle->GetEta())< fEtaCutMin){
 811                                 if(fHistograms != NULL){
 812                                         fHistograms->FillHistogram("ESD_CutEta_InvMass",GetMotherCandidateMass());
 813                                         fHistograms->FillHistogram("ESD_CutEta_Pt",GetMotherCandidatePt());
 814                                 }
 815                                 fCurrentV0IndexNumber++;
 816                                 continue;
 817                         }
 818                 }
 819
 820                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_Pt_Qt",GetMotherCandidatePt(),armenterosQtAlpha[0]);
 821
 822                 if(fDoMC){
 823                         if ( HasSameMCMother() == kTRUE){
 824                                 GetArmenterosQtAlpha(fNegativeMCParticle,
 825                                         fPositiveMCParticle,
 826                                         fMotherMCParticle,
 827                                         armenterosQtAlphaMC);
 828                         }
 829                         GetArmenterosQtAlpha(fNegativeMCParticle,
 830                                 fPositiveMCParticle,
 831                                 GetMotherCandidateKFCombination(),
 832                                 armenterosQtAlphaESDMC );
 833                 }
 834
 835                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_goodtracks_alfa_qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 836                 if( fCurrentNegativeKFParticle->GetPt()> 0.150 &&       fCurrentPositiveKFParticle->GetPt()> 0.150){
 837                         fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_minPt_GT_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 838                 }
 839                 if(fDoMC){
 840                         fHistograms->FillHistogram("ESD_TrueConvAllV0s_ESDMother_Alpha_Qt",armenterosQtAlphaESDMC[1],armenterosQtAlphaESDMC[0]);
 841                         if ( HasSameMCMother() == kTRUE){
 842                                 fHistograms->FillHistogram("ESD_TrueConvSameMother_ESDMother_Alpha_Qt",armenterosQtAlphaESDMC[1],armenterosQtAlphaESDMC[0]);
 843                                 fHistograms->FillHistogram("ESD_TrueConvSameMother_MCMother_Alpha_Qt",armenterosQtAlphaMC[1],armenterosQtAlphaMC[0]);
 844                                 if (fMotherMCParticle->GetPdgCode() == 22 ){
 845                                         fHistograms->FillHistogram("ESD_TrueConvGamma_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 846                                         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Qt",GetMotherCandidatePt(),armenterosQtAlpha[0]);
 847                                 } else if ( fMotherMCParticle->GetPdgCode() == 310 ){
 848                                         fHistograms->FillHistogram("ESD_TrueConvK0s_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 849                                 } else if ( fMotherMCParticle->GetPdgCode() == 113 ){
 850                                         fHistograms->FillHistogram("ESD_TrueConvRho0_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 851                                 } else if ( fMotherMCParticle->GetPdgCode() == 333 ){
 852                                         fHistograms->FillHistogram("ESD_TrueConvPhi_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 853                                 } else if ( (fMotherMCParticle->GetPdgCode() == 3122 || fMotherMCParticle->GetPdgCode() == -3122) ){
 854                                         fHistograms->FillHistogram("ESD_TrueConvLambda_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 855                                 } else if ( (fMotherMCParticle->GetPdgCode() == 2114 || fMotherMCParticle->GetPdgCode() == -2114) ){
 856                                         fHistograms->FillHistogram("ESD_TrueConvDelta_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 857                                 } else if ( (fMotherMCParticle->GetPdgCode() == 313 ||
 858                                                                 fMotherMCParticle->GetPdgCode() == 323 ||
 859                                                                 fMotherMCParticle->GetPdgCode() == -323 ) ){
 860                                         fHistograms->FillHistogram("ESD_TrueConvKStar_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 861                                 } else {
 862                                         fHistograms->FillHistogram("ESD_TrueConvUnknown_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 863                                         fHistograms->FillHistogram("ESD_TrueConvUnknown_Qt_PDG",fMotherMCParticle->GetPdgCode());
 864                                 //      cout << "unidentfied mother: pdg-C mother " << fMotherMCParticle->GetPdgCode() << " daughters " << fNegativeMCParticle->GetPdgCode() << "\t" << fPositiveMCParticle->GetPdgCode() << endl;
 865                                 }
 866                         }       else {
 867                                 fHistograms->FillHistogram("ESD_TrueConvComb_Alpha_Qt",armenterosQtAlpha[1],armenterosQtAlpha[0]);
 868                         }
 869                 }
 870
 871                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_E_dEdxP",fCurrentNegativeESDTrack->P(),fCurrentNegativeESDTrack->GetTPCsignal());
 872                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_P_dEdxP",fCurrentPositiveESDTrack->P(),fCurrentPositiveESDTrack->GetTPCsignal());
 873                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_E_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron));
 874                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_P_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron));
 875                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_PiPl_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion));
 876                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_PiMi_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion));
 877                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_KPl_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kKaon));
 878                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_KMi_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kKaon));
 879                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_PPl_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kProton));
 880                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_PMi_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kProton));
 881                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_MuPl_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kMuon));
 882                 fHistograms->FillHistogram("ESD_AllV0sCurrentFinder_MuMi_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kMuon));
 883
 884                 if(fDodEdxSigmaCut == kTRUE){
 885                         if( fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
 886                         fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ||
 887                         fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
 888                         fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine ){
 889                                 //iResult=kFALSE;
 890                                 if(fHistograms != NULL ){
 891                                         fHistograms->FillHistogram("ESD_CutdEdxSigmaElectronLine_InvMass",GetMotherCandidateMass());
 892                                         fHistograms->FillHistogram("ESD_CutdEdxSigmaElectronLine_Pt",GetMotherCandidatePt());
 893                                         // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 894                                         // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 895                                         //fUpdateV0AlreadyCalled = kTRUE;
 896                                 }
 897                                 fCurrentV0IndexNumber++;
 898                                 continue;
 899                         }
 900                         if(fDoCF){
 901                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepdEdxElectronselection);                                                     // for CF
 902                         }
 903
 904                         fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_PiPl_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion));
 905                         fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_PiMi_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion));
 906                         fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_MuPl_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kMuon));
 907                         fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_MuMi_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kMuon));
 908
 909
 910                         if( fCurrentPositiveESDTrack->P()>fPIDMinPnSigmaAbovePionLine && fCurrentPositiveESDTrack->P()<fPIDMaxPnSigmaAbovePionLine ){
 911                                 if(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
 912                                         fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
 913                                         fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
 914                                         //              iResult=kFALSE;
 915                                         if(fHistograms != NULL){
 916                                                 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
 917                                                 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_Pt",GetMotherCandidatePt());
 918                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 919                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 920                                                 //fUpdateV0AlreadyCalled = kTRUE;
 921                                         }
 922                                         fCurrentV0IndexNumber++;
 923                                         continue;
 924                                 }
 925                         }
 926
 927                         if( fCurrentNegativeESDTrack->P()>fPIDMinPnSigmaAbovePionLine && fCurrentNegativeESDTrack->P()<fPIDMaxPnSigmaAbovePionLine){
 928                                 if(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
 929                                         fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
 930                                         fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
 931                                         //              iResult=kFALSE;
 932                                         if(fHistograms != NULL){
 933                                                 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
 934                                                 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_Pt",GetMotherCandidatePt());
 935                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 936                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 937                                                 //fUpdateV0AlreadyCalled = kTRUE;
 938                                         }
 939                                         fCurrentV0IndexNumber++;
 940                                         continue;
 941                                 }
 942                         }
 943
 944                         // High Pt
 945                         if( fCurrentPositiveESDTrack->P()>fPIDMaxPnSigmaAbovePionLine ){
 946                                 if(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
 947                                         fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
 948                                         fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLineHighPt){
 949                                         //              iResult=kFALSE;
 950                                         if(fHistograms != NULL){
 951                                                 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
 952                                                 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_Pt",GetMotherCandidatePt());
 953                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 954                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 955                                                 //fUpdateV0AlreadyCalled = kTRUE;
 956                                         }
 957                                         fCurrentV0IndexNumber++;
 958                                         continue;
 959                                 }
 960                         }
 961
 962                         if( fCurrentNegativeESDTrack->P()>fPIDMaxPnSigmaAbovePionLine){
 963                                 if(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
 964                                         fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
 965                                         fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion)<fPIDnSigmaAbovePionLineHighPt){
 966                         //              iResult=kFALSE;
 967                                         if(fHistograms != NULL){
 968                                                 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_InvMass",GetMotherCandidateMass());
 969                                                 fHistograms->FillHistogram("ESD_CutdEdxSigmaPionLine_Pt",GetMotherCandidatePt());
 970                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 971                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 972                                                 //fUpdateV0AlreadyCalled = kTRUE;
 973                                         }
 974                                         fCurrentV0IndexNumber++;
 975                                         continue;
 976                                 }
 977                         }
 978
 979
 980                         if(fDoCF){
 981                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepdEdxPionrejection);                                                         // for CF
 982                         }
 983
 984                 }
 985
 986                 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_KPl_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kKaon));
 987                 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_KMi_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kKaon));
 988
 989                 if(fDoKaonRejectionLowP == kTRUE){
 990                         if( fCurrentNegativeESDTrack->P()<fPIDMinPKaonRejectionLowP ){
 991                                 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kKaon))<fPIDnSigmaAtLowPAroundKaonLine){
 992                                         if(fHistograms != NULL){
 993                                                 fHistograms->FillHistogram("ESD_CutKaonRejectionLowP_InvMass",GetMotherCandidateMass());
 994                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
 995                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
 996                                                 //fUpdateV0AlreadyCalled = kTRUE;
 997                                         }
 998                                         fCurrentV0IndexNumber++;
 999                                         continue;
1000                                 }
1001                         }
1002                         if( fCurrentPositiveESDTrack->P()<fPIDMinPKaonRejectionLowP ){
1003                                 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kKaon))<fPIDnSigmaAtLowPAroundKaonLine){
1004                                         if(fHistograms != NULL){
1005                                                 fHistograms->FillHistogram("ESD_CutKaonRejectionLowP_InvMass",GetMotherCandidateMass());
1006                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1007                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1008                                                 //fUpdateV0AlreadyCalled = kTRUE;
1009                                         }
1010                                         fCurrentV0IndexNumber++;
1011                                         continue;
1012                                 }
1013                         }
1014                 }
1015
1016                 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_PPl_SigdEdxP",fCurrentPositiveESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kProton));
1017                 fHistograms->FillHistogram("ESD_ConvGammaBeforeCorresCut_PMi_SigdEdxP",fCurrentNegativeESDTrack->P(),fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kKaon));
1018
1019                 if(fDoProtonRejectionLowP == kTRUE){
1020                         if( fCurrentNegativeESDTrack->P()<fPIDMinPProtonRejectionLowP){
1021                                 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kProton))<fPIDnSigmaAtLowPAroundProtonLine){
1022                                         if(fHistograms != NULL){
1023                                                 fHistograms->FillHistogram("ESD_CutProtonRejectionLowP_InvMass",GetMotherCandidateMass());
1024                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1025                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1026                                                 //fUpdateV0AlreadyCalled = kTRUE;
1027                                         }
1028                                         fCurrentV0IndexNumber++;
1029                                         continue;
1030                                 }
1031                         }
1032                         if( fCurrentPositiveESDTrack->P()<fPIDMinPProtonRejectionLowP ){
1033                                 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kProton))<fPIDnSigmaAtLowPAroundProtonLine){
1034                                         if(fHistograms != NULL){
1035                                                 fHistograms->FillHistogram("ESD_CutProtonRejectionLowP_InvMass",GetMotherCandidateMass());
1036                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1037                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1038                                                 //fUpdateV0AlreadyCalled = kTRUE;
1039                                         }
1040                                         fCurrentV0IndexNumber++;
1041                                         continue;
1042                                 }
1043                         }
1044                 }
1045
1046                 if(fDoPionRejectionLowP == kTRUE){
1047                         if( fCurrentNegativeESDTrack->P()<fPIDMinPPionRejectionLowP ){
1048                                 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentNegativeESDTrack,AliPID::kPion))<fPIDnSigmaAtLowPAroundPionLine){
1049                                         if(fHistograms != NULL){
1050                                                 fHistograms->FillHistogram("ESD_CutPionRejectionLowP_InvMass",GetMotherCandidateMass());
1051                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1052                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1053                                                 //fUpdateV0AlreadyCalled = kTRUE;
1054                                         }
1055                                         fCurrentV0IndexNumber++;
1056                                         continue;
1057                                 }
1058                         }
1059                         if( fCurrentPositiveESDTrack->P()<fPIDMinPPionRejectionLowP ){
1060                                 if( TMath::Abs(fgESDpid->NumberOfSigmasTPC(fCurrentPositiveESDTrack,AliPID::kPion))<fPIDnSigmaAtLowPAroundPionLine){
1061                                         if(fHistograms != NULL){
1062                                                 fHistograms->FillHistogram("ESD_CutPionRejectionLowP_InvMass",GetMotherCandidateMass());
1063                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1064                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1065                                                 //fUpdateV0AlreadyCalled = kTRUE;
1066                                         }
1067                                         fCurrentV0IndexNumber++;
1068                                         continue;
1069                                 }
1070                         }
1071                 }
1072
1073
1074                 Double_t psiPair  = -1;
1075                 psiPair = GetPsiPair(fCurrentV0);
1076
1077                 if(psiPair > fPsiPairCut){
1078                    if(fHistograms != NULL ){
1079                       fHistograms->FillHistogram("ESD_CutPsiPair_InvMass",GetMotherCandidateMass());
1080                       fHistograms->FillHistogram("ESD_CutPsiPair_Pt",GetMotherCandidatePt());
1081                       // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1082                       // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1083                       //        fUpdateV0AlreadyCalled = kTRUE;
1084                    }
1085
1086                    fCurrentV0IndexNumber++;
1087                    continue;
1088                 }
1089
1090
1091                 Double_t cosineOfPointingAngle  = -1;
1092                 cosineOfPointingAngle = GetV0CosineOfPointingAngle(fCurrentXValue,fCurrentYValue,fCurrentZValue);
1093                 if(cosineOfPointingAngle < fCosinePointCut){
1094                    if(fHistograms != NULL ){
1095                         fHistograms->FillHistogram("ESD_CutCosinePoint_InvMass",GetMotherCandidateMass());
1096                         fHistograms->FillHistogram("ESD_CutCosinePoint_Pt",GetMotherCandidatePt());
1097                         // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1098                         // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1099                         //      fUpdateV0AlreadyCalled = kTRUE;
1100                         }
1101                         fCurrentV0IndexNumber++;
1102                         continue;
1103                 }
1104
1105
1106                 if( fDoTOFsigmaCut == kTRUE ){ // RRnewTOF start /////////////////////////////////////////////////////////////////////////////
1107                         Bool_t PosTrackNotTOFelec = kFALSE;
1108                         Bool_t NegTrackNotTOFelec = kFALSE;
1109                         if( (fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTOFpid) && !(fCurrentPositiveESDTrack->GetStatus() & AliESDtrack::kTOFmismatch) ){
1110                                 Double_t t0pos = fgESDpid->GetTOFResponse().GetStartTime(fCurrentPositiveESDTrack->P());
1111                                 Double_t fnSigmaPos = fgESDpid->NumberOfSigmasTOF(fCurrentPositiveESDTrack, AliPID::kElectron, t0pos);
1112                                 if( (fnSigmaPos>fTofPIDnSigmaAboveElectronLine) || (fnSigmaPos<fTofPIDnSigmaBelowElectronLine) ) PosTrackNotTOFelec = kTRUE;
1113                         }
1114                         if( (fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTOFpid) && !(fCurrentNegativeESDTrack->GetStatus() & AliESDtrack::kTOFmismatch) ){
1115                                 Double_t t0neg = fgESDpid->GetTOFResponse().GetStartTime(fCurrentNegativeESDTrack->P());
1116                                 Double_t fnSigmaNeg = fgESDpid->NumberOfSigmasTOF(fCurrentNegativeESDTrack, AliPID::kElectron, t0neg);
1117                                 if( (fnSigmaNeg>fTofPIDnSigmaAboveElectronLine) || (fnSigmaNeg<fTofPIDnSigmaBelowElectronLine) ) NegTrackNotTOFelec = kTRUE;
1118                         }
1119                         if( (PosTrackNotTOFelec==kTRUE) || (NegTrackNotTOFelec==kTRUE) ){
1120                                 if(fHistograms != NULL){
1121                                         fHistograms->FillHistogram("ESD_CutTOFsigmaElec_InvMass",GetMotherCandidateMass());
1122                                         fHistograms->FillHistogram("ESD_CutTOFsigmaElec_Pt",GetMotherCandidatePt());
1123                                 }
1124                                 fCurrentV0IndexNumber++;
1125                                 continue;
1126                         }
1127                 } /////////////////////////////// RRnewTOF end ///////////////////////////////////////////////////////////////////////////////
1128
1129
1130                 // Gamma selection based on QT from Armenteros
1131                 if(fDoQtGammaSelection == kTRUE){ // RRnew start : apply different qT-cut above/below
1132                         if(fDoHighPtQtGammaSelection){
1133                                 if(GetMotherCandidatePt() < fPtBorderForQt){
1134                                         if(armenterosQtAlpha[0]>fQtMax){
1135                                                 if(fHistograms != NULL){
1136                                                         fHistograms->FillHistogram("ESD_CutQt_InvMass",GetMotherCandidateMass());
1137                                                         fHistograms->FillHistogram("ESD_CutQt_Pt",GetMotherCandidatePt());
1138                                                 }
1139                                                 fCurrentV0IndexNumber++;
1140                                                 continue;
1141                                         }
1142                                 } else {
1143                                         if(armenterosQtAlpha[0]>fHighPtQtMax)   {
1144                                                 if(fHistograms != NULL){
1145                                                         fHistograms->FillHistogram("ESD_CutQt_InvMass",GetMotherCandidateMass());
1146                                                         fHistograms->FillHistogram("ESD_CutQt_Pt",GetMotherCandidatePt());
1147                                                 }
1148                                                 fCurrentV0IndexNumber++;
1149                                                 continue;
1150                                         }
1151                                 }
1152                         } else {
1153                                 if(armenterosQtAlpha[0]>fQtMax){
1154                                         if(fHistograms != NULL){
1155                                                 fHistograms->FillHistogram("ESD_CutQt_InvMass",GetMotherCandidateMass());
1156                                                 fHistograms->FillHistogram("ESD_CutQt_Pt",GetMotherCandidatePt());
1157                                         }
1158                                         fCurrentV0IndexNumber++;
1159                                         continue;
1160                                 }
1161                         }
1162                 }        // RRnew end
1163
1164                 if(fDoPhotonAsymmetryCut == kTRUE){
1165                         if( fNegativeTrackLorentzVector->P()>fMinPPhotonAsymmetryCut ){
1166                                 Double_t trackNegAsy=0;
1167                                 if (fCurrentMotherKFCandidate->GetP()!=0.){
1168                                         trackNegAsy= fNegativeTrackLorentzVector->P()/fMotherCandidateLorentzVector->P();
1169                                 }
1170                                 if( trackNegAsy<fMinPhotonAsymmetry ||trackNegAsy>(1.- fMinPhotonAsymmetry)){
1171                                         if(fHistograms != NULL){
1172                                                 fHistograms->FillHistogram("ESD_CutPhotonAsymmetry_InvMass",GetMotherCandidateMass());
1173                                                 fHistograms->FillHistogram("ESD_CutPhotonAsymmetry_Pt",GetMotherCandidatePt());
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                                         fCurrentV0IndexNumber++;
1179                                         continue;
1180                                 }
1181                         }
1182
1183                         if( fPositiveTrackLorentzVector->P()>fMinPPhotonAsymmetryCut ){
1184                                 Double_t trackPosAsy=0;
1185                                 if (fCurrentMotherKFCandidate->GetP()!=0.){
1186                                         trackPosAsy= fPositiveTrackLorentzVector->P()/fMotherCandidateLorentzVector->P();
1187                                 }
1188                                 if( trackPosAsy<fMinPhotonAsymmetry ||trackPosAsy>(1.- fMinPhotonAsymmetry)){
1189                                         if(fHistograms != NULL){
1190                                                 fHistograms->FillHistogram("ESD_CutPhotonAsymmetry_InvMass",GetMotherCandidateMass());
1191                                                 fHistograms->FillHistogram("ESD_CutPhotonAsymmetry_Pt",GetMotherCandidatePt());
1192                                                 // to avoid filling the other cut histograms. So in this case fUpdateV0AlreadyCalled also serves as a flag for the histogram filling
1193                                                 // it will anyway be set to true at the end of the UpdateV0Information function, and there are no return until the end
1194                                                 //fUpdateV0AlreadyCalled = kTRUE;
1195                                         }
1196                                         fCurrentV0IndexNumber++;
1197                                         continue;
1198                                 }
1199                         }
1200                 }
1201                 //checks if we have a prim vertex
1202                 //if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {
1203                 if(GetNumberOfContributorsVtx()<=0) {
1204                         if(fHistograms != NULL){
1205                                 fHistograms->FillHistogram("ESD_CutNContributors_InvMass",GetMotherCandidateMass());
1206                                 fHistograms->FillHistogram("ESD_CutNContributors_Pt",GetMotherCandidatePt());
1207                         }
1208                         fCurrentV0IndexNumber++;
1209                         continue;
1210                 }
1211                 if(fDoCF){
1212                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepNContributors);            // for CF
1213                 }
1214
1215                 //Check the pid probability
1216                 if(CheckPIDProbability(fPIDProbabilityCutNegativeParticle,fPIDProbabilityCutPositiveParticle)==kFALSE){
1217                         if(fHistograms != NULL){
1218                                 fHistograms->FillHistogram("ESD_CutPIDProb_InvMass",GetMotherCandidateMass());
1219                                 fHistograms->FillHistogram("ESD_CutPIDProb_Pt",GetMotherCandidatePt());
1220                         }
1221                         fCurrentV0IndexNumber++;
1222                         continue;
1223                 }
1224                 if(fDoCF){
1225                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCPID);                   // for CF
1226                 }
1227
1228
1229                 /* Moved further up so corr framework can work
1230                          if(UpdateV0Information() == kFALSE){
1231                          fCurrentV0IndexNumber++;
1232                          continue;
1233                          }
1234                 */
1235                 if(fCurrentNegativeESDTrack->GetNcls(1) < fMinClsTPC || fCurrentPositiveESDTrack->GetNcls(1) < fMinClsTPC ){
1236                         if(fHistograms != NULL){
1237                                 fHistograms->FillHistogram("ESD_CutMinNClsTPC_InvMass",GetMotherCandidateMass());
1238                                 fHistograms->FillHistogram("ESD_CutMinNClsTPC_Pt",GetMotherCandidatePt());
1239                         }
1240                         fCurrentV0IndexNumber++;
1241                         continue;
1242                 }
1243                 if(fDoCF){
1244                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepMinClsTPC);                // for CF
1245                 }
1246                 Double_t negclsToF = 0.;
1247                 if (!fUseCorrectedTPCClsInfo ){
1248                         if(fCurrentNegativeESDTrack->GetTPCNclsF()!=0   ){
1249                                 negclsToF = (Double_t)fCurrentNegativeESDTrack->GetNcls(1)/(Double_t)fCurrentNegativeESDTrack->GetTPCNclsF();
1250                         }
1251                 } else {
1252                         negclsToF = fCurrentNegativeESDTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
1253                 }
1254
1255                 Double_t posclsToF = 0.;
1256                 if (!fUseCorrectedTPCClsInfo ){
1257                         if(fCurrentPositiveESDTrack->GetTPCNclsF()!=0   ){
1258                                 posclsToF = (Double_t)fCurrentPositiveESDTrack->GetNcls(1)/(Double_t)fCurrentPositiveESDTrack->GetTPCNclsF();
1259                         }
1260                 }else{
1261                         posclsToF = fCurrentPositiveESDTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
1262                 }
1263
1264                 if( negclsToF < fMinClsTPCToF ||        posclsToF < fMinClsTPCToF ){
1265                         if(fHistograms != NULL){
1266                                 fHistograms->FillHistogram("ESD_CutMinNClsTPCToF_InvMass",GetMotherCandidateMass());
1267                                 fHistograms->FillHistogram("ESD_CutMinNClsTPCToF_Pt",GetMotherCandidatePt());
1268                         }
1269                         fCurrentV0IndexNumber++;
1270                         continue;
1271                 }
1272
1273
1274
1275
1276                 if(fUseKFParticle){
1277
1278
1279                         if( fCurrentNegativeKFParticle->GetPt()< fSinglePtCut ||        fCurrentPositiveKFParticle->GetPt()< fSinglePtCut){
1280                                 if(fHistograms != NULL){
1281                                         fHistograms->FillHistogram("ESD_CutSinglePt_InvMass",GetMotherCandidateMass());
1282                                         fHistograms->FillHistogram("ESD_CutSinglePt_Pt",GetMotherCandidatePt());
1283                                 }
1284                                 fCurrentV0IndexNumber++;
1285                                 continue;
1286                         }
1287                         if(fDoCF){
1288                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepSinglePt);         // for CF
1289                         }
1290
1291
1292                         if(fCurrentMotherKFCandidate->GetNDF()<=0){
1293                                 if(fHistograms != NULL){
1294                                         fHistograms->FillHistogram("ESD_CutNDF_InvMass",GetMotherCandidateMass());
1295                                         fHistograms->FillHistogram("ESD_CutNDF_Pt",GetMotherCandidatePt());
1296                                 }
1297                                 fCurrentV0IndexNumber++;
1298                                 continue;
1299                         }
1300                         if(fDoCF){
1301                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepNDF);              // for CF
1302                         }
1303
1304                         Double_t chi2V0 = fCurrentMotherKFCandidate->GetChi2()/fCurrentMotherKFCandidate->GetNDF();
1305                         if(chi2V0 > fChi2CutConversion || chi2V0 <=0){
1306                                 if(fHistograms != NULL){
1307                                         fHistograms->FillHistogram("ESD_CutChi2_InvMass",GetMotherCandidateMass());
1308                                         fHistograms->FillHistogram("ESD_CutChi2_Pt",GetMotherCandidatePt());
1309                                 }
1310                                 fCurrentV0IndexNumber++;
1311                                 continue;
1312                         }
1313                         if(fDoCF){
1314                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepChi2);                     // for CF
1315                         }
1316
1317                         if(fDoCF){
1318                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepEta);                      // for CF
1319                         }
1320
1321                         if(fMotherCandidateLorentzVector->Pt()<fPtCut){
1322                                 if(fHistograms != NULL){
1323                                         fHistograms->FillHistogram("ESD_CutPt_InvMass",GetMotherCandidateMass());
1324                                         fHistograms->FillHistogram("ESD_CutPt_Pt",GetMotherCandidatePt());
1325                                 }
1326                                 fCurrentV0IndexNumber++;
1327                                 continue;
1328                         }
1329                         if(fDoCF){
1330                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepPt);                       // for CF
1331                         }
1332
1333                 }
1334                 else if(fUseESDTrack){
1335                         //TODO
1336                 }
1337
1338                 if(fHistograms != NULL){
1339                         fHistograms->FillHistogram("ESD_GoodV0s_InvMass",GetMotherCandidateMass());
1340                 }
1341
1342                 //              fCurrentEventGoodV0s.push_back(*fCurrentMotherKFCandidate);
1343
1344                 if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
1345                         fCurrentMotherKFCandidate->E()=fCurrentMotherKFCandidate->GetP();
1346                 }
1347
1348                 if(fDoMC&& fUseMCPSmearing>0){
1349                         SmearKFParticle(fCurrentMotherKFCandidate);
1350
1351                 }
1352
1353                 new((*fCurrentEventGoodV0s)[fCurrentEventGoodV0s->GetEntriesFast()])    AliKFParticle(*fCurrentMotherKFCandidate);
1354                 fV0Pindex.push_back(fCurrentV0->GetPindex());
1355                 fV0Nindex.push_back(fCurrentV0->GetNindex());
1356
1357                 iResult=kTRUE;//means we have a v0 who survived all the cuts applied
1358
1359                 fNumberOfGoodV0s++;
1360
1361                 fCurrentV0IndexNumber++;
1362
1363                 break;
1364         }
1365         return iResult;
1366 }
1367
1368 Bool_t AliV0Reader::UpdateV0Information(){
1369         //see header file for documentation
1370
1371         const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0);
1372         const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0);
1373
1374         Bool_t iResult=kTRUE;                                           // for taking out not refitted, kinks and like sign tracks
1375
1376         Bool_t switchTracks = kFALSE;
1377
1378         fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1379         fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1380
1381
1382         if(fCurrentPositiveESDTrack->GetSign() == -1 && fCurrentNegativeESDTrack->GetSign() == 1){      // switch wrong signed tracks
1383                 fCurrentNegativeESDTrack = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1384                 fCurrentPositiveESDTrack = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1385                 switchTracks = kTRUE;
1386         }
1387
1388
1389         if(fCurrentNegativeKFParticle != NULL){
1390            delete fCurrentNegativeKFParticle;
1391         }
1392         if(switchTracks == kFALSE){
1393                 fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fNegativeTrackPID);
1394         }
1395         else{
1396                 fCurrentNegativeKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fNegativeTrackPID);
1397         }
1398
1399         if(fCurrentPositiveKFParticle != NULL){
1400                 delete fCurrentPositiveKFParticle;
1401         }
1402         if(switchTracks == kFALSE){
1403                 fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamP()),fPositiveTrackPID);
1404         }
1405         else{
1406                 fCurrentPositiveKFParticle = new AliKFParticle(*(fCurrentV0->GetParamN()),fPositiveTrackPID);
1407         }
1408
1409         if(fCurrentMotherKFCandidate != NULL){
1410                 delete fCurrentMotherKFCandidate;
1411         }
1412
1413         if(fUseConstructGamma==kTRUE){
1414                 fCurrentMotherKFCandidate = new AliKFParticle();//(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
1415                 fCurrentMotherKFCandidate->ConstructGamma(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
1416         }else{
1417                 fCurrentMotherKFCandidate = new AliKFParticle(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
1418                 if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
1419                         fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);
1420                 }
1421         }
1422         if(fUseImprovedVertex == kTRUE){
1423                 AliKFVertex primaryVertexImproved(*GetPrimaryVertex());
1424                 primaryVertexImproved+=*fCurrentMotherKFCandidate;
1425                 fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);
1426         }
1427
1428         fCurrentMotherKFCandidate->GetMass(fMotherCandidateKFMass,fMotherCandidateKFWidth);
1429
1430         if(fNegativeTrackLorentzVector != NULL){
1431                 delete fNegativeTrackLorentzVector;
1432         }
1433         if(fUseKFParticle){
1434                 fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeKFParticle->Px(),fCurrentNegativeKFParticle->Py(),fCurrentNegativeKFParticle->Pz());
1435         }
1436         else    { //if(fUseESDTrack){
1437                 fNegativeTrackLorentzVector = new TLorentzVector(fCurrentNegativeESDTrack->Px(),fCurrentNegativeESDTrack->Py(),fCurrentNegativeESDTrack->Pz());
1438         }
1439
1440         if(fPositiveTrackLorentzVector != NULL){
1441                 delete fPositiveTrackLorentzVector;
1442         }
1443         if(fUseKFParticle){
1444                 fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveKFParticle->Px(),fCurrentPositiveKFParticle->Py(),fCurrentPositiveKFParticle->Pz());
1445         }
1446         else    { // if(fUseESDTrack){  fPositiveTrackLorentzVector must be reinitialized, so assuming use ESD if not kfparticle Svein.
1447                 fPositiveTrackLorentzVector = new TLorentzVector(fCurrentPositiveESDTrack->Px(),fCurrentPositiveESDTrack->Py(),fCurrentPositiveESDTrack->Pz());
1448         }
1449
1450         if(fMotherCandidateLorentzVector != NULL){
1451                 delete fMotherCandidateLorentzVector;
1452         }
1453
1454         fMotherCandidateLorentzVector = new TLorentzVector(*fNegativeTrackLorentzVector + *fPositiveTrackLorentzVector);
1455
1456         if(fPositiveTrackPID==-11 && fNegativeTrackPID==11){
1457                 fMotherCandidateLorentzVector->SetXYZM(fMotherCandidateLorentzVector->Px() ,fMotherCandidateLorentzVector->Py(),fMotherCandidateLorentzVector->Pz(),0.);
1458         }
1459
1460
1461         if(fDoMC == kTRUE){
1462                 fMotherMCParticle= NULL;
1463                 if(switchTracks == kFALSE){
1464                         fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
1465                         fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
1466                 }else{
1467                         fNegativeMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetPindex())->GetLabel()));
1468                         fPositiveMCParticle = fMCStack->Particle(TMath::Abs(fESDEvent->GetTrack(fCurrentV0->GetNindex())->GetLabel()));
1469                 }
1470
1471                 if(fPositiveMCParticle->GetMother(0)>-1){
1472                         fMotherMCParticle = fMCStack->Particle(fPositiveMCParticle->GetMother(0));
1473                 }
1474         }
1475
1476
1477
1478
1479
1480         // for CF
1481 //       Double_t containerInput[3];
1482 //       if(fDoCF){
1483 //               containerInput[0] = GetMotherCandidatePt();
1484 //               containerInput[1] = GetMotherCandidateEta();
1485 //               containerInput[2] = GetMotherCandidateMass();
1486
1487 //               fCFManager->GetParticleContainer()->Fill(containerInput,kStepLikeSign);                // for CF
1488 //               fCFManager->GetParticleContainer()->Fill(containerInput,kStepTPCRefit);                // for CF
1489 //               fCFManager->GetParticleContainer()->Fill(containerInput,kStepKinks);           // for CF
1490 //       }
1491
1492
1493         if(fUseOwnXYZCalculation == kFALSE){
1494                 if(fUseConstructGamma == kFALSE){
1495                         fCurrentV0->GetXYZ(fCurrentXValue,fCurrentYValue,fCurrentZValue);
1496                 }else{
1497                         fCurrentXValue=GetMotherCandidateKFCombination()->GetX();
1498                         fCurrentYValue=GetMotherCandidateKFCombination()->GetY();
1499                         fCurrentZValue=GetMotherCandidateKFCombination()->GetZ();
1500                 }
1501         }
1502         else{
1503                 Double_t convpos[3]={0,0,0};
1504                 GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos);
1505 //              fCurrentMotherKF->SetConversionPoint(convpos);
1506
1507 //              Double_t convpos[2];
1508 //              convpos[0]=0;
1509 //              convpos[1]=0;
1510 //
1511 //              GetConvPosXY(GetPositiveESDTrack(),GetNegativeESDTrack(),GetMagneticField(),convpos);
1512 //
1513                 fCurrentXValue = convpos[0];
1514                 fCurrentYValue = convpos[1];
1515                 fCurrentZValue = convpos[2];
1516         }
1517         fUpdateV0AlreadyCalled = kTRUE;
1518
1519         return iResult;
1520 }
1521
1522
1523
1524 Bool_t AliV0Reader::HasSameMCMother(){
1525         //see header file for documentation
1526
1527         Bool_t iResult = kFALSE;
1528         if(fDoMC == kTRUE){
1529                 if(fNegativeMCParticle != NULL && fPositiveMCParticle != NULL){
1530                         if(fNegativeMCParticle->GetMother(0) == fPositiveMCParticle->GetMother(0))
1531         if(fMotherMCParticle){
1532                 iResult = kTRUE;
1533         }
1534                 }
1535         }
1536         return iResult;
1537 }
1538
1539 Bool_t AliV0Reader::CheckPIDProbability(Double_t negProbCut, Double_t posProbCut){
1540         //see header file for documentation
1541
1542         Bool_t iResult=kFALSE;
1543
1544         //      Double_t *posProbArray = new Double_t[10];
1545         //      Double_t *negProbArray = new Double_t[10];
1546         //-AM The TPCpid method expects an array of length kSPECIES that is 5 not 10
1547
1548         Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1549         Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1550
1551         AliESDtrack* negTrack   = GetNegativeESDTrack();
1552         AliESDtrack* posTrack   = GetPositiveESDTrack();
1553         //fESDEvent->GetTrack(fCurrentV0->GetNindex());
1554                 //fESDEvent->GetTrack(fCurrentV0->GetPindex());
1555         //-AM for switchtracks==true the above is a bug
1556
1557         if(negProbArray && posProbArray){
1558
1559                 negTrack->GetTPCpid(negProbArray);
1560                 posTrack->GetTPCpid(posProbArray);
1561
1562                 //      if(negProbArray != NULL && posProbArray != NULL){ // this is not allowed anymore for some reason(RC19)
1563                 if(negProbArray[GetSpeciesIndex(-1)]>=negProbCut && posProbArray[GetSpeciesIndex(1)]>=posProbCut){
1564                         iResult=kTRUE;
1565                 }
1566         }
1567         delete [] posProbArray;
1568         delete [] negProbArray;
1569         return iResult;
1570 }
1571
1572 void AliV0Reader::GetPIDProbability(Double_t &negPIDProb,Double_t & posPIDProb){
1573         // see header file for documentation
1574
1575         //Double_t *posProbArray = new Double_t[10];
1576         // Double_t *negProbArray = new Double_t[10];
1577         //-AM The TPCpid method expects an array of length kSPECIES that is 5 not 10
1578         Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1579         Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1580
1581 //       AliESDtrack* negTrack  = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1582 //       AliESDtrack* posTrack  = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1583         //-AM for switchtracks the above is a bug
1584         AliESDtrack* negTrack   = GetNegativeESDTrack();
1585         AliESDtrack* posTrack   = GetPositiveESDTrack();
1586
1587         if(negProbArray && posProbArray){
1588                 negTrack->GetTPCpid(negProbArray);
1589                 posTrack->GetTPCpid(posProbArray);
1590
1591                 //      if(negProbArray!=NULL && posProbArray!=NULL){ // this is not allowed anymore for some reason(RC19)
1592                 negPIDProb = negProbArray[GetSpeciesIndex(-1)];
1593                 posPIDProb = posProbArray[GetSpeciesIndex(1)];
1594         }
1595         delete [] posProbArray;
1596         delete [] negProbArray;
1597 }
1598
1599 void AliV0Reader::GetPIDProbabilityMuonPion(Double_t &negPIDProb,Double_t & posPIDProb){
1600         // see header file for documentation
1601
1602
1603         Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1604         Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1605
1606         // AliESDtrack* negTrack        = fESDEvent->GetTrack(fCurrentV0->GetNindex());
1607         // AliESDtrack* posTrack        = fESDEvent->GetTrack(fCurrentV0->GetPindex());
1608         //-AM for switchtracks the above is a bug
1609
1610         AliESDtrack* negTrack   = GetNegativeESDTrack();
1611         AliESDtrack* posTrack   = GetPositiveESDTrack();
1612
1613         if(negProbArray && posProbArray){
1614                 negTrack->GetTPCpid(negProbArray);
1615                 posTrack->GetTPCpid(posProbArray);
1616
1617                 //      if(negProbArray!=NULL && posProbArray!=NULL){ // this is not allowed anymore for some reason(RC19)
1618
1619                 negPIDProb = negProbArray[1]+negProbArray[2];
1620                 posPIDProb = posProbArray[1]+posProbArray[2];
1621         }
1622         delete [] posProbArray;
1623         delete [] negProbArray;
1624 }
1625
1626 void AliV0Reader::UpdateEventByEventData(){
1627         //see header file for documentation
1628         if(fCurrentEventGoodV0s->GetEntriesFast() >0 ){
1629                 if(fCalculateBackground){
1630                         if(fUseChargedTrackMultiplicityForBG == kTRUE){
1631         fBGEventHandler->AddEvent(fCurrentEventGoodV0s,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
1632         //filling z and multiplicity histograms
1633         fHistograms->FillHistogram("ESD_Z_distribution",fESDEvent->GetPrimaryVertex()->GetZ());
1634         fHistograms->FillHistogram("ESD_multiplicity_distribution",CountESDTracks());
1635         fHistograms->FillHistogram("ESD_ZvsMultiplicity",fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
1636                         }
1637                         else{ // means we use #V0s for multiplicity
1638         fBGEventHandler->AddEvent(fCurrentEventGoodV0s,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
1639         //filling z and multiplicity histograms
1640         fHistograms->FillHistogram("ESD_Z_distribution",fESDEvent->GetPrimaryVertex()->GetZ());
1641         fHistograms->FillHistogram("ESD_multiplicity_distribution",fNumberOfGoodV0s);
1642         fHistograms->FillHistogram("ESD_ZvsMultiplicity",fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
1643                         }
1644                 }
1645         }
1646         fCurrentEventGoodV0s->Delete();
1647         fCurrentV0IndexNumber=0;
1648         fNumberOfESDTracks=0;
1649
1650         fV0Pindex.clear();
1651         fV0Nindex.clear();
1652
1653         delete fCurrentEventGoodV0s;
1654         fCurrentEventGoodV0s = NULL;
1655         delete fBrem;
1656         fBrem = NULL;
1657         delete fCurrentNegativeKFParticle;
1658         fCurrentNegativeKFParticle = NULL;
1659         delete fCurrentPositiveKFParticle;
1660         fCurrentPositiveKFParticle = NULL;
1661         delete fCurrentMotherKFCandidate;
1662         fCurrentMotherKFCandidate = NULL;
1663         delete fNegativeTrackLorentzVector;
1664         fNegativeTrackLorentzVector = NULL;
1665         delete fPositiveTrackLorentzVector;
1666         fPositiveTrackLorentzVector = NULL;
1667         delete fMotherCandidateLorentzVector;
1668         fMotherCandidateLorentzVector = NULL;
1669
1670         //      fBGEventHandler->PrintBGArray(); // for debugging
1671 }
1672
1673
1674 Double_t AliV0Reader::GetNegativeTrackPhi() const{
1675         //see header file for documentation
1676
1677         Double_t offset=0;
1678         if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){
1679                 offset = -2*TMath::Pi();
1680         }
1681         return fNegativeTrackLorentzVector->Phi()+offset;
1682 }
1683
1684 Double_t AliV0Reader::GetPositiveTrackPhi() const{
1685         //see header file for documentation
1686
1687         Double_t offset=0;
1688         if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){
1689                 offset = -2*TMath::Pi();
1690         }
1691         return fPositiveTrackLorentzVector->Phi()+offset;
1692 }
1693
1694 Double_t AliV0Reader::GetMotherCandidatePhi() const{
1695         //see header file for documentation
1696
1697         Double_t offset=0;
1698         if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){
1699                 offset = -2*TMath::Pi();
1700         }
1701         return fMotherCandidateLorentzVector->Phi()+offset;
1702 }
1703
1704
1705 Double_t AliV0Reader::GetMotherCandidateRapidity() const{
1706         //see header file for documentation
1707
1708         Double_t rapidity=0;
1709         if(fMotherCandidateLorentzVector->Energy() - fMotherCandidateLorentzVector->Pz() == 0 || fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz() == 0) rapidity=0;
1710         else rapidity = 0.5*(TMath::Log((fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz()) / (fMotherCandidateLorentzVector->Energy()-fMotherCandidateLorentzVector->Pz())));
1711         return rapidity;
1712
1713 }
1714
1715
1716
1717
1718
1719 Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){
1720         //see header file for documentation
1721
1722         Int_t iResult = 10; // Unknown particle
1723
1724         if(chargeOfTrack==-1){ //negative track
1725                 switch(abs(fNegativeTrackPID)){
1726                 case 11:                         //electron
1727                         iResult = 0;
1728                         break;
1729                 case 13:                         //muon
1730                         iResult = 1;
1731                         break;
1732                 case 211:                       //pion
1733                         iResult = 2;
1734                         break;
1735                 case 321:                       //kaon
1736                         iResult = 3;
1737                         break;
1738                 case 2212:               //proton
1739                         iResult = 4;
1740                         break;
1741                 case 22:                         //photon
1742                         iResult = 5;
1743                         break;
1744                 case 111:                       //pi0
1745                         iResult = 6;
1746                         break;
1747                 case 2112:               //neutron
1748                         iResult = 7;
1749                         break;
1750                 case 311:                       //K0
1751                         iResult = 8;
1752                         break;
1753
1754                         //Put in here for kSPECIES::kEleCon     ????
1755                 }
1756         }
1757         else if(chargeOfTrack==1){ //positive track
1758                 switch(abs(fPositiveTrackPID)){
1759                 case 11:                         //electron
1760                         iResult = 0;
1761                         break;
1762                 case 13:                         //muon
1763                         iResult = 1;
1764                         break;
1765                 case 211:                       //pion
1766                         iResult = 2;
1767                         break;
1768                 case 321:                       //kaon
1769                         iResult = 3;
1770                         break;
1771                 case 2212:               //proton
1772                         iResult = 4;
1773                         break;
1774                 case 22:                         //photon
1775                         iResult = 5;
1776                         break;
1777                 case 111:                       //pi0
1778                         iResult = 6;
1779                         break;
1780                 case 2112:               //neutron
1781                         iResult = 7;
1782                         break;
1783                 case 311:                       //K0
1784                         iResult = 8;
1785                         break;
1786
1787                         //Put in here for kSPECIES::kEleCon     ????
1788                 }
1789         }
1790         else{
1791                 //Wrong parameter.. Print warning
1792         }
1793         return iResult;
1794 }
1795
1796 //_____________________________________________________________________________________________________________________
1797 Bool_t AliV0Reader::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
1798         // see header file for documentation
1799
1800         Double_t        helix[6];
1801         track->GetHelixParameters(helix,b);
1802
1803         Double_t xpos = helix[5];
1804         Double_t ypos = helix[0];
1805         Double_t radius = TMath::Abs(1./helix[4]);
1806         Double_t phi = helix[2];
1807
1808         if(phi < 0){
1809                 phi = phi + 2*TMath::Pi();
1810         }
1811
1812         phi -= TMath::Pi()/2.;
1813         Double_t xpoint =       radius * TMath::Cos(phi);
1814         Double_t ypoint =       radius * TMath::Sin(phi);
1815
1816         if(b<0){
1817                 if(charge > 0){
1818                         xpoint = - xpoint;
1819                         ypoint = - ypoint;
1820                 }
1821
1822                 if(charge < 0){
1823                         xpoint =        xpoint;
1824                         ypoint =        ypoint;
1825                 }
1826         }
1827         if(b>0){
1828                 if(charge > 0){
1829                         xpoint =        xpoint;
1830                         ypoint =        ypoint;
1831                 }
1832
1833                 if(charge < 0){
1834                         xpoint = - xpoint;
1835                         ypoint = - ypoint;
1836                 }
1837         }
1838         center[0] =     xpos + xpoint;
1839         center[1] =     ypos + ypoint;
1840
1841         return 1;
1842 }
1843
1844 //_________________________________________________________________________________________________________
1845 // Bool_t AliV0Reader::GetHelixCenter(AliESDtrack* track, Double_t b,Int_t charge, Double_t center[2]){
1846 // // Bool_t AliV0Reader::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
1847 //      // see header file for documentation
1848 //
1849 //      Double_t pi = 3.14159265358979323846;
1850 //
1851 //      Double_t        helix[6];
1852 //      track->GetHelixParameters(helix,b);
1853 //
1854 //      Double_t xpos = helix[5];
1855 //      Double_t ypos = helix[0];
1856 //      Double_t radius = TMath::Abs(1./helix[4]);
1857 //      Double_t phi = helix[2];
1858 //
1859 //      if(phi < 0){
1860 //              phi = phi + 2*pi;
1861 //      }
1862 //
1863 //      phi -= pi/2.;
1864 //      Double_t xpoint =       radius * TMath::Cos(phi);
1865 //      Double_t ypoint =       radius * TMath::Sin(phi);
1866 //
1867 //      if(b<0){
1868 //              if(charge > 0){
1869 //                      xpoint = - xpoint;
1870 //                      ypoint = - ypoint;
1871 //              }
1872 //
1873 //              if(charge < 0){
1874 //                      xpoint =        xpoint;
1875 //                      ypoint =        ypoint;
1876 //              }
1877 //      }
1878 //      if(b>0){
1879 //              if(charge > 0){
1880 //                      xpoint =        xpoint;
1881 //                      ypoint =        ypoint;
1882 //              }
1883 //
1884 //              if(charge < 0){
1885 //                      xpoint = - xpoint;
1886 //                      ypoint = - ypoint;
1887 //              }
1888 //      }
1889 //      center[0] =     xpos + xpoint;
1890 //      center[1] =     ypos + ypoint;
1891 //
1892 //      return 1;
1893 // }
1894
1895 //____________________________________________________________________________________________________________________________
1896 Bool_t AliV0Reader::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3]){
1897
1898     if(!pparam||!nparam)return kFALSE;
1899
1900         Double_t helixcenterpos[2];
1901         GetHelixCenter(pparam,GetMagneticField(),pparam->Charge(),helixcenterpos);
1902
1903         Double_t helixcenterneg[2];
1904         GetHelixCenter(nparam,GetMagneticField(),nparam->Charge(),helixcenterneg);
1905
1906         Double_t helixpos[6];
1907         pparam->GetHelixParameters(helixpos,GetMagneticField());
1908         Double_t posradius = TMath::Abs(1./helixpos[4]);
1909
1910         Double_t helixneg[6];
1911         nparam->GetHelixParameters(helixneg,GetMagneticField());
1912         Double_t negradius = TMath::Abs(1./helixneg[4]);
1913
1914         // Calculate xy-position
1915
1916         Double_t xpos = helixcenterpos[0];
1917         Double_t ypos = helixcenterpos[1];
1918         Double_t xneg = helixcenterneg[0];
1919         Double_t yneg = helixcenterneg[1];
1920
1921         convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1922         convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
1923
1924
1925         // Calculate z-position
1926
1927         Double_t deltaXPos = convpos[0] -       xpos;
1928         Double_t deltaYPos = convpos[1] -       ypos;
1929
1930         Double_t deltaXNeg = convpos[0] -       xneg;
1931         Double_t deltaYNeg = convpos[1] -       yneg;
1932
1933         Double_t alphaPos =     TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
1934         Double_t alphaNeg =     TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
1935
1936         Double_t vertexXNeg =   xneg +  TMath::Abs(negradius)*TMath::Cos(alphaNeg);
1937         Double_t vertexYNeg =   yneg +  TMath::Abs(negradius)*TMath::Sin(alphaNeg);
1938
1939         Double_t vertexXPos =   xpos +  TMath::Abs(posradius)*TMath::Cos(alphaPos);
1940         Double_t vertexYPos =   ypos +  TMath::Abs(posradius)*TMath::Sin(alphaPos);
1941
1942         Double_t b = fESDEvent->GetMagneticField();
1943
1944         AliExternalTrackParam p(*pparam);
1945         AliExternalTrackParam n(*nparam);
1946
1947         TVector2 vertexPos(vertexXPos,vertexYPos);
1948         TVector2 vertexNeg(vertexXNeg,vertexYNeg);
1949
1950         // Convert to local coordinate system
1951         vertexPos=vertexPos.Rotate(-p.GetAlpha());
1952         vertexNeg=vertexNeg.Rotate(-p.GetAlpha());
1953
1954         // Propagate Track Params to Vertex
1955         p.PropagateTo(vertexPos.X(),b);
1956         n.PropagateTo(vertexNeg.X(),b);
1957
1958         convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
1959
1960         return kTRUE;
1961 }
1962 //
1963 // //__________________________________________________________________________________________________________
1964 Bool_t AliV0Reader::GetConvPosXY(AliESDtrack* ptrack, AliESDtrack* ntrack, Double_t b, Double_t convpos[2]){
1965         //see header file for documentation
1966
1967         Double_t helixcenterpos[2];
1968         GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
1969
1970         Double_t helixcenterneg[2];
1971         GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
1972
1973         Double_t        poshelix[6];
1974         ptrack->GetHelixParameters(poshelix,b);
1975         Double_t posradius = TMath::Abs(1./poshelix[4]);
1976
1977         Double_t        neghelix[6];
1978         ntrack->GetHelixParameters(neghelix,b);
1979         Double_t negradius = TMath::Abs(1./neghelix[4]);
1980
1981         Double_t xpos = helixcenterpos[0];
1982         Double_t ypos = helixcenterpos[1];
1983         Double_t xneg = helixcenterneg[0];
1984         Double_t yneg = helixcenterneg[1];
1985
1986         convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1987         convpos[1] = (ypos*negradius+   yneg*posradius)/(negradius+posradius);
1988
1989         return 1;
1990 }
1991 //
1992 //
1993 //
1994 Double_t AliV0Reader::GetConvPosZ(AliESDtrack* ptrack,AliESDtrack* ntrack, Double_t b){
1995         //see header file for documentation
1996
1997         Double_t        helixpos[6];
1998         ptrack->GetHelixParameters(helixpos,b);
1999
2000         Double_t        helixneg[6];
2001         ntrack->GetHelixParameters(helixneg,b);
2002
2003         Double_t negtrackradius =       TMath::Abs(1./helixneg[4]);
2004         Double_t postrackradius =       TMath::Abs(1./helixpos[4]);
2005
2006         Double_t pi = 3.14159265358979323846;
2007
2008         Double_t convpos[2];
2009         GetConvPosXY(ptrack,ntrack,b,convpos);
2010
2011          Double_t convposx = convpos[0];
2012          Double_t convposy = convpos[1];
2013
2014          Double_t helixcenterpos[2];
2015          GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
2016
2017          Double_t helixcenterneg[2];
2018          GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
2019
2020          Double_t xpos = helixcenterpos[0];
2021          Double_t ypos = helixcenterpos[1];
2022          Double_t xneg = helixcenterneg[0];
2023          Double_t yneg = helixcenterneg[1];
2024
2025          Double_t deltaXPos = convposx -        xpos;
2026          Double_t deltaYPos = convposy -        ypos;
2027
2028          Double_t deltaXNeg = convposx -        xneg;
2029          Double_t deltaYNeg = convposy -        yneg;
2030
2031          Double_t alphaPos =    pi + TMath::ATan2(-deltaYPos,-deltaXPos);
2032          Double_t alphaNeg =    pi + TMath::ATan2(-deltaYNeg,-deltaXNeg);
2033
2034          Double_t vertexXNeg =  xneg +  TMath::Abs(negtrackradius)*
2035          TMath::Cos(alphaNeg);
2036          Double_t vertexYNeg =  yneg +  TMath::Abs(negtrackradius)*
2037          TMath::Sin(alphaNeg);
2038
2039          Double_t vertexXPos =  xpos +  TMath::Abs(postrackradius)*
2040          TMath::Cos(alphaPos);
2041          Double_t vertexYPos =  ypos +  TMath::Abs(postrackradius)*
2042          TMath::Sin(alphaPos);
2043
2044          Double_t x0neg =        helixneg[5];
2045          Double_t y0neg =        helixneg[0];
2046
2047          Double_t x0pos =        helixpos[5];
2048          Double_t y0pos =        helixpos[0];
2049
2050          Double_t dNeg = TMath::Sqrt((vertexXNeg -      x0neg)*(vertexXNeg - x0neg)
2051                                                                                                                          +(vertexYNeg - y0neg)*(vertexYNeg - y0neg));
2052
2053          Double_t dPos = TMath::Sqrt((vertexXPos -      x0pos)*(vertexXPos - x0pos)
2054                                                                                                                          +(vertexYPos - y0pos)*(vertexYPos - y0pos));
2055
2056          Double_t rNeg =        TMath::Sqrt(negtrackradius*negtrackradius -
2057          dNeg*dNeg/4.);
2058
2059          Double_t rPos = TMath::Sqrt(postrackradius*postrackradius -
2060          dPos*dPos/4.);
2061
2062          Double_t deltabetaNeg =        2*(pi +  TMath::ATan2(-dNeg/2.,-rNeg));
2063          Double_t deltabetaPos = 2*(pi + TMath::ATan2(-dPos/2.,-rPos));
2064
2065          Double_t deltaUNeg = negtrackradius*deltabetaNeg;
2066          Double_t deltaUPos = postrackradius*deltabetaPos;
2067
2068          Double_t zphaseNeg = ntrack->GetZ() +  deltaUNeg * ntrack->GetTgl();
2069          Double_t zphasePos = ptrack->GetZ() +  deltaUPos * ptrack->GetTgl();
2070
2071          Double_t convposz = (zphasePos*negtrackradius+zphaseNeg*postrackradius)/(negtrackradius+postrackradius);
2072
2073          return convposz;
2074 }
2075 //
2076 AliGammaConversionKFVector* AliV0Reader::GetBGGoodV0s(Int_t /*event*/) const{
2077         /*
2078         if(fUseChargedTrackMultiplicityForBG == kTRUE){
2079                 return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
2080         }
2081         else{ // means we use #v0s as multiplicity
2082                 return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
2083         }
2084         */
2085         return NULL;
2086 }
2087
2088 Int_t AliV0Reader::CountESDTracks(){
2089         // see header file for documentation
2090         if(fNumberOfESDTracks == 0){ // count the good esd tracks
2091                 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2092                         AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
2093                         if(!curTrack){
2094         continue;
2095                         }
2096                         if(fEsdTrackCuts->AcceptTrack(curTrack) ){
2097         fNumberOfESDTracks++;
2098                         }
2099                 }
2100         }
2101
2102         return fNumberOfESDTracks;
2103 }
2104
2105 Bool_t AliV0Reader::CheckIfPi0IsMother(Int_t label){
2106         // see headerfile for documentation
2107         Bool_t iResult=kFALSE;
2108         //      cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
2109         if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 111){
2110                 iResult=kTRUE;
2111         }
2112         return iResult;
2113 }
2114
2115 Bool_t AliV0Reader::CheckIfEtaIsMother(Int_t label){
2116         // see headerfile for documentation
2117         Bool_t iResult=kFALSE;
2118         //      cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
2119         if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 221){
2120                 iResult=kTRUE;
2121         }
2122         return iResult;
2123 }
2124
2125
2126
2127
2128 Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliKFParticle* negativeKFParticle, const AliKFParticle * positiveKFParticle, Double_t armenterosQtAlpha[2], Int_t kfProductionMethod){
2129
2130    AliKFParticle PosParticle = *positiveKFParticle;
2131    AliKFParticle NegParticle = *negativeKFParticle;
2132    AliKFParticle Gamma;
2133    if(kfProductionMethod < 3)
2134       Gamma.ConstructGamma(PosParticle, NegParticle);
2135    else if(kfProductionMethod == 3){
2136       Gamma += PosParticle;
2137       Gamma += NegParticle;
2138    }
2139
2140    Double_t VertexGamma[3] = {Gamma.GetX(), Gamma.GetY(), Gamma.GetZ()};
2141    PosParticle.TransportToPoint(VertexGamma);
2142    NegParticle.TransportToPoint(VertexGamma);
2143
2144    AliKFParticle::GetArmenterosPodolanski(PosParticle, NegParticle, armenterosQtAlpha);
2145
2146    return 1;
2147 }
2148
2149
2150
2151
2152 Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliKFParticle* negativeKFParticle, const AliKFParticle * positiveKFParticle, const AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlpha[2] ){
2153         //see header file for documentation
2154
2155         TVector3 momentumVectorPositiveKF(positiveKFParticle->GetPx(),positiveKFParticle->GetPy(),positiveKFParticle->GetPz());
2156         TVector3 momentumVectorNegativeKF(negativeKFParticle->GetPx(),negativeKFParticle->GetPy(),negativeKFParticle->GetPz());
2157         TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
2158
2159         Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2160         Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2161
2162         Float_t alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2163                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2164
2165
2166         Float_t qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2167
2168         armenterosQtAlpha[0]=qt;
2169         armenterosQtAlpha[1]=alfa;
2170
2171         return 1;
2172
2173 }
2174
2175 Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliESDv0* v0, Double_t armenterosQtAlpha[2])
2176 {       //see header file for documentation
2177
2178    Double_t mn[3] = {0,0,0};
2179    Double_t mp[3] = {0,0,0};
2180    Double_t mm[3] = {0,0,0};
2181
2182
2183    Int_t pIndex = 0, nIndex = 0;
2184    pIndex = v0->GetPindex();
2185    nIndex = v0->GetNindex();
2186
2187    AliESDtrack* d[2];
2188    d[0] = dynamic_cast<AliESDtrack*>(fESDEvent->GetTrack(pIndex));
2189    d[1] = dynamic_cast<AliESDtrack*>(fESDEvent->GetTrack(nIndex));
2190
2191    Int_t sign[2];
2192    sign[0] = (int)d[0]->GetSign();
2193    sign[1] = (int)d[1]->GetSign();
2194
2195    Bool_t correct = kFALSE;
2196
2197
2198    if(-1 == sign[0] && 1 == sign[1]){
2199       correct = kFALSE;
2200    }
2201    else{
2202       correct = kTRUE;
2203    }
2204
2205
2206    if(correct){
2207       v0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
2208       v0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
2209    }
2210    else{
2211       v0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
2212       v0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
2213    }
2214    v0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
2215
2216    TVector3 vecN(mn[0],mn[1],mn[2]);
2217    TVector3 vecP(mp[0],mp[1],mp[2]);
2218    TVector3 vecM(mm[0],mm[1],mm[2]);
2219
2220    Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
2221    Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
2222
2223    Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
2224       ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
2225    Double_t qt = vecP.Mag()*sin(thetaP);
2226
2227    armenterosQtAlpha[0]=qt;
2228    armenterosQtAlpha[1]=alfa;
2229
2230    return 1;
2231
2232 }
2233
2234
2235
2236 Bool_t AliV0Reader::GetArmenterosQtAlpha(const TParticle* negativeParticle, const TParticle * positiveParticle, const AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlpha[2] ){
2237         //see header file for documentation
2238
2239         TVector3 momentumVectorPositiveKF(positiveParticle->Px(),positiveParticle->Py(),positiveParticle->Pz());
2240         TVector3 momentumVectorNegativeKF(negativeParticle->Px(),negativeParticle->Py(),negativeParticle->Pz());
2241         TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
2242
2243         Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2244         Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2245
2246         Float_t alfa;
2247         Float_t qt;
2248         if ( positiveParticle->GetPdgCode() == 11 || positiveParticle->GetPdgCode() == 13 || positiveParticle->GetPdgCode() == 15){
2249                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2250                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2251                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2252         } else if ( negativeParticle->GetPdgCode() == -11 || negativeParticle->GetPdgCode() == -13 || negativeParticle->GetPdgCode() == -15){
2253                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2254                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2255                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2256         } else if (positiveParticle->GetPdgCode() < 0 && positiveParticle->GetPdgCode() != -11 && positiveParticle->GetPdgCode() != -13 && positiveParticle->GetPdgCode() != -15 ){
2257                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2258                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2259                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2260         } else {
2261                 alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2262                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2263                 qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2264         }
2265
2266         armenterosQtAlpha[0]=qt;
2267         armenterosQtAlpha[1]=alfa;
2268         return 1;
2269
2270 }
2271
2272 Bool_t AliV0Reader::GetArmenterosQtAlpha(const TParticle* negativeParticle, const TParticle * positiveParticle, const TParticle * gammaCandidate, Double_t armenterosQtAlpha[2] ){
2273         //see header file for documentation
2274
2275         TVector3 momentumVectorPositiveKF(positiveParticle->Px(),positiveParticle->Py(),positiveParticle->Pz());
2276         TVector3 momentumVectorNegativeKF(negativeParticle->Px(),negativeParticle->Py(),negativeParticle->Pz());
2277         TVector3 vecV0(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
2278
2279         Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2280         Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2281
2282         Float_t alfa;
2283         Float_t qt;
2284         if ( positiveParticle->GetPdgCode() == 11 || positiveParticle->GetPdgCode() == 13 || positiveParticle->GetPdgCode() == 15){
2285                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2286                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2287                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2288         } else if ( negativeParticle->GetPdgCode() == -11 || negativeParticle->GetPdgCode() == -13 || negativeParticle->GetPdgCode() == -15){
2289                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2290                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2291                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2292         } else if (positiveParticle->GetPdgCode() < 0 && positiveParticle->GetPdgCode() != -11 && positiveParticle->GetPdgCode() != -13 && positiveParticle->GetPdgCode() != -15 ){
2293                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2294                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2295                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2296         } else {
2297                 alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2298                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2299                 qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2300         }
2301
2302         armenterosQtAlpha[0]=qt;
2303         armenterosQtAlpha[1]=alfa;
2304         return 1;
2305
2306 }
2307
2308
2309 Int_t AliV0Reader::GetFirstTPCRow(Double_t radius){
2310
2311
2312         Int_t firstTPCRow=0;
2313         Double_t radiusI        =       84.8;
2314         Double_t radiusO        = 134.6;
2315         Double_t radiusOB = 198.;
2316         Double_t rSizeI  = 0.75;
2317         Double_t rSizeO  = 1.;
2318         Double_t rSizeOB        = 1.5;
2319         Int_t nClsI=63;
2320         Int_t nClsIO=127;
2321
2322         if(radius <= radiusI){
2323                 return firstTPCRow;
2324         }
2325         if(radius>radiusI && radius<=radiusO){
2326                 firstTPCRow = (Int_t)((radius-radiusI)/rSizeI);
2327         }
2328         if(radius>radiusO && radius<=radiusOB){
2329                 firstTPCRow = (Int_t)(nClsI+(radius-radiusO)/rSizeO);
2330         }
2331
2332         if(radius>radiusOB){
2333                 firstTPCRow =(Int_t)(nClsIO+(radius-radiusOB)/rSizeOB);
2334         }
2335
2336
2337         return firstTPCRow;
2338 }
2339 void AliV0Reader::SmearKFParticle(AliKFParticle * kfParticle)
2340 {
2341         Double_t facPBrem = 1.;
2342         Double_t facPSig = 0.;
2343
2344         Double_t phi=0.;
2345         Double_t theta=0.;
2346         Double_t P=0.;
2347
2348         P=kfParticle->GetP();
2349         phi=kfParticle->GetPhi();
2350         if( kfParticle->GetP()!=0){
2351                 theta=acos( kfParticle->Pz()/ kfParticle->GetP());
2352         }
2353
2354         if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
2355                 facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
2356         }
2357
2358         if( fPBremSmearing != 1.){
2359                 if(fBrem!=NULL){
2360                         facPBrem = fBrem->GetRandom();
2361                 }
2362         }
2363
2364         kfParticle->Px() = facPBrem* (1+facPSig)* P*sin(theta)*cos(phi) ;
2365         kfParticle->Py() = facPBrem* (1+facPSig)* P*sin(theta)*sin(phi) ;
2366         kfParticle->Pz() = facPBrem* (1+facPSig)* P*cos(theta) ;
2367         kfParticle->E() = kfParticle->GetP();
2368 }
2369
2370 ///________________________________________________________________________
2371 const AliExternalTrackParam *AliV0Reader::GetExternalTrackParam(AliESDv0 *v0,Int_t charge){
2372
2373     if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
2374
2375     Int_t label;
2376     if(charge>0)label=0;
2377     else label=1;
2378     // Check for sign flip
2379
2380     if(v0){
2381         if(!v0->GetParamN()||!v0->GetParamP())return 0x0;
2382         if(!GetTrack(v0->GetNindex())||!fESDEvent->GetTrack(v0->GetPindex()))return 0x0;
2383         if((GetTrack(v0->GetPindex()))->Charge()==charge){
2384 //          fCurrentTrackLabels[label]=v0->GetPindex();
2385             return v0->GetParamP();}
2386         if((GetTrack(v0->GetNindex()))->Charge()==charge){
2387 //          fCurrentTrackLabels[label]=v0->GetNindex();
2388             return v0->GetParamN();}
2389     }
2390     return 0x0;
2391 }
2392
2393 ///________________________________________________________________________
2394 AliVTrack *AliV0Reader::GetTrack(Int_t label){
2395     if(fESDEvent){
2396                 return (AliESDtrack*)fESDEvent->GetTrack(label);
2397     }
2398 //     if(fAODEvent)return (AliAODTrack*)GetAODTrack(label);
2399         return 0x0;
2400 }
2401
2402 Double_t AliV0Reader::GetV0CosineOfPointingAngle(Double_t V0PointX, Double_t V0PointY, Double_t V0PointZ){
2403    // calculates the pointing angle of the recalculated V0
2404
2405    Double_t momV0[3]; //momentum of the V0
2406    fCurrentV0->GetPxPyPz(momV0[0],momV0[1],momV0[2]);
2407
2408    Double_t PosV0[3]; //Recalculated V0 Position vector
2409
2410    PosV0[0] = V0PointX - fESDEvent->GetPrimaryVertex()->GetX();
2411    PosV0[1] = V0PointY - fESDEvent->GetPrimaryVertex()->GetY();
2412    PosV0[2] = V0PointZ - fESDEvent->GetPrimaryVertex()->GetZ();
2413
2414    Double_t momV02 = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
2415    Double_t PosV02 = PosV0[0]*PosV0[0] + PosV0[1]*PosV0[1] + PosV0[2]*PosV0[2];
2416
2417    Double_t cosinePointingAngle = (PosV0[0]*momV0[0] +  PosV0[1]*momV0[1] + PosV0[2]*momV0[2] ) / TMath::Sqrt(momV02 * PosV02);
2418
2419    return cosinePointingAngle;
2420 }
2421
2422
2423 Double_t AliV0Reader::GetPsiPair(AliESDv0* v0)
2424 {
2425    //
2426    // Angle between daughter momentum plane and plane
2427    //
2428    Float_t magField = fESDEvent->GetMagneticField();
2429
2430    Double_t xyz[3] = {0.,0.,0.};
2431    v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
2432
2433    Double_t mn[3] = {0,0,0};
2434    Double_t mp[3] = {0,0,0};
2435
2436    v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
2437    v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
2438
2439    Double_t deltat = 1.;
2440    deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) -  TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
2441    Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
2442
2443    Double_t momPosProp[3] = {0,0,0};
2444    Double_t momNegProp[3] = {0,0,0};
2445
2446    AliExternalTrackParam nt = *(v0->GetParamN());
2447    AliExternalTrackParam pt = *(v0->GetParamP());
2448
2449    Double_t psiPair = 4.;
2450    if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
2451
2452    if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
2453
2454    pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
2455    nt.GetPxPyPz(momNegProp);
2456
2457   Double_t pEle =
2458      TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
2459
2460   Double_t pPos =
2461      TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
2462
2463   Double_t scalarproduct =
2464      momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
2465
2466   Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
2467
2468   psiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
2469
2470   return psiPair;
2471 }
2472
2473 Int_t AliV0Reader::GetNumberOfHijingPlusPythiaPrimeries(Int_t excludeHeaderType){
2474
2475    // Calculate NPrimaries for LHC11a10b_*
2476
2477    Int_t nproduced = 0;
2478    AliGenCocktailEventHeader *cHeader = dynamic_cast<AliGenCocktailEventHeader*>(fMCEvent->GenEventHeader());
2479    if(cHeader){
2480       TList *genHeaders = cHeader->GetHeaders();
2481       AliGenEventHeader* gh = 0;
2482       for(Int_t i = 0; i<genHeaders->GetEntries();i++){
2483          gh = (AliGenEventHeader*)genHeaders->At(i);
2484          TString GeneratorName = gh->GetName();
2485
2486          if(GeneratorName.CompareTo("Hijing") == 0){
2487             nproduced = nproduced + gh->NProduced();
2488             //            cout<<i<<"   "<<GeneratorName<<"   "<<gh->NProduced()<<endl;
2489          }
2490          else if(GeneratorName.CompareTo("Pythia") == 0 && excludeHeaderType == 1){
2491             nproduced = nproduced + gh->NProduced();
2492             //            cout<<i<<"   "<<GeneratorName<<"   "<<gh->NProduced()<<endl;
2493          }
2494       }
2495    }
2496    if(!cHeader){
2497       nproduced = fMCStack->GetNprimary();
2498    }
2499
2500    //  cout<<fMCStack->GetNprimary()-nproduced<<endl;
2501
2502    return nproduced;
2503 }
2504
2505
2506 Bool_t AliV0Reader::IsParticleFromBGEvent(Int_t index){
2507
2508    Bool_t particleFromBG = kFALSE;
2509
2510    if(index == -1) return kFALSE;
2511    if(index > fNumberOfPrimerisFromHijingAndPythia && index < fMCStack->GetNprimary()){
2512       //      cout<<fMCEvent->IsFromBGEvent(index)<<endl;
2513       // cout<<index<<"   "<<fNumberOfPrimerisFromHijingAndPythia<<endl;
2514       return kTRUE;
2515    }
2516    //  else cout<<"Passt Noch "<<index<<"   "<<fNumberOfPrimerisFromHijingAndPythia<<endl;
2517    // cout<<fMCEvent->IsFromBGEvent(index)<<endl;
2518    TParticle *BGParticle = fMCStack->Particle(index);
2519    if(BGParticle->GetMother(0) > -1) return kFALSE;
2520    Int_t indexMother = fMCStack->Particle(index)->GetMother(0);
2521    particleFromBG = IsParticleFromBGEvent(indexMother);
2522
2523    return kFALSE;
2524 }