PWG4/GammaConv/AliV0Reader.cxx

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