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[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[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();
 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
1654
1655         //      fBGEventHandler->PrintBGArray(); // for debugging
1656 }
1657
1658
1659 Double_t AliV0Reader::GetNegativeTrackPhi() const{
1660         //see header file for documentation
1661
1662         Double_t offset=0;
1663         if(fNegativeTrackLorentzVector->Phi()> TMath::Pi()){
1664                 offset = -2*TMath::Pi();
1665         }
1666         return fNegativeTrackLorentzVector->Phi()+offset;
1667 }
1668
1669 Double_t AliV0Reader::GetPositiveTrackPhi() const{
1670         //see header file for documentation
1671
1672         Double_t offset=0;
1673         if(fPositiveTrackLorentzVector->Phi()> TMath::Pi()){
1674                 offset = -2*TMath::Pi();
1675         }
1676         return fPositiveTrackLorentzVector->Phi()+offset;
1677 }
1678
1679 Double_t AliV0Reader::GetMotherCandidatePhi() const{
1680         //see header file for documentation
1681
1682         Double_t offset=0;
1683         if(fMotherCandidateLorentzVector->Phi()> TMath::Pi()){
1684                 offset = -2*TMath::Pi();
1685         }
1686         return fMotherCandidateLorentzVector->Phi()+offset;
1687 }
1688
1689
1690 Double_t AliV0Reader::GetMotherCandidateRapidity() const{
1691         //see header file for documentation
1692
1693         Double_t rapidity=0;
1694         if(fMotherCandidateLorentzVector->Energy() - fMotherCandidateLorentzVector->Pz() == 0 || fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz() == 0) rapidity=0;
1695         else rapidity = 0.5*(TMath::Log((fMotherCandidateLorentzVector->Energy() + fMotherCandidateLorentzVector->Pz()) / (fMotherCandidateLorentzVector->Energy()-fMotherCandidateLorentzVector->Pz())));
1696         return rapidity;
1697
1698 }
1699
1700
1701
1702
1703
1704 Int_t AliV0Reader::GetSpeciesIndex(Int_t chargeOfTrack){
1705         //see header file for documentation
1706
1707         Int_t iResult = 10; // Unknown particle
1708
1709         if(chargeOfTrack==-1){ //negative track
1710                 switch(abs(fNegativeTrackPID)){
1711                 case 11:                         //electron
1712                         iResult = 0;
1713                         break;
1714                 case 13:                         //muon
1715                         iResult = 1;
1716                         break;
1717                 case 211:                       //pion
1718                         iResult = 2;
1719                         break;
1720                 case 321:                       //kaon
1721                         iResult = 3;
1722                         break;
1723                 case 2212:               //proton
1724                         iResult = 4;
1725                         break;
1726                 case 22:                         //photon
1727                         iResult = 5;
1728                         break;
1729                 case 111:                       //pi0
1730                         iResult = 6;
1731                         break;
1732                 case 2112:               //neutron
1733                         iResult = 7;
1734                         break;
1735                 case 311:                       //K0
1736                         iResult = 8;
1737                         break;
1738
1739                         //Put in here for kSPECIES::kEleCon     ????
1740                 }
1741         }
1742         else if(chargeOfTrack==1){ //positive track
1743                 switch(abs(fPositiveTrackPID)){
1744                 case 11:                         //electron
1745                         iResult = 0;
1746                         break;
1747                 case 13:                         //muon
1748                         iResult = 1;
1749                         break;
1750                 case 211:                       //pion
1751                         iResult = 2;
1752                         break;
1753                 case 321:                       //kaon
1754                         iResult = 3;
1755                         break;
1756                 case 2212:               //proton
1757                         iResult = 4;
1758                         break;
1759                 case 22:                         //photon
1760                         iResult = 5;
1761                         break;
1762                 case 111:                       //pi0
1763                         iResult = 6;
1764                         break;
1765                 case 2112:               //neutron
1766                         iResult = 7;
1767                         break;
1768                 case 311:                       //K0
1769                         iResult = 8;
1770                         break;
1771
1772                         //Put in here for kSPECIES::kEleCon     ????
1773                 }
1774         }
1775         else{
1776                 //Wrong parameter.. Print warning
1777         }
1778         return iResult;
1779 }
1780
1781 //_____________________________________________________________________________________________________________________
1782 Bool_t AliV0Reader::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
1783         // see header file for documentation
1784
1785         Double_t        helix[6];
1786         track->GetHelixParameters(helix,b);
1787
1788         Double_t xpos = helix[5];
1789         Double_t ypos = helix[0];
1790         Double_t radius = TMath::Abs(1./helix[4]);
1791         Double_t phi = helix[2];
1792
1793         if(phi < 0){
1794                 phi = phi + 2*TMath::Pi();
1795         }
1796
1797         phi -= TMath::Pi()/2.;
1798         Double_t xpoint =       radius * TMath::Cos(phi);
1799         Double_t ypoint =       radius * TMath::Sin(phi);
1800
1801         if(b<0){
1802                 if(charge > 0){
1803                         xpoint = - xpoint;
1804                         ypoint = - ypoint;
1805                 }
1806
1807                 if(charge < 0){
1808                         xpoint =        xpoint;
1809                         ypoint =        ypoint;
1810                 }
1811         }
1812         if(b>0){
1813                 if(charge > 0){
1814                         xpoint =        xpoint;
1815                         ypoint =        ypoint;
1816                 }
1817
1818                 if(charge < 0){
1819                         xpoint = - xpoint;
1820                         ypoint = - ypoint;
1821                 }
1822         }
1823         center[0] =     xpos + xpoint;
1824         center[1] =     ypos + ypoint;
1825
1826         return 1;
1827 }
1828
1829 //_________________________________________________________________________________________________________
1830 // Bool_t AliV0Reader::GetHelixCenter(AliESDtrack* track, Double_t b,Int_t charge, Double_t center[2]){
1831 // // Bool_t AliV0Reader::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
1832 //      // see header file for documentation
1833 //
1834 //      Double_t pi = 3.14159265358979323846;
1835 //
1836 //      Double_t        helix[6];
1837 //      track->GetHelixParameters(helix,b);
1838 //
1839 //      Double_t xpos = helix[5];
1840 //      Double_t ypos = helix[0];
1841 //      Double_t radius = TMath::Abs(1./helix[4]);
1842 //      Double_t phi = helix[2];
1843 //
1844 //      if(phi < 0){
1845 //              phi = phi + 2*pi;
1846 //      }
1847 //
1848 //      phi -= pi/2.;
1849 //      Double_t xpoint =       radius * TMath::Cos(phi);
1850 //      Double_t ypoint =       radius * TMath::Sin(phi);
1851 //
1852 //      if(b<0){
1853 //              if(charge > 0){
1854 //                      xpoint = - xpoint;
1855 //                      ypoint = - ypoint;
1856 //              }
1857 //
1858 //              if(charge < 0){
1859 //                      xpoint =        xpoint;
1860 //                      ypoint =        ypoint;
1861 //              }
1862 //      }
1863 //      if(b>0){
1864 //              if(charge > 0){
1865 //                      xpoint =        xpoint;
1866 //                      ypoint =        ypoint;
1867 //              }
1868 //
1869 //              if(charge < 0){
1870 //                      xpoint = - xpoint;
1871 //                      ypoint = - ypoint;
1872 //              }
1873 //      }
1874 //      center[0] =     xpos + xpoint;
1875 //      center[1] =     ypos + ypoint;
1876 //
1877 //      return 1;
1878 // }
1879
1880 //____________________________________________________________________________________________________________________________
1881 Bool_t AliV0Reader::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3]){
1882
1883     if(!pparam||!nparam)return kFALSE;
1884
1885         Double_t helixcenterpos[2];
1886         GetHelixCenter(pparam,GetMagneticField(),pparam->Charge(),helixcenterpos);
1887
1888         Double_t helixcenterneg[2];
1889         GetHelixCenter(nparam,GetMagneticField(),nparam->Charge(),helixcenterneg);
1890
1891         Double_t helixpos[6];
1892         pparam->GetHelixParameters(helixpos,GetMagneticField());
1893         Double_t posradius = TMath::Abs(1./helixpos[4]);
1894
1895         Double_t helixneg[6];
1896         nparam->GetHelixParameters(helixneg,GetMagneticField());
1897         Double_t negradius = TMath::Abs(1./helixneg[4]);
1898
1899         // Calculate xy-position
1900
1901         Double_t xpos = helixcenterpos[0];
1902         Double_t ypos = helixcenterpos[1];
1903         Double_t xneg = helixcenterneg[0];
1904         Double_t yneg = helixcenterneg[1];
1905
1906         convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1907         convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
1908
1909
1910         // Calculate z-position
1911
1912         Double_t deltaXPos = convpos[0] -       xpos;
1913         Double_t deltaYPos = convpos[1] -       ypos;
1914
1915         Double_t deltaXNeg = convpos[0] -       xneg;
1916         Double_t deltaYNeg = convpos[1] -       yneg;
1917
1918         Double_t alphaPos =     TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
1919         Double_t alphaNeg =     TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
1920
1921         Double_t vertexXNeg =   xneg +  TMath::Abs(negradius)*TMath::Cos(alphaNeg);
1922         Double_t vertexYNeg =   yneg +  TMath::Abs(negradius)*TMath::Sin(alphaNeg);
1923
1924         Double_t vertexXPos =   xpos +  TMath::Abs(posradius)*TMath::Cos(alphaPos);
1925         Double_t vertexYPos =   ypos +  TMath::Abs(posradius)*TMath::Sin(alphaPos);
1926
1927         Double_t b = fESDEvent->GetMagneticField();
1928
1929         AliExternalTrackParam p(*pparam);
1930         AliExternalTrackParam n(*nparam);
1931
1932         TVector2 vertexPos(vertexXPos,vertexYPos);
1933         TVector2 vertexNeg(vertexXNeg,vertexYNeg);
1934
1935         // Convert to local coordinate system
1936         vertexPos=vertexPos.Rotate(-p.GetAlpha());
1937         vertexNeg=vertexNeg.Rotate(-p.GetAlpha());
1938
1939         // Propagate Track Params to Vertex
1940         p.PropagateTo(vertexPos.X(),b);
1941         n.PropagateTo(vertexNeg.X(),b);
1942
1943         convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
1944
1945         return kTRUE;
1946 }
1947 //
1948 // //__________________________________________________________________________________________________________
1949 Bool_t AliV0Reader::GetConvPosXY(AliESDtrack* ptrack, AliESDtrack* ntrack, Double_t b, Double_t convpos[2]){
1950         //see header file for documentation
1951
1952         Double_t helixcenterpos[2];
1953         GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
1954
1955         Double_t helixcenterneg[2];
1956         GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
1957
1958         Double_t        poshelix[6];
1959         ptrack->GetHelixParameters(poshelix,b);
1960         Double_t posradius = TMath::Abs(1./poshelix[4]);
1961
1962         Double_t        neghelix[6];
1963         ntrack->GetHelixParameters(neghelix,b);
1964         Double_t negradius = TMath::Abs(1./neghelix[4]);
1965
1966         Double_t xpos = helixcenterpos[0];
1967         Double_t ypos = helixcenterpos[1];
1968         Double_t xneg = helixcenterneg[0];
1969         Double_t yneg = helixcenterneg[1];
1970
1971         convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1972         convpos[1] = (ypos*negradius+   yneg*posradius)/(negradius+posradius);
1973
1974         return 1;
1975 }
1976 //
1977 //
1978 //
1979 Double_t AliV0Reader::GetConvPosZ(AliESDtrack* ptrack,AliESDtrack* ntrack, Double_t b){
1980         //see header file for documentation
1981
1982         Double_t        helixpos[6];
1983         ptrack->GetHelixParameters(helixpos,b);
1984
1985         Double_t        helixneg[6];
1986         ntrack->GetHelixParameters(helixneg,b);
1987
1988         Double_t negtrackradius =       TMath::Abs(1./helixneg[4]);
1989         Double_t postrackradius =       TMath::Abs(1./helixpos[4]);
1990
1991         Double_t pi = 3.14159265358979323846;
1992
1993         Double_t convpos[2];
1994         GetConvPosXY(ptrack,ntrack,b,convpos);
1995
1996          Double_t convposx = convpos[0];
1997          Double_t convposy = convpos[1];
1998
1999          Double_t helixcenterpos[2];
2000          GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
2001
2002          Double_t helixcenterneg[2];
2003          GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
2004
2005          Double_t xpos = helixcenterpos[0];
2006          Double_t ypos = helixcenterpos[1];
2007          Double_t xneg = helixcenterneg[0];
2008          Double_t yneg = helixcenterneg[1];
2009
2010          Double_t deltaXPos = convposx -        xpos;
2011          Double_t deltaYPos = convposy -        ypos;
2012
2013          Double_t deltaXNeg = convposx -        xneg;
2014          Double_t deltaYNeg = convposy -        yneg;
2015
2016          Double_t alphaPos =    pi + TMath::ATan2(-deltaYPos,-deltaXPos);
2017          Double_t alphaNeg =    pi + TMath::ATan2(-deltaYNeg,-deltaXNeg);
2018
2019          Double_t vertexXNeg =  xneg +  TMath::Abs(negtrackradius)*
2020          TMath::Cos(alphaNeg);
2021          Double_t vertexYNeg =  yneg +  TMath::Abs(negtrackradius)*
2022          TMath::Sin(alphaNeg);
2023
2024          Double_t vertexXPos =  xpos +  TMath::Abs(postrackradius)*
2025          TMath::Cos(alphaPos);
2026          Double_t vertexYPos =  ypos +  TMath::Abs(postrackradius)*
2027          TMath::Sin(alphaPos);
2028
2029          Double_t x0neg =        helixneg[5];
2030          Double_t y0neg =        helixneg[0];
2031
2032          Double_t x0pos =        helixpos[5];
2033          Double_t y0pos =        helixpos[0];
2034
2035          Double_t dNeg = TMath::Sqrt((vertexXNeg -      x0neg)*(vertexXNeg - x0neg)
2036                                                                                                                          +(vertexYNeg - y0neg)*(vertexYNeg - y0neg));
2037
2038          Double_t dPos = TMath::Sqrt((vertexXPos -      x0pos)*(vertexXPos - x0pos)
2039                                                                                                                          +(vertexYPos - y0pos)*(vertexYPos - y0pos));
2040
2041          Double_t rNeg =        TMath::Sqrt(negtrackradius*negtrackradius -
2042          dNeg*dNeg/4.);
2043
2044          Double_t rPos = TMath::Sqrt(postrackradius*postrackradius -
2045          dPos*dPos/4.);
2046
2047          Double_t deltabetaNeg =        2*(pi +  TMath::ATan2(-dNeg/2.,-rNeg));
2048          Double_t deltabetaPos = 2*(pi + TMath::ATan2(-dPos/2.,-rPos));
2049
2050          Double_t deltaUNeg = negtrackradius*deltabetaNeg;
2051          Double_t deltaUPos = postrackradius*deltabetaPos;
2052
2053          Double_t zphaseNeg = ntrack->GetZ() +  deltaUNeg * ntrack->GetTgl();
2054          Double_t zphasePos = ptrack->GetZ() +  deltaUPos * ptrack->GetTgl();
2055
2056          Double_t convposz = (zphasePos*negtrackradius+zphaseNeg*postrackradius)/(negtrackradius+postrackradius);
2057
2058          return convposz;
2059 }
2060 //
2061 AliGammaConversionKFVector* AliV0Reader::GetBGGoodV0s(Int_t /*event*/) const{
2062         /*
2063         if(fUseChargedTrackMultiplicityForBG == kTRUE){
2064                 return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),CountESDTracks());
2065         }
2066         else{ // means we use #v0s as multiplicity
2067                 return fBGEventHandler->GetBGGoodV0s(event,fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfGoodV0s);
2068         }
2069         */
2070         return NULL;
2071 }
2072
2073 Int_t AliV0Reader::CountESDTracks(){
2074         // see header file for documentation
2075         if(fNumberOfESDTracks == 0){ // count the good esd tracks
2076                 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2077                         AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
2078                         if(!curTrack){
2079         continue;
2080                         }
2081                         if(fEsdTrackCuts->AcceptTrack(curTrack) ){
2082         fNumberOfESDTracks++;
2083                         }
2084                 }
2085         }
2086
2087         return fNumberOfESDTracks;
2088 }
2089
2090 Bool_t AliV0Reader::CheckIfPi0IsMother(Int_t label){
2091         // see headerfile for documentation
2092         Bool_t iResult=kFALSE;
2093         //      cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
2094         if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 111){
2095                 iResult=kTRUE;
2096         }
2097         return iResult;
2098 }
2099
2100 Bool_t AliV0Reader::CheckIfEtaIsMother(Int_t label){
2101         // see headerfile for documentation
2102         Bool_t iResult=kFALSE;
2103         //      cout<<"Checking particle label, particle is: "<<fMCStack->Particle(TMath::Abs(label))->GetName()<<endl;
2104         if(fMCStack->Particle(TMath::Abs(label))->GetPdgCode() == 221){
2105                 iResult=kTRUE;
2106         }
2107         return iResult;
2108 }
2109
2110
2111
2112
2113 Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliKFParticle* negativeKFParticle, const AliKFParticle * positiveKFParticle, Double_t armenterosQtAlpha[2], Int_t kfProductionMethod){
2114
2115    AliKFParticle PosParticle = *positiveKFParticle;
2116    AliKFParticle NegParticle = *negativeKFParticle;
2117    AliKFParticle Gamma;
2118    if(kfProductionMethod < 3)
2119       Gamma.ConstructGamma(PosParticle, NegParticle);
2120    else if(kfProductionMethod == 3){
2121       Gamma += PosParticle;
2122       Gamma += NegParticle;
2123    }
2124
2125    Double_t VertexGamma[3] = {Gamma.GetX(), Gamma.GetY(), Gamma.GetZ()};
2126    PosParticle.TransportToPoint(VertexGamma);
2127    NegParticle.TransportToPoint(VertexGamma);
2128
2129    AliKFParticle::GetArmenterosPodolanski(PosParticle, NegParticle, armenterosQtAlpha);
2130
2131    return 1;
2132 }
2133
2134
2135
2136
2137 Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliKFParticle* negativeKFParticle, const AliKFParticle * positiveKFParticle, const AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlpha[2] ){
2138         //see header file for documentation
2139
2140         TVector3 momentumVectorPositiveKF(positiveKFParticle->GetPx(),positiveKFParticle->GetPy(),positiveKFParticle->GetPz());
2141         TVector3 momentumVectorNegativeKF(negativeKFParticle->GetPx(),negativeKFParticle->GetPy(),negativeKFParticle->GetPz());
2142         TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
2143
2144         Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2145         Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2146
2147         Float_t alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2148                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2149
2150
2151         Float_t qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2152
2153         armenterosQtAlpha[0]=qt;
2154         armenterosQtAlpha[1]=alfa;
2155
2156         return 1;
2157
2158 }
2159
2160 Bool_t AliV0Reader::GetArmenterosQtAlpha(const AliESDv0* v0, Double_t armenterosQtAlpha[2])
2161 {       //see header file for documentation
2162
2163    Double_t mn[3] = {0,0,0};
2164    Double_t mp[3] = {0,0,0};
2165    Double_t mm[3] = {0,0,0};
2166
2167
2168    Int_t pIndex = 0, nIndex = 0;
2169    pIndex = v0->GetPindex();
2170    nIndex = v0->GetNindex();
2171
2172    AliESDtrack* d[2];
2173    d[0] = dynamic_cast<AliESDtrack*>(fESDEvent->GetTrack(pIndex));
2174    d[1] = dynamic_cast<AliESDtrack*>(fESDEvent->GetTrack(nIndex));
2175
2176    Int_t sign[2];
2177    sign[0] = (int)d[0]->GetSign();
2178    sign[1] = (int)d[1]->GetSign();
2179
2180    Bool_t correct = kFALSE;
2181
2182
2183    if(-1 == sign[0] && 1 == sign[1]){
2184       correct = kFALSE;
2185    }
2186    else{
2187       correct = kTRUE;
2188    }
2189
2190
2191    if(correct){
2192       v0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
2193       v0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
2194    }
2195    else{
2196       v0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
2197       v0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
2198    }
2199    v0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
2200
2201    TVector3 vecN(mn[0],mn[1],mn[2]);
2202    TVector3 vecP(mp[0],mp[1],mp[2]);
2203    TVector3 vecM(mm[0],mm[1],mm[2]);
2204
2205    Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
2206    Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
2207
2208    Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
2209       ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
2210    Double_t qt = vecP.Mag()*sin(thetaP);
2211
2212    armenterosQtAlpha[0]=qt;
2213    armenterosQtAlpha[1]=alfa;
2214
2215    return 1;
2216
2217 }
2218
2219
2220
2221 Bool_t AliV0Reader::GetArmenterosQtAlpha(const TParticle* negativeParticle, const TParticle * positiveParticle, const AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlpha[2] ){
2222         //see header file for documentation
2223
2224         TVector3 momentumVectorPositiveKF(positiveParticle->Px(),positiveParticle->Py(),positiveParticle->Pz());
2225         TVector3 momentumVectorNegativeKF(negativeParticle->Px(),negativeParticle->Py(),negativeParticle->Pz());
2226         TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
2227
2228         Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2229         Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2230
2231         Float_t alfa;
2232         Float_t qt;
2233         if ( positiveParticle->GetPdgCode() == 11 || positiveParticle->GetPdgCode() == 13 || positiveParticle->GetPdgCode() == 15){
2234                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2235                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2236                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2237         } else if ( negativeParticle->GetPdgCode() == -11 || negativeParticle->GetPdgCode() == -13 || negativeParticle->GetPdgCode() == -15){
2238                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2239                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2240                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2241         } else if (positiveParticle->GetPdgCode() < 0 && positiveParticle->GetPdgCode() != -11 && positiveParticle->GetPdgCode() != -13 && positiveParticle->GetPdgCode() != -15 ){
2242                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2243                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2244                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2245         } else {
2246                 alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2247                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2248                 qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2249         }
2250
2251         armenterosQtAlpha[0]=qt;
2252         armenterosQtAlpha[1]=alfa;
2253         return 1;
2254
2255 }
2256
2257 Bool_t AliV0Reader::GetArmenterosQtAlpha(const TParticle* negativeParticle, const TParticle * positiveParticle, const TParticle * gammaCandidate, Double_t armenterosQtAlpha[2] ){
2258         //see header file for documentation
2259
2260         TVector3 momentumVectorPositiveKF(positiveParticle->Px(),positiveParticle->Py(),positiveParticle->Pz());
2261         TVector3 momentumVectorNegativeKF(negativeParticle->Px(),negativeParticle->Py(),negativeParticle->Pz());
2262         TVector3 vecV0(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
2263
2264         Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2265         Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2266
2267         Float_t alfa;
2268         Float_t qt;
2269         if ( positiveParticle->GetPdgCode() == 11 || positiveParticle->GetPdgCode() == 13 || positiveParticle->GetPdgCode() == 15){
2270                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2271                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2272                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2273         } else if ( negativeParticle->GetPdgCode() == -11 || negativeParticle->GetPdgCode() == -13 || negativeParticle->GetPdgCode() == -15){
2274                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2275                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2276                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2277         } else if (positiveParticle->GetPdgCode() < 0 && positiveParticle->GetPdgCode() != -11 && positiveParticle->GetPdgCode() != -13 && positiveParticle->GetPdgCode() != -15 ){
2278                 alfa =((momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)-(momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos))/
2279                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2280                 qt = momentumVectorNegativeKF.Mag()*TMath::Sin(thetaV0neg);
2281         } else {
2282                 alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2283                 ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2284                 qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2285         }
2286
2287         armenterosQtAlpha[0]=qt;
2288         armenterosQtAlpha[1]=alfa;
2289         return 1;
2290
2291 }
2292
2293
2294 Int_t AliV0Reader::GetFirstTPCRow(Double_t radius){
2295
2296
2297         Int_t firstTPCRow=0;
2298         Double_t radiusI        =       84.8;
2299         Double_t radiusO        = 134.6;
2300         Double_t radiusOB = 198.;
2301         Double_t rSizeI  = 0.75;
2302         Double_t rSizeO  = 1.;
2303         Double_t rSizeOB        = 1.5;
2304         Int_t nClsI=63;
2305         Int_t nClsIO=127;
2306
2307         if(radius <= radiusI){
2308                 return firstTPCRow;
2309         }
2310         if(radius>radiusI && radius<=radiusO){
2311                 firstTPCRow = (Int_t)((radius-radiusI)/rSizeI);
2312         }
2313         if(radius>radiusO && radius<=radiusOB){
2314                 firstTPCRow = (Int_t)(nClsI+(radius-radiusO)/rSizeO);
2315         }
2316
2317         if(radius>radiusOB){
2318                 firstTPCRow =(Int_t)(nClsIO+(radius-radiusOB)/rSizeOB);
2319         }
2320
2321
2322         return firstTPCRow;
2323 }
2324 void AliV0Reader::SmearKFParticle(AliKFParticle * kfParticle)
2325 {
2326         Double_t facPBrem = 1.;
2327         Double_t facPSig = 0.;
2328
2329         Double_t phi=0.;
2330         Double_t theta=0.;
2331         Double_t P=0.;
2332
2333         P=kfParticle->GetP();
2334         phi=kfParticle->GetPhi();
2335         if( kfParticle->GetP()!=0){
2336                 theta=acos( kfParticle->Pz()/ kfParticle->GetP());
2337         }
2338
2339         if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
2340                 facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
2341         }
2342
2343         if( fPBremSmearing != 1.){
2344                 if(fBrem!=NULL){
2345                         facPBrem = fBrem->GetRandom();
2346                 }
2347         }
2348
2349         kfParticle->Px() = facPBrem* (1+facPSig)* P*sin(theta)*cos(phi) ;
2350         kfParticle->Py() = facPBrem* (1+facPSig)* P*sin(theta)*sin(phi) ;
2351         kfParticle->Pz() = facPBrem* (1+facPSig)* P*cos(theta) ;
2352         kfParticle->E() = kfParticle->GetP();
2353 }
2354
2355 ///________________________________________________________________________
2356 const AliExternalTrackParam *AliV0Reader::GetExternalTrackParam(AliESDv0 *v0,Int_t charge){
2357
2358     if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
2359
2360     Int_t label;
2361     if(charge>0)label=0;
2362     else label=1;
2363     // Check for sign flip
2364
2365     if(v0){
2366         if(!v0->GetParamN()||!v0->GetParamP())return 0x0;
2367         if(!GetTrack(v0->GetNindex())||!fESDEvent->GetTrack(v0->GetPindex()))return 0x0;
2368         if((GetTrack(v0->GetPindex()))->Charge()==charge){
2369 //          fCurrentTrackLabels[label]=v0->GetPindex();
2370             return v0->GetParamP();}
2371         if((GetTrack(v0->GetNindex()))->Charge()==charge){
2372 //          fCurrentTrackLabels[label]=v0->GetNindex();
2373             return v0->GetParamN();}
2374     }
2375     return 0x0;
2376 }
2377
2378 ///________________________________________________________________________
2379 AliVTrack *AliV0Reader::GetTrack(Int_t label){
2380     if(fESDEvent){
2381                 return (AliESDtrack*)fESDEvent->GetTrack(label);
2382     }
2383 //     if(fAODEvent)return (AliAODTrack*)GetAODTrack(label);
2384         return 0x0;
2385 }
2386
2387 Double_t AliV0Reader::GetV0CosineOfPointingAngle(Double_t V0PointX, Double_t V0PointY, Double_t V0PointZ){
2388    // calculates the pointing angle of the recalculated V0
2389
2390    Double_t momV0[3]; //momentum of the V0
2391    fCurrentV0->GetPxPyPz(momV0[0],momV0[1],momV0[2]);
2392
2393    Double_t PosV0[3]; //Recalculated V0 Position vector
2394
2395    PosV0[0] = V0PointX - fESDEvent->GetPrimaryVertex()->GetX();
2396    PosV0[1] = V0PointY - fESDEvent->GetPrimaryVertex()->GetY();
2397    PosV0[2] = V0PointZ - fESDEvent->GetPrimaryVertex()->GetZ();
2398
2399    Double_t momV02 = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
2400    Double_t PosV02 = PosV0[0]*PosV0[0] + PosV0[1]*PosV0[1] + PosV0[2]*PosV0[2];
2401
2402    Double_t cosinePointingAngle = (PosV0[0]*momV0[0] +  PosV0[1]*momV0[1] + PosV0[2]*momV0[2] ) / TMath::Sqrt(momV02 * PosV02);
2403
2404    return cosinePointingAngle;
2405 }
2406
2407
2408 Double_t AliV0Reader::GetPsiPair(AliESDv0* v0)
2409 {
2410    //
2411    // Angle between daughter momentum plane and plane
2412    //
2413    Float_t magField = fESDEvent->GetMagneticField();
2414
2415    Double_t xyz[3] = {0.,0.,0.};
2416    v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
2417
2418    Double_t mn[3] = {0,0,0};
2419    Double_t mp[3] = {0,0,0};
2420
2421    v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
2422    v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
2423
2424    Double_t deltat = 1.;
2425    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
2426    Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
2427
2428    Double_t momPosProp[3] = {0,0,0};
2429    Double_t momNegProp[3] = {0,0,0};
2430
2431    AliExternalTrackParam nt = *(v0->GetParamN());
2432    AliExternalTrackParam pt = *(v0->GetParamP());
2433
2434    Double_t psiPair = 4.;
2435    if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
2436
2437    if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
2438
2439    pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
2440    nt.GetPxPyPz(momNegProp);
2441
2442   Double_t pEle =
2443      TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
2444
2445   Double_t pPos =
2446      TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
2447
2448   Double_t scalarproduct =
2449      momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
2450
2451   Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
2452
2453   psiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
2454
2455   return psiPair;
2456 }
2457
2458 Int_t AliV0Reader::GetNumberOfHijingPlusPythiaPrimeries(){
2459
2460    // Calculate NPrimaries for LHC11a10b_*
2461
2462    Int_t nproduced = 0;
2463    AliGenCocktailEventHeader *cHeader = dynamic_cast<AliGenCocktailEventHeader*>(fMCEvent->GenEventHeader());
2464    if(cHeader){
2465       TList *genHeaders = cHeader->GetHeaders();
2466       AliGenEventHeader* gh = 0;
2467       for(Int_t i = 0; i<genHeaders->GetEntries();i++){
2468          gh = (AliGenEventHeader*)genHeaders->At(i);
2469          TString GeneratorName = gh->GetName();
2470          //  cout<<i<<"   "<<GeneratorName<<"   "<<gh->NProduced()<<endl;
2471          if(GeneratorName.CompareTo("Hijing") == 0){
2472             nproduced = nproduced + gh->NProduced();
2473          }
2474          else if(GeneratorName.CompareTo("Pythia") == 0){
2475             nproduced = nproduced + gh->NProduced();
2476          }
2477       }
2478    }
2479    if(!cHeader){
2480       nproduced = fMCStack->GetNprimary();
2481    }
2482
2483    //  cout<<fMCStack->GetNprimary()-nproduced<<endl;
2484
2485    return nproduced;
2486 }
2487
2488
2489 Bool_t AliV0Reader::IsParticleFromBGEvent(Int_t index){
2490
2491    Bool_t particleFromBG = kFALSE;
2492
2493    if(index == -1) return kFALSE;
2494    if(index > fNumberOfPrimerisFromHijingAndPythia && index < fMCStack->GetNprimary()){
2495       // cout<<index<<"   "<<fNumberOfPrimerisFromHijingAndPythia<<endl;
2496       return kTRUE;
2497    }
2498    //  else cout<<"Passt Noch "<<index<<"   "<<fNumberOfPrimerisFromHijingAndPythia<<endl;
2499    // cout<<fMCEvent->IsFromBGEvent(index)<<endl;
2500    TParticle *BGParticle = fMCStack->Particle(index);
2501    if(BGParticle->GetMother(0) > -1)  return kFALSE;
2502    Int_t indexMother = fMCStack->Particle(index)->GetMother(0);
2503    particleFromBG = IsParticleFromBGEvent(indexMother);
2504
2505    return kFALSE;
2506 }