PWG4/GammaConv/AliAnalysisTaskGammaConversion.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.1                                                            *
   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
  23 #include <TChain.h>
  24
  25 // analysis
  26 #include "AliAnalysisTaskGammaConversion.h"
  27 #include "AliStack.h"
  28 #include "AliLog.h"
  29 #include "AliESDtrackCuts.h"
  30 #include "TNtuple.h"
  31 //#include "AliCFManager.h"  // for CF
  32 //#include "AliCFContainer.h"   // for CF
  33 #include "AliGammaConversionAODObject.h"
  34 #include "AliGammaConversionBGHandler.h"
  35
  36 class AliCFContainer;
  37 class AliCFManager;
  38 class AliKFVertex;
  39 class AliAODHandler;
  40 class AliAODEvent;
  41 class ALiESDEvent;
  42 class AliMCEvent;
  43 class AliMCEventHandler;
  44 class AliESDInputHandler;
  45 class AliAnalysisManager;
  46 class Riostream;
  47 class TFile;
  48 class TInterpreter;
  49 class TSystem;
  50 class TROOT;
  51
  52 ClassImp(AliAnalysisTaskGammaConversion)
  53
  54
  55 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
  56 AliAnalysisTaskSE(),
  57   fV0Reader(NULL),
  58   fStack(NULL),
  59   fMCTruth(NULL),    // for CF
  60   fGCMCEvent(NULL),    // for CF
  61   fESDEvent(NULL),
  62   fOutputContainer(NULL),
  63   fCFManager(0x0),   // for CF
  64   fHistograms(NULL),
  65   fTriggerCINT1B(kFALSE),
  66   fDoMCTruth(kFALSE),
  67   fDoNeutralMeson(kFALSE),
  68   fDoJet(kFALSE),
  69   fDoChic(kFALSE),
  70   fKFReconstructedGammasTClone(NULL),
  71   fCurrentEventPosElectronTClone(NULL),
  72   fCurrentEventNegElectronTClone(NULL),
  73   fKFReconstructedGammasCutTClone(NULL),
  74   fPreviousEventTLVNegElectronTClone(NULL),
  75   fPreviousEventTLVPosElectronTClone(NULL),
  76   fElectronv1(),
  77   fElectronv2(),
  78   fElectronMass(-1),
  79   fGammaMass(-1),
  80   fPi0Mass(-1),
  81   fEtaMass(-1),
  82   fGammaWidth(-1),
  83   fPi0Width(-1),
  84   fEtaWidth(-1),
  85   fMinOpeningAngleGhostCut(0.),
  86   fEsdTrackCuts(NULL),
  87   fCalculateBackground(kFALSE),
  88   fWriteNtuple(kFALSE),
  89   fGammaNtuple(NULL),
  90   fNeutralMesonNtuple(NULL),
  91   fTotalNumberOfAddedNtupleEntries(0),
  92   fChargedParticles(NULL),
  93   fChargedParticlesId(),
  94   fGammaPtHighest(0.),
  95   fMinPtForGammaJet(1.),
  96   fMinIsoConeSize(0.2),
  97   fMinPtIsoCone(0.7),
  98   fMinPtGamChargedCorr(0.5),
  99   fMinPtJetCone(0.5),
 100   fLeadingChargedIndex(-1),
 101   fLowPtMapping(1.),
 102   fHighPtMapping(3.),
 103   fDoCF(kFALSE),
 104   fAODBranch(NULL),
 105   fAODBranchName("GammaConv"),
 106   fDoNeutralMesonV0MCCheck(kFALSE),
 107   fKFReconstructedGammasV0Index()
 108 {
 109   // Default constructor
 110
 111   /*   Kenneth: the default constructor should not have any define input/output or the call to SetESDtrackCuts
 112   // Common I/O in slot 0
 113   DefineInput (0, TChain::Class());
 114   DefineOutput(0, TTree::Class());
 115
 116   // Your private output
 117   DefineOutput(1, TList::Class());
 118
 119   // Define standard ESD track cuts for Gamma-hadron correlation
 120   SetESDtrackCuts();
 121   */
 122 }
 123
 124 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
 125   AliAnalysisTaskSE(name),
 126   fV0Reader(NULL),
 127   fStack(NULL),
 128   fMCTruth(NULL),    // for CF
 129   fGCMCEvent(NULL),    // for CF
 130   fESDEvent(NULL),
 131   fOutputContainer(0x0),
 132   fCFManager(0x0),   // for CF
 133   fHistograms(NULL),
 134   fTriggerCINT1B(kFALSE),
 135   fDoMCTruth(kFALSE),
 136   fDoNeutralMeson(kFALSE),
 137   fDoJet(kFALSE),
 138   fDoChic(kFALSE),
 139   fKFReconstructedGammasTClone(NULL),
 140   fCurrentEventPosElectronTClone(NULL),
 141   fCurrentEventNegElectronTClone(NULL),
 142   fKFReconstructedGammasCutTClone(NULL),
 143   fPreviousEventTLVNegElectronTClone(NULL),
 144   fPreviousEventTLVPosElectronTClone(NULL),
 145   fElectronv1(),
 146   fElectronv2(),
 147   fElectronMass(-1),
 148   fGammaMass(-1),
 149   fPi0Mass(-1),
 150   fEtaMass(-1),
 151   fGammaWidth(-1),
 152   fPi0Width(-1),
 153   fEtaWidth(-1),
 154   fMinOpeningAngleGhostCut(0.),
 155   fEsdTrackCuts(NULL),
 156   fCalculateBackground(kFALSE),
 157   fWriteNtuple(kFALSE),
 158   fGammaNtuple(NULL),
 159   fNeutralMesonNtuple(NULL),
 160   fTotalNumberOfAddedNtupleEntries(0),
 161   fChargedParticles(NULL),
 162   fChargedParticlesId(),
 163   fGammaPtHighest(0.),
 164   fMinPtForGammaJet(1.),
 165   fMinIsoConeSize(0.2),
 166   fMinPtIsoCone(0.7),
 167   fMinPtGamChargedCorr(0.5),
 168   fMinPtJetCone(0.5),
 169   fLeadingChargedIndex(-1),
 170   fLowPtMapping(1.),
 171   fHighPtMapping(3.),
 172   fDoCF(kFALSE),
 173   fAODBranch(NULL),
 174   fAODBranchName("GammaConv"),
 175   fDoNeutralMesonV0MCCheck(kFALSE),
 176   fKFReconstructedGammasV0Index()
 177 {
 178   // Common I/O in slot 0
 179   DefineInput (0, TChain::Class());
 180   DefineOutput(0, TTree::Class());
 181
 182   // Your private output
 183   DefineOutput(1, TList::Class());
 184   DefineOutput(2, AliCFContainer::Class());  // for CF
 185
 186
 187   // Define standard ESD track cuts for Gamma-hadron correlation
 188   SetESDtrackCuts();
 189 }
 190
 191 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
 192 {
 193   // Remove all pointers
 194
 195   if(fOutputContainer){
 196     fOutputContainer->Clear() ;
 197     delete fOutputContainer ;
 198   }
 199   if(fHistograms){
 200     delete fHistograms;
 201   }
 202   if(fV0Reader){
 203     delete fV0Reader;
 204   }
 205
 206   // for CF
 207   if(fCFManager){
 208     delete fCFManager;
 209   }
 210
 211   if(fEsdTrackCuts){
 212     delete fEsdTrackCuts;
 213   }
 214
 215   if (fAODBranch) {
 216     fAODBranch->Clear();
 217     delete fAODBranch ;
 218   }
 219 }
 220
 221
 222 void AliAnalysisTaskGammaConversion::Init()
 223 {
 224   // Initialization
 225   // AliLog::SetGlobalLogLevel(AliLog::kError);
 226 }
 227 void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
 228 {
 229   // SetESDtrackCuts
 230   if (fEsdTrackCuts!=NULL){
 231     delete fEsdTrackCuts;
 232   }
 233   fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
 234   //standard cuts from:
 235   //http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
 236
 237   // Cuts used up to 3rd of March
 238
 239  //  fEsdTrackCuts->SetMinNClustersTPC(50);
 240 //   fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 241 //   fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
 242 //   fEsdTrackCuts->SetRequireITSRefit(kTRUE);
 243 //   fEsdTrackCuts->SetMaxNsigmaToVertex(3);
 244 //   fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
 245
 246   //------- To be tested-----------
 247
 248    Double_t minNClustersTPC = 50;
 249    Double_t maxChi2PerClusterTPC = 4.0;
 250    Double_t maxDCAtoVertexXY = 2.4; // cm
 251    Double_t maxDCAtoVertexZ  = 3.2; // cm
 252    fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 253    fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
 254    fEsdTrackCuts->SetRequireITSRefit(kTRUE);
 255    //   fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
 256    fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 257    fEsdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
 258    fEsdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
 259    fEsdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
 260    fEsdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
 261    fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
 262    fEsdTrackCuts->SetEtaRange(-0.9, 0.9);
 263    fEsdTrackCuts->SetPtRange(0.15);
 264
 265
 266
 267
 268   //-----  From Jacek 10.03.03 ------------------/
 269 //     minNClustersTPC = 70;
 270 //     maxChi2PerClusterTPC = 4.0;
 271 //     maxDCAtoVertexXY = 2.4; // cm
 272 //     maxDCAtoVertexZ  = 3.2; // cm
 273
 274 //     esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 275 //     esdTrackCuts->SetRequireTPCRefit(kFALSE);
 276 //     esdTrackCuts->SetRequireTPCStandAlone(kTRUE);
 277 //     esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 278 //     esdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
 279 //     esdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
 280 //     esdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
 281 //     esdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
 282 //     esdTrackCuts->SetDCAToVertex2D(kTRUE);
 283
 284
 285
 286   //  fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 287   //  fV0Reader->SetESDtrackCuts(fEsdTrackCuts);
 288 }
 289
 290 void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
 291 {
 292   // Execute analysis for current event
 293
 294   //  ConnectInputData("");
 295
 296   //Each event needs an empty branch
 297   //  fAODBranch->Clear();
 298   fAODBranch->Delete();
 299
 300   if(fKFReconstructedGammasTClone == NULL){
 301     fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
 302   }
 303   if(fCurrentEventPosElectronTClone == NULL){
 304     fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
 305   }
 306   if(fCurrentEventNegElectronTClone == NULL){
 307     fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
 308   }
 309   if(fKFReconstructedGammasCutTClone == NULL){
 310     fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
 311   }
 312   if(fPreviousEventTLVNegElectronTClone == NULL){
 313     fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
 314   }
 315   if(fPreviousEventTLVPosElectronTClone == NULL){
 316     fPreviousEventTLVPosElectronTClone  = new TClonesArray("TLorentzVector",0);
 317   }
 318   if(fChargedParticles == NULL){
 319     fChargedParticles = new TClonesArray("AliESDtrack",0);
 320   }
 321
 322   //clear TClones
 323   fKFReconstructedGammasTClone->Delete();
 324   fCurrentEventPosElectronTClone->Delete();
 325   fCurrentEventNegElectronTClone->Delete();
 326   fKFReconstructedGammasCutTClone->Delete();
 327   fPreviousEventTLVNegElectronTClone->Delete();
 328   fPreviousEventTLVPosElectronTClone->Delete();
 329
 330   //clear vectors
 331   fElectronv1.clear();
 332   fElectronv2.clear();
 333
 334   fChargedParticles->Delete();
 335
 336   fChargedParticlesId.clear();
 337
 338   fKFReconstructedGammasV0Index.clear();
 339
 340   //Clear the data in the v0Reader
 341   //  fV0Reader->UpdateEventByEventData();
 342
 343   //Take Only events with proper trigger
 344   /*
 345   if(fTriggerCINT1B){
 346     if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
 347   }
 348   */
 349
 350   // Process the MC information
 351   if(fDoMCTruth){
 352     ProcessMCData();
 353   }
 354
 355   //Process the v0 information with no cuts
 356   ProcessV0sNoCut();
 357
 358   // Process the v0 information
 359   ProcessV0s();
 360
 361
 362   //Fill Gamma AOD
 363   FillAODWithConversionGammas() ;
 364
 365   // Process reconstructed gammas
 366   if(fDoNeutralMeson == kTRUE){
 367     ProcessGammasForNeutralMesonAnalysis();
 368   }
 369
 370   if(fDoMCTruth == kTRUE){
 371     CheckV0Efficiency();
 372   }
 373   //Process reconstructed gammas electrons for Chi_c Analysis
 374   if(fDoChic == kTRUE){
 375     ProcessGammaElectronsForChicAnalysis();
 376   }
 377   // Process reconstructed gammas for gamma Jet/hadron correlations
 378   if(fDoJet == kTRUE){
 379     ProcessGammasForGammaJetAnalysis();
 380   }
 381
 382   //calculate background if flag is set
 383   if(fCalculateBackground){
 384     CalculateBackground();
 385   }
 386
 387   //Clear the data in the v0Reader
 388   fV0Reader->UpdateEventByEventData();
 389
 390   PostData(1, fOutputContainer);
 391   PostData(2, fCFManager->GetParticleContainer());  // for CF
 392
 393 }
 394
 395 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t *option){
 396   // see header file for documentation
 397
 398   AliAnalysisTaskSE::ConnectInputData(option);
 399
 400   if(fV0Reader == NULL){
 401     // Write warning here cuts and so on are default if this ever happens
 402   }
 403   fV0Reader->Initialize();
 404   fDoMCTruth = fV0Reader->GetDoMCTruth();
 405 }
 406
 407
 408
 409 void AliAnalysisTaskGammaConversion::ProcessMCData(){
 410   // see header file for documentation
 411
 412   fStack = fV0Reader->GetMCStack();
 413   fMCTruth = fV0Reader->GetMCTruth();  // for CF
 414   fGCMCEvent = fV0Reader->GetMCEvent();  // for CF
 415
 416
 417   // for CF
 418   Double_t containerInput[3];
 419   if(fDoCF){
 420     if(!fGCMCEvent) cout << "NO MC INFO FOUND" << endl;
 421     fCFManager->SetEventInfo(fGCMCEvent);
 422   }
 423   // end for CF
 424
 425
 426   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
 427     return; // aborts if the primary vertex does not have contributors.
 428   }
 429
 430   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
 431     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
 432
 433     if (!particle) {
 434       //print warning here
 435       continue;
 436     }
 437
 438     ///////////////////////Begin Chic Analysis/////////////////////////////
 439
 440     if(particle->GetPdgCode() == 443){//Is JPsi
 441       if(particle->GetNDaughters()==2){
 442         if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
 443            TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
 444
 445           TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
 446           TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
 447           if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
 448             fHistograms->FillTable("Table_Electrons",3);//e+ e-  from J/Psi inside acceptance
 449
 450           if( TMath::Abs(daug0->Eta()) < 0.9){
 451             if(daug0->GetPdgCode() == -11)
 452               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
 453             else
 454               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
 455
 456           }
 457           if(TMath::Abs(daug1->Eta()) < 0.9){
 458             if(daug1->GetPdgCode() == -11)
 459               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
 460             else
 461               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
 462           }
 463         }
 464       }
 465     }
 466     //              const int CHI_C0   = 10441;
 467     //              const int CHI_C1   = 20443;
 468     //              const int CHI_C2   = 445
 469     if(particle->GetPdgCode() == 22){//gamma from JPsi
 470       if(particle->GetMother(0) > -1){
 471         if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
 472            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
 473            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
 474           if(TMath::Abs(particle->Eta()) < 1.2)
 475             fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
 476         }
 477       }
 478     }
 479     if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
 480       if( particle->GetNDaughters() == 2){
 481         TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
 482         TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
 483
 484         if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
 485           if( daug0->GetPdgCode() == 443){
 486             TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
 487             TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
 488             if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
 489               fHistograms->FillTable("Table_Electrons",18);
 490
 491           }//if
 492           else if (daug1->GetPdgCode() == 443){
 493             TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
 494             TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
 495             if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
 496               fHistograms->FillTable("Table_Electrons",18);
 497           }//else if
 498         }//gamma o Jpsi
 499       }//GetNDaughters
 500     }
 501
 502
 503     /////////////////////End Chic Analysis////////////////////////////
 504
 505
 506     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )    continue;
 507
 508     if(particle->R()>fV0Reader->GetMaxRCut())   continue; // cuts on distance from collision point
 509
 510     Double_t tmpPhi=particle->Phi();
 511
 512     if(particle->Phi()> TMath::Pi()){
 513       tmpPhi = particle->Phi()-(2*TMath::Pi());
 514     }
 515
 516     Double_t rapidity;
 517     if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
 518       rapidity=0;
 519     }
 520     else{
 521       rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
 522     }
 523
 524     //process the gammas
 525     if (particle->GetPdgCode() == 22){
 526
 527
 528       if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
 529         continue; // no photon as mothers!
 530       }
 531
 532       if(particle->GetMother(0) >= fStack->GetNprimary()){
 533         continue; // the gamma has a mother, and it is not a primary particle
 534       }
 535
 536       fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
 537       fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
 538       fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
 539       fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
 540       fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
 541
 542       // for CF
 543       if(fDoCF){
 544         containerInput[0] = particle->Pt();
 545         containerInput[1] = particle->Eta();
 546         if(particle->GetMother(0) >=0){
 547           containerInput[2] = fStack->Particle(particle->GetMother(0))->GetMass();
 548         }
 549         else{
 550           containerInput[2]=-1;
 551         }
 552
 553         fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);                                        // generated gamma
 554       }
 555       if(particle->GetMother(0) < 0){   // direct gamma
 556         fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
 557         fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
 558         fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
 559         fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
 560         fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);
 561       }
 562
 563       // looking for conversion (electron + positron from pairbuilding (= 5) )
 564       TParticle* ePos = NULL;
 565       TParticle* eNeg = NULL;
 566
 567       if(particle->GetNDaughters() >= 2){
 568         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
 569           TParticle *tmpDaughter = fStack->Particle(daughterIndex);
 570           if(tmpDaughter->GetUniqueID() == 5){
 571             if(tmpDaughter->GetPdgCode() == 11){
 572               eNeg = tmpDaughter;
 573             }
 574             else if(tmpDaughter->GetPdgCode() == -11){
 575               ePos = tmpDaughter;
 576             }
 577           }
 578         }
 579       }
 580
 581
 582       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
 583         continue;
 584       }
 585
 586
 587       Double_t ePosPhi = ePos->Phi();
 588       if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
 589
 590       Double_t eNegPhi = eNeg->Phi();
 591       if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
 592
 593
 594       if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
 595         continue; // no reconstruction below the Pt cut
 596       }
 597
 598       if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
 599         continue;
 600       }
 601
 602       if(ePos->R()>fV0Reader->GetMaxRCut()){
 603         continue; // cuts on distance from collision point
 604       }
 605
 606       if(TMath::Abs(ePos->Vz()) > fV0Reader->GetMaxZCut()){
 607         continue;   // outside material
 608       }
 609
 610
 611       if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  > ePos->R()){
 612         continue;               // line cut to exclude regions where we do not reconstruct
 613       }
 614
 615
 616       // for CF
 617       if(fDoCF){
 618         fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructable);  // reconstructable gamma
 619       }
 620       fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
 621       fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
 622       fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
 623       fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
 624       fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
 625       fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
 626
 627       fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
 628       fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
 629       fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
 630       fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
 631
 632       fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
 633       fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
 634       fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
 635       fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
 636
 637
 638       // begin Mapping
 639       Int_t rBin    = fHistograms->GetRBin(ePos->R());
 640       Int_t zBin    = fHistograms->GetZBin(ePos->Vz());
 641       Int_t phiBin  = fHistograms->GetPhiBin(particle->Phi());
 642
 643       TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());
 644
 645       TString nameMCMappingPhiR="";
 646       nameMCMappingPhiR.Form("MC_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
 647       // fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
 648
 649       TString nameMCMappingPhi="";
 650       nameMCMappingPhi.Form("MC_Conversion_Mapping_Phi%02d",phiBin);
 651       //      fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
 652       //fHistograms->FillHistogram(nameMCMappingPhi, ePos->Vz(), particle->Eta());
 653
 654       TString nameMCMappingR="";
 655       nameMCMappingR.Form("MC_Conversion_Mapping_R%02d",rBin);
 656       //      fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
 657       //fHistograms->FillHistogram(nameMCMappingR,ePos->Vz(), particle->Eta());
 658
 659       TString nameMCMappingPhiInR="";
 660       nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_in_R_%02d",rBin);
 661       //      fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
 662       fHistograms->FillHistogram(nameMCMappingPhiInR, vtxPos.Phi());
 663
 664       TString nameMCMappingZInR="";
 665       nameMCMappingZInR.Form("MC_Conversion_Mapping_Z_in_R_%02d",rBin);
 666       fHistograms->FillHistogram(nameMCMappingZInR,ePos->Vz() );
 667
 668
 669       TString nameMCMappingPhiInZ="";
 670       nameMCMappingPhiInZ.Form("MC_Conversion_Mapping_Phi_in_Z_%02d",zBin);
 671       //      fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
 672       fHistograms->FillHistogram(nameMCMappingPhiInZ, vtxPos.Phi());
 673
 674       TString nameMCMappingRInZ="";
 675       nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
 676       fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
 677
 678       if(particle->Pt() > fLowPtMapping && particle->Pt()< fHighPtMapping){
 679         TString nameMCMappingMidPtPhiInR="";
 680         nameMCMappingMidPtPhiInR.Form("MC_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
 681         fHistograms->FillHistogram(nameMCMappingMidPtPhiInR, vtxPos.Phi());
 682
 683         TString nameMCMappingMidPtZInR="";
 684         nameMCMappingMidPtZInR.Form("MC_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
 685         fHistograms->FillHistogram(nameMCMappingMidPtZInR,ePos->Vz() );
 686
 687
 688         TString nameMCMappingMidPtPhiInZ="";
 689         nameMCMappingMidPtPhiInZ.Form("MC_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
 690         fHistograms->FillHistogram(nameMCMappingMidPtPhiInZ, vtxPos.Phi());
 691
 692         TString nameMCMappingMidPtRInZ="";
 693         nameMCMappingMidPtRInZ.Form("MC_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
 694         fHistograms->FillHistogram(nameMCMappingMidPtRInZ,ePos->R() );
 695
 696       }
 697
 698       //end mapping
 699
 700       fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
 701       fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
 702       fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
 703       fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
 704       fHistograms->FillHistogram("MC_ConvGamma_E_AsymmetryP",particle->P(),eNeg->P()/particle->P());
 705       fHistograms->FillHistogram("MC_ConvGamma_P_AsymmetryP",particle->P(),ePos->P()/particle->P());
 706
 707
 708       if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
 709         fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
 710         fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
 711         fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
 712         fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
 713         fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
 714
 715       } // end direct gamma
 716       else{   // mother exits
 717         /*      if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0
 718                 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
 719                 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2
 720                 ){
 721                 fMCGammaChic.push_back(particle);
 722                 }
 723         */
 724       }  // end if mother exits
 725     } // end if particle is a photon
 726
 727
 728
 729     // process motherparticles (2 gammas as daughters)
 730     // the motherparticle had already to pass the R and the eta cut, but no line cut.
 731     // the line cut is just valid for the conversions!
 732
 733     if(particle->GetNDaughters() == 2){
 734
 735       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
 736       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
 737
 738       if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
 739
 740       // Check the acceptance for both gammas
 741       Bool_t gammaEtaCut = kTRUE;
 742       if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut()  ) gammaEtaCut = kFALSE;
 743       Bool_t gammaRCut = kTRUE;
 744       if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut()  ) gammaRCut = kFALSE;
 745
 746
 747
 748       // check for conversions now -> have to pass eta, R and line cut!
 749       Bool_t daughter0Electron = kFALSE;
 750       Bool_t daughter0Positron = kFALSE;
 751       Bool_t daughter1Electron = kFALSE;
 752       Bool_t daughter1Positron = kFALSE;
 753
 754       if(daughter0->GetNDaughters() >= 2){  // first gamma
 755         for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
 756           TParticle *tmpDaughter = fStack->Particle(TrackIndex);
 757           if(tmpDaughter->GetUniqueID() == 5){
 758             if(tmpDaughter->GetPdgCode() == 11){
 759               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
 760                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
 761                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
 762                     daughter0Electron = kTRUE;
 763                   }
 764                 }
 765
 766               }
 767             }
 768             else if(tmpDaughter->GetPdgCode() == -11){
 769               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
 770                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
 771                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
 772                     daughter0Positron = kTRUE;
 773                   }
 774                 }
 775               }
 776             }
 777           }
 778         }
 779       }
 780
 781
 782       if(daughter1->GetNDaughters() >= 2){   // second gamma
 783         for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
 784           TParticle *tmpDaughter = fStack->Particle(TrackIndex);
 785           if(tmpDaughter->GetUniqueID() == 5){
 786             if(tmpDaughter->GetPdgCode() == 11){
 787               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
 788                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
 789                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
 790                     daughter1Electron = kTRUE;
 791                   }
 792                 }
 793
 794               }
 795             }
 796             else if(tmpDaughter->GetPdgCode() == -11){
 797               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
 798                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
 799                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
 800                     daughter1Positron = kTRUE;
 801                   }
 802                 }
 803
 804               }
 805
 806             }
 807           }
 808         }
 809       }
 810
 811
 812       if(particle->GetPdgCode()==111){     //Pi0
 813         if( iTracks >= fStack->GetNprimary()){
 814           fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
 815           fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
 816           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
 817           fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
 818           fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
 819           fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
 820           fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
 821           fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
 822           fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
 823
 824           if(gammaEtaCut && gammaRCut){
 825             //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
 826             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
 827             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
 828             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
 829               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
 830               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
 831             }
 832           }
 833         }
 834         else{
 835           fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());
 836           fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
 837           fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
 838           fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
 839           fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
 840           fHistograms->FillHistogram("MC_Pi0_R", particle->R());
 841           fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
 842           fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
 843           fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
 844           if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_Fiducial", particle->Pt());
 845
 846           if(gammaEtaCut && gammaRCut){
 847             //    if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
 848             fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
 849             fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
 850             if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_withinAcceptance_Fiducial", particle->Pt());
 851
 852             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
 853               fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
 854               fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
 855               fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
 856               if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_ConvGamma_withinAcceptance_Fiducial", particle->Pt());
 857             }
 858           }
 859         }
 860       }
 861
 862       if(particle->GetPdgCode()==221){   //Eta
 863         fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
 864         fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
 865         fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
 866         fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
 867         fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
 868         fHistograms->FillHistogram("MC_Eta_R", particle->R());
 869         fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
 870         fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
 871         fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
 872
 873         if(gammaEtaCut && gammaRCut){
 874           //    if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
 875           fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
 876           fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
 877           if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
 878             fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
 879             fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
 880             fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
 881           }
 882
 883         }
 884
 885       }
 886
 887       // all motherparticles with 2 gammas as daughters
 888       fHistograms->FillHistogram("MC_Mother_R", particle->R());
 889       fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
 890       fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
 891       fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
 892       fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
 893       fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
 894       fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
 895       fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
 896       fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
 897       fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
 898       fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());
 899
 900       if(gammaEtaCut && gammaRCut){
 901         //      if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
 902         fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
 903         fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
 904         fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());
 905         if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
 906           fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
 907           fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
 908           fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());
 909
 910         }
 911
 912
 913       }  // end passed R and eta cut
 914
 915     } // end if(particle->GetNDaughters() == 2)
 916
 917   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
 918
 919 } // end ProcessMCData
 920
 921
 922
 923 void AliAnalysisTaskGammaConversion::FillNtuple(){
 924   //Fills the ntuple with the different values
 925
 926   if(fGammaNtuple == NULL){
 927     return;
 928   }
 929   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
 930   for(Int_t i=0;i<numberOfV0s;i++){
 931     Float_t values[27] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 932     AliESDv0 * cV0 = fV0Reader->GetV0(i);
 933     Double_t negPID=0;
 934     Double_t posPID=0;
 935     fV0Reader->GetPIDProbability(negPID,posPID);
 936     values[0]=cV0->GetOnFlyStatus();
 937     values[1]=fV0Reader->CheckForPrimaryVertex();
 938     values[2]=negPID;
 939     values[3]=posPID;
 940     values[4]=fV0Reader->GetX();
 941     values[5]=fV0Reader->GetY();
 942     values[6]=fV0Reader->GetZ();
 943     values[7]=fV0Reader->GetXYRadius();
 944     values[8]=fV0Reader->GetMotherCandidateNDF();
 945     values[9]=fV0Reader->GetMotherCandidateChi2();
 946     values[10]=fV0Reader->GetMotherCandidateEnergy();
 947     values[11]=fV0Reader->GetMotherCandidateEta();
 948     values[12]=fV0Reader->GetMotherCandidatePt();
 949     values[13]=fV0Reader->GetMotherCandidateMass();
 950     values[14]=fV0Reader->GetMotherCandidateWidth();
 951     //    values[15]=fV0Reader->GetMotherMCParticle()->Pt();   MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
 952     values[16]=fV0Reader->GetOpeningAngle();
 953     values[17]=fV0Reader->GetNegativeTrackEnergy();
 954     values[18]=fV0Reader->GetNegativeTrackPt();
 955     values[19]=fV0Reader->GetNegativeTrackEta();
 956     values[20]=fV0Reader->GetNegativeTrackPhi();
 957     values[21]=fV0Reader->GetPositiveTrackEnergy();
 958     values[22]=fV0Reader->GetPositiveTrackPt();
 959     values[23]=fV0Reader->GetPositiveTrackEta();
 960     values[24]=fV0Reader->GetPositiveTrackPhi();
 961     values[25]=fV0Reader->HasSameMCMother();
 962     if(values[25] != 0){
 963       values[26]=fV0Reader->GetMotherMCParticlePDGCode();
 964       values[15]=fV0Reader->GetMotherMCParticle()->Pt();
 965     }
 966     fTotalNumberOfAddedNtupleEntries++;
 967     fGammaNtuple->Fill(values);
 968   }
 969   fV0Reader->ResetV0IndexNumber();
 970
 971 }
 972
 973 void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
 974   // Process all the V0's without applying any cuts to it
 975
 976   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
 977   for(Int_t i=0;i<numberOfV0s;i++){
 978     /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
 979
 980     if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
 981       continue;
 982     }
 983
 984     //    if( !fV0Reader->GetV0(i)->GetOnFlyStatus()){
 985     if( !fV0Reader->CheckV0FinderStatus(i)){
 986       continue;
 987     }
 988
 989
 990     if( !((fV0Reader->GetNegativeESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ||
 991         !((fV0Reader->GetPositiveESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ){
 992       continue;
 993     }
 994
 995     if( fV0Reader->GetNegativeESDTrack()->GetSign()== fV0Reader->GetPositiveESDTrack()->GetSign()){
 996       continue;
 997     }
 998
 999     if( fV0Reader->GetNegativeESDTrack()->GetKinkIndex(0) > 0 ||
1000         fV0Reader->GetPositiveESDTrack()->GetKinkIndex(0) > 0) {
1001       continue;
1002     }
1003
1004     if(fDoMCTruth){
1005
1006       if(fV0Reader->HasSameMCMother() == kFALSE){
1007         continue;
1008       }
1009
1010       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1011       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1012
1013       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1014         continue;
1015       }
1016       if(negativeMC->GetPdgCode() == positiveMC->GetPdgCode()){
1017         continue;
1018       }
1019
1020       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
1021         continue;
1022       }
1023
1024       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1025
1026         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1027         fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1028         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1029         fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1030         fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1031         fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1032         fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1033         fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1034         fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1035         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1036
1037         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1038         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1039
1040         fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1041         fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
1042         fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1043         fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1044         fHistograms->FillHistogram("ESD_NoCutConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1045         fHistograms->FillHistogram("ESD_NoCutConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1046         fHistograms->FillHistogram("ESD_NoCutConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1047         fHistograms->FillHistogram("ESD_NoCutConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1048
1049         fHistograms->FillHistogram("ESD_NoCutConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1050         fHistograms->FillHistogram("ESD_NoCutConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1051         fHistograms->FillHistogram("ESD_NoCutConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1052         fHistograms->FillHistogram("ESD_NoCutConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1053
1054         //store MCTruth properties
1055         fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1056         fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1057         fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1058       }
1059     }
1060   }
1061   fV0Reader->ResetV0IndexNumber();
1062 }
1063
1064 void AliAnalysisTaskGammaConversion::ProcessV0s(){
1065   // see header file for documentation
1066
1067
1068   if(fWriteNtuple == kTRUE){
1069     FillNtuple();
1070   }
1071
1072   Int_t nSurvivingV0s=0;
1073   while(fV0Reader->NextV0()){
1074     nSurvivingV0s++;
1075
1076
1077     //-------------------------- filling v0 information -------------------------------------
1078     fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
1079     fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1080     fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1081     fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1082
1083     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
1084     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
1085     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
1086     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
1087     fHistograms->FillHistogram("ESD_E_nTPCClusters", fV0Reader->GetNegativeTracknTPCClusters());
1088     fHistograms->FillHistogram("ESD_E_nITSClusters", fV0Reader->GetNegativeTracknITSClusters());
1089
1090     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
1091     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
1092     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
1093     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
1094     fHistograms->FillHistogram("ESD_P_nTPCClusters", fV0Reader->GetPositiveTracknTPCClusters());
1095     fHistograms->FillHistogram("ESD_P_nITSClusters", fV0Reader->GetPositiveTracknITSClusters());
1096
1097     fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1098     fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1099     fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1100     fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1101     fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1102     fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1103     fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1104     fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1105     fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1106     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1107
1108     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1109     fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1110
1111     fHistograms->FillHistogram("ESD_ConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1112     fHistograms->FillHistogram("ESD_ConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1113     fHistograms->FillHistogram("ESD_ConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1114     fHistograms->FillHistogram("ESD_ConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1115
1116     fHistograms->FillHistogram("ESD_ConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1117     fHistograms->FillHistogram("ESD_ConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1118     fHistograms->FillHistogram("ESD_ConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1119     fHistograms->FillHistogram("ESD_ConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1120
1121
1122     // begin mapping
1123     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());
1124     Int_t zBin    = fHistograms->GetZBin(fV0Reader->GetZ());
1125     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
1126     TVector3 vtxConv(fV0Reader->GetX(),fV0Reader->GetY(), fV0Reader->GetZ());
1127
1128     //    Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
1129
1130     TString nameESDMappingPhiR="";
1131     nameESDMappingPhiR.Form("ESD_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
1132     //fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
1133
1134     TString nameESDMappingPhi="";
1135     nameESDMappingPhi.Form("ESD_Conversion_Mapping_Phi%02d",phiBin);
1136     //fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
1137
1138     TString nameESDMappingR="";
1139     nameESDMappingR.Form("ESD_Conversion_Mapping_R%02d",rBin);
1140     //fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);
1141
1142     TString nameESDMappingPhiInR="";
1143     nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_in_R_%02d",rBin);
1144     //    fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1145     fHistograms->FillHistogram(nameESDMappingPhiInR, vtxConv.Phi());
1146
1147     TString nameESDMappingZInR="";
1148     nameESDMappingZInR.Form("ESD_Conversion_Mapping_Z_in_R_%02d",rBin);
1149     fHistograms->FillHistogram(nameESDMappingZInR, fV0Reader->GetZ());
1150
1151     TString nameESDMappingPhiInZ="";
1152     nameESDMappingPhiInZ.Form("ESD_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1153     //    fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1154     fHistograms->FillHistogram(nameESDMappingPhiInZ, vtxConv.Phi());
1155
1156     TString nameESDMappingRInZ="";
1157     nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
1158     fHistograms->FillHistogram(nameESDMappingRInZ, fV0Reader->GetXYRadius());
1159
1160     if(fV0Reader->GetMotherCandidatePt() > fLowPtMapping && fV0Reader->GetMotherCandidatePt()< fHighPtMapping){
1161       TString nameESDMappingMidPtPhiInR="";
1162       nameESDMappingMidPtPhiInR.Form("ESD_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1163       fHistograms->FillHistogram(nameESDMappingMidPtPhiInR, vtxConv.Phi());
1164
1165       TString nameESDMappingMidPtZInR="";
1166       nameESDMappingMidPtZInR.Form("ESD_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1167       fHistograms->FillHistogram(nameESDMappingMidPtZInR, fV0Reader->GetZ());
1168
1169       TString nameESDMappingMidPtPhiInZ="";
1170       nameESDMappingMidPtPhiInZ.Form("ESD_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1171       fHistograms->FillHistogram(nameESDMappingMidPtPhiInZ, vtxConv.Phi());
1172
1173       TString nameESDMappingMidPtRInZ="";
1174       nameESDMappingMidPtRInZ.Form("ESD_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1175       fHistograms->FillHistogram(nameESDMappingMidPtRInZ, fV0Reader->GetXYRadius());
1176
1177     }
1178
1179
1180     // end mapping
1181
1182     new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()])  AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
1183     fKFReconstructedGammasV0Index.push_back(fV0Reader->GetCurrentV0IndexNumber()-1);
1184     //    fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
1185     fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
1186     fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
1187
1188     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
1189     if(fDoMCTruth){
1190
1191       if(fV0Reader->HasSameMCMother() == kFALSE){
1192         continue;
1193       }
1194       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1195       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1196
1197       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1198         continue;
1199       }
1200
1201       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1202         continue;
1203       }
1204       if(negativeMC->GetUniqueID() == 4 && positiveMC->GetUniqueID() ==4){// fill r distribution for Dalitz decays
1205         if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 111){ //pi0
1206           fHistograms->FillHistogram("ESD_TrueDalitzContamination_R", fV0Reader->GetXYRadius());
1207         }
1208       }
1209
1210       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
1211         continue;
1212       }
1213
1214       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1215
1216         if(fDoCF){
1217           Double_t containerInput[3];
1218           containerInput[0] = fV0Reader->GetMotherCandidatePt();
1219           containerInput[1] = fV0Reader->GetMotherCandidateEta();
1220           containerInput[2] = fV0Reader->GetMotherCandidateMass();
1221           fCFManager->GetParticleContainer()->Fill(containerInput,kStepTrueGamma); // for CF
1222         }
1223
1224         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1225         fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1226         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1227         fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1228         fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1229         fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1230         fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1231         fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1232         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1233         fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1234         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
1235         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
1236         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1237         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1238
1239         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1240         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1241
1242
1243         fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1244         fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
1245         fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1246         fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1247
1248         fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1249         fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1250         fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1251         fHistograms->FillHistogram("ESD_TrueConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1252
1253         fHistograms->FillHistogram("ESD_TrueConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1254         fHistograms->FillHistogram("ESD_TrueConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1255         fHistograms->FillHistogram("ESD_TrueConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1256         fHistograms->FillHistogram("ESD_TrueConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1257
1258
1259
1260         //store MCTruth properties
1261         fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1262         fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1263         fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1264
1265         //resolution
1266         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();
1267         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();
1268         Double_t resdPt = 0;
1269         if(mcpt > 0){
1270           resdPt = ((esdpt - mcpt)/mcpt)*100;
1271         }
1272         else if(mcpt < 0){
1273           cout<<"Pt of MC particle is negative, this will cause wrong calculation of resPt"<<endl;
1274         }
1275
1276         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);
1277         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
1278         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
1279
1280         Double_t resdZ = 0;
1281         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
1282           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;
1283         }
1284         Double_t resdZAbs = 0;
1285         resdZAbs = (fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz());
1286
1287         fHistograms->FillHistogram("Resolution_dZAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdZAbs);
1288         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
1289         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
1290         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
1291
1292         Double_t resdR = 0;
1293         if(fV0Reader->GetNegativeMCParticle()->R() != 0){
1294           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;
1295         }
1296         Double_t resdRAbs = 0;
1297         resdRAbs = (fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R());
1298
1299         fHistograms->FillHistogram("Resolution_dRAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdRAbs);
1300         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
1301         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
1302         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
1303         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);
1304
1305         Double_t resdPhiAbs=0;
1306         resdPhiAbs=0;
1307         resdPhiAbs= (fV0Reader->GetMotherCandidatePhi()-fV0Reader->GetNegativeMCParticle()->Phi());
1308         fHistograms->FillHistogram("Resolution_dPhiAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdPhiAbs);
1309
1310       }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
1311     }//if(fDoMCTruth)
1312   }//while(fV0Reader->NextV0)
1313   fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
1314   fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
1315   fHistograms->FillHistogram("ESD_NumberOfContributorsVtx", fV0Reader->GetNumberOfContributorsVtx());
1316
1317   fV0Reader->ResetV0IndexNumber();
1318
1319 }
1320
1321 void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
1322   // Fill AOD with reconstructed Gamma
1323
1324   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
1325     //  for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
1326     //Create AOD particle object from AliKFParticle
1327
1328     /*    AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
1329     //You could add directly AliKFParticle objects to the AOD, avoiding dependences with PartCorr
1330     //but this means that I have to work a little bit more in my side.
1331     //AODPWG4Particle objects are simpler and lighter, I think
1332     AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
1333     gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
1334     gamma.SetCaloLabel(-1,-1); //How to get the MC label of the 2 electrons that form the gamma?
1335     gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
1336     gamma.SetPdg(AliCaloPID::kPhotonConv); //photon id
1337     gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
1338
1339     //Add it to the aod list
1340     Int_t i = fAODBranch->GetEntriesFast();
1341     new((*fAODBranch)[i])  AliAODPWG4Particle(gamma);
1342     */
1343     //    AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
1344     AliKFParticle * gammakf = (AliKFParticle *)fKFReconstructedGammasTClone->At(gammaIndex);
1345     AliGammaConversionAODObject aodObject;
1346     aodObject.SetPx(gammakf->GetPx());
1347     aodObject.SetPy(gammakf->GetPy());
1348     aodObject.SetPz(gammakf->GetPz());
1349     aodObject.SetLabel1(fElectronv1[gammaIndex]);
1350     aodObject.SetLabel2(fElectronv2[gammaIndex]);
1351     Int_t i = fAODBranch->GetEntriesFast();
1352     new((*fAODBranch)[i])  AliGammaConversionAODObject(aodObject);
1353   }
1354
1355 }
1356
1357
1358 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
1359   // see header file for documentation
1360
1361   //  for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
1362   //    for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
1363
1364   if(fKFReconstructedGammasTClone->GetEntriesFast()>fV0Reader->GetNumberOfV0s()){
1365     cout<<"Warning, number of entries in the tclone is bigger than number of v0s"<<endl;
1366   }
1367
1368   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1369     for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
1370
1371       //      AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
1372       //      AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
1373
1374       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1375       AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
1376
1377       if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1378         continue;
1379       }
1380       if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1381         continue;
1382       }
1383
1384       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
1385
1386       Double_t massTwoGammaCandidate = 0.;
1387       Double_t widthTwoGammaCandidate = 0.;
1388       Double_t chi2TwoGammaCandidate =10000.;
1389       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
1390       if(twoGammaCandidate->GetNDF()>0){
1391         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
1392
1393         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){
1394
1395           TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
1396           TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
1397
1398           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
1399           Double_t rapidity;
1400           if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){
1401             rapidity=0;
1402           }
1403           else{
1404             rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
1405           }
1406
1407           if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut){
1408             delete twoGammaCandidate;
1409             continue;   // minimum opening angle to avoid using ghosttracks
1410           }
1411
1412           fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
1413           fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
1414           fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
1415           fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
1416           fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);
1417           fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
1418           fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
1419           fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!
1420           fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
1421           fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
1422           fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1423           fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
1424
1425           //      if(fDoNeutralMesonV0MCCheck){
1426           if(fDoMCTruth){
1427             //Kenneth: Checking the eta of the gamma to check the difference between 0.9 and 1.2
1428             Int_t indexKF1 = fKFReconstructedGammasV0Index.at(firstGammaIndex);
1429             if(indexKF1<fV0Reader->GetNumberOfV0s()){
1430               fV0Reader->GetV0(indexKF1);//updates to the correct v0
1431               Double_t eta1 = fV0Reader->GetMotherCandidateEta();
1432               Bool_t isRealPi0=kFALSE;
1433               Int_t gamma1MotherLabel=-1;
1434               if(fV0Reader->HasSameMCMother() == kTRUE){
1435                 //cout<<"This v0 is a real v0!!!!"<<endl;
1436                 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1437                 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1438                 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1439                   if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
1440                     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1441                       gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
1442                     }
1443                   }
1444                 }
1445               }
1446               Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
1447               if(indexKF1 == indexKF2){
1448                 cout<<"index of the two KF particles are the same.... should not happen"<<endl;
1449               }
1450               if(indexKF2<fV0Reader->GetNumberOfV0s()){
1451                 fV0Reader->GetV0(indexKF2);
1452                 Double_t eta2 = fV0Reader->GetMotherCandidateEta();
1453                 Int_t gamma2MotherLabel=-1;
1454                 if(fV0Reader->HasSameMCMother() == kTRUE){
1455                   TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1456                   TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1457                   if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1458                     if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
1459                       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1460                         gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
1461                       }
1462                     }
1463                   }
1464                 }
1465                 if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
1466                   if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
1467                     isRealPi0=kTRUE;
1468                   }
1469                 }
1470
1471                 if(TMath::Abs(eta1)>0.9 && TMath::Abs(eta2)>0.9){
1472                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
1473                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1474                   if(isRealPi0){
1475                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_1212",massTwoGammaCandidate);
1476                     fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_1212",openingAngleTwoGammaCandidate);
1477                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1478                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1479                     fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
1480                   }
1481                 }
1482                 else if(TMath::Abs(eta1)>0.9 || TMath::Abs(eta2)>0.9){
1483                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
1484                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1485                   if(isRealPi0){
1486                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_0912",massTwoGammaCandidate);
1487                     fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0912",openingAngleTwoGammaCandidate);
1488                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1489                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1490                     fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
1491                   }
1492                 }
1493                 else{
1494                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
1495                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1496                   if(isRealPi0){
1497                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_0909",massTwoGammaCandidate);
1498                     fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0909",openingAngleTwoGammaCandidate);
1499                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1500                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1501                     fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
1502                   }
1503                 }
1504               }
1505             }
1506           }
1507           if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 &&  TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
1508             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1509             fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
1510           }
1511
1512           if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
1513             fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
1514             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1515           }
1516           else if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 || TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
1517             fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
1518             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1519           }
1520           else{
1521             fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
1522             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1523           }
1524         }
1525       }
1526       delete twoGammaCandidate;
1527     }
1528   }
1529 }
1530
1531 void AliAnalysisTaskGammaConversion::CalculateBackground(){
1532   // see header file for documentation
1533
1534
1535   TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
1536
1537   for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
1538     AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
1539     for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
1540       AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
1541       for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1542         AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
1543
1544         AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
1545
1546         Double_t massBG =0.;
1547         Double_t widthBG = 0.;
1548         Double_t chi2BG =10000.;
1549         backgroundCandidate->GetMass(massBG,widthBG);
1550         if(backgroundCandidate->GetNDF()>0){
1551           chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
1552           if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){
1553
1554             TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
1555             TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
1556
1557             Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
1558
1559             Double_t rapidity;
1560             if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
1561             else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
1562
1563
1564
1565
1566             if(openingAngleBG < fMinOpeningAngleGhostCut ){
1567               delete backgroundCandidate;
1568               continue;   // minimum opening angle to avoid using ghosttracks
1569             }
1570
1571             // original
1572             fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
1573             fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
1574             fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
1575             fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
1576             fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
1577             fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
1578             fHistograms->FillHistogram("ESD_Background_Mass", massBG);
1579             fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
1580             fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
1581             fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
1582             fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
1583             fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
1584
1585             if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
1586               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
1587               fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
1588             }
1589
1590
1591             // test
1592             AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
1593
1594             Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
1595             Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1596
1597             fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
1598             fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
1599             fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
1600             fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
1601             fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
1602             fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
1603             fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
1604             fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
1605             fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
1606             fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
1607             fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
1608             fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
1609
1610             if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
1611               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
1612               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
1613             }
1614           }
1615         }
1616         delete backgroundCandidate;
1617       }
1618     }
1619   }
1620 }
1621
1622
1623 void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
1624   //ProcessGammasForGammaJetAnalysis
1625
1626   Double_t distIsoMin;
1627
1628   CreateListOfChargedParticles();
1629
1630
1631   //  for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
1632   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
1633     AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
1634     TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
1635
1636     if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
1637       distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
1638
1639       if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
1640         CalculateJetCone(gammaIndex);
1641       }
1642     }
1643   }
1644 }
1645
1646 void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
1647   // CreateListOfChargedParticles
1648
1649   fESDEvent = fV0Reader->GetESDEvent();
1650   Int_t numberOfESDTracks=0;
1651   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1652     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1653
1654     if(!curTrack){
1655       continue;
1656     }
1657
1658     if(fEsdTrackCuts->AcceptTrack(curTrack) ){
1659       new((*fChargedParticles)[fChargedParticles->GetEntriesFast()])  AliESDtrack(*curTrack);
1660       //      fChargedParticles.push_back(curTrack);
1661       fChargedParticlesId.push_back(iTracks);
1662       numberOfESDTracks++;
1663     }
1664   }
1665   fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
1666 }
1667 void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
1668   // CaculateJetCone
1669
1670   Double_t cone;
1671   Double_t coneSize=0.3;
1672   Double_t ptJet=0;
1673
1674   //  AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
1675   AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
1676
1677   TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
1678
1679   AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
1680
1681   Double_t momLeadingCharged[3];
1682   leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
1683
1684   TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
1685
1686   Double_t phi1=momentumVectorLeadingCharged.Phi();
1687   Double_t eta1=momentumVectorLeadingCharged.Eta();
1688   Double_t phi3=momentumVectorCurrentGamma.Phi();
1689
1690   for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
1691     AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
1692     Int_t chId = fChargedParticlesId[iCh];
1693     if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
1694     Double_t mom[3];
1695     curTrack->GetConstrainedPxPyPz(mom);
1696     TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
1697     Double_t phi2=momentumVectorChargedParticle.Phi();
1698     Double_t eta2=momentumVectorChargedParticle.Eta();
1699
1700
1701     cone=100.;
1702     if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
1703       cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
1704     }else{
1705       if( (phi2 - phi1)> TMath::TwoPi()-coneSize ){
1706         cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-TMath::TwoPi()-phi1),2) );
1707       }
1708       if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
1709         cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
1710       }
1711     }
1712
1713     if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
1714       ptJet+= momentumVectorChargedParticle.Pt();
1715       Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
1716       Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
1717       fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
1718       fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
1719
1720     }
1721
1722     Double_t dphiHdrGam=phi3-phi2;
1723     if ( dphiHdrGam < (-TMath::PiOver2())){
1724       dphiHdrGam+=(TMath::TwoPi());
1725     }
1726
1727     if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
1728       dphiHdrGam-=(TMath::TwoPi());
1729     }
1730
1731     if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
1732       fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
1733     }
1734   }//track loop
1735
1736
1737 }
1738
1739 Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
1740   // GetMinimumDistanceToCharge
1741
1742   Double_t fIsoMin=100.;
1743   Double_t ptLeadingCharged=-1.;
1744
1745   fLeadingChargedIndex=-1;
1746
1747   AliKFParticle * gammaHighestPt = (AliKFParticle*)fKFReconstructedGammasTClone->At(indexHighestPtGamma);
1748   TVector3 momentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
1749
1750   Double_t phi1=momentumVectorgammaHighestPt.Phi();
1751   Double_t eta1=momentumVectorgammaHighestPt.Eta();
1752
1753   for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
1754     AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
1755     Int_t chId = fChargedParticlesId[iCh];
1756     if(fElectronv1[indexHighestPtGamma]==chId || fElectronv2[indexHighestPtGamma]==chId) continue;
1757     Double_t mom[3];
1758     curTrack->GetConstrainedPxPyPz(mom);
1759     TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
1760     Double_t phi2=momentumVectorChargedParticle.Phi();
1761     Double_t eta2=momentumVectorChargedParticle.Eta();
1762     Double_t iso=pow(  (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
1763
1764     if(momentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
1765       if (iso<fIsoMin){
1766         fIsoMin=iso;
1767       }
1768     }
1769
1770     Double_t dphiHdrGam=phi1-phi2;
1771     if ( dphiHdrGam < (-TMath::PiOver2())){
1772       dphiHdrGam+=(TMath::TwoPi());
1773     }
1774
1775     if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
1776       dphiHdrGam-=(TMath::TwoPi());
1777     }
1778     if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
1779       fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
1780     }
1781
1782     if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
1783       if (momentumVectorChargedParticle.Pt()> ptLeadingCharged && momentumVectorChargedParticle.Pt()>0.1*momentumVectorgammaHighestPt.Pt()){
1784         ptLeadingCharged=momentumVectorChargedParticle.Pt();
1785         fLeadingChargedIndex=iCh;
1786       }
1787     }
1788
1789   }//track loop
1790   fHistograms->FillHistogram("ESD_MinimumIsoDistance",fIsoMin);
1791   return fIsoMin;
1792
1793 }
1794
1795 Int_t  AliAnalysisTaskGammaConversion::GetIndexHighestPtGamma(){
1796   //GetIndexHighestPtGamma
1797
1798   Int_t indexHighestPtGamma=-1;
1799   //Double_t
1800   fGammaPtHighest = -100.;
1801
1802   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1803     AliKFParticle * gammaHighestPtCandidate = (AliKFParticle*)fKFReconstructedGammasTClone->At(firstGammaIndex);
1804     TVector3 momentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
1805     if (momentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
1806       fGammaPtHighest=momentumVectorgammaHighestPtCandidate.Pt();
1807       //gammaHighestPt = gammaHighestPtCandidate;
1808       indexHighestPtGamma=firstGammaIndex;
1809     }
1810   }
1811
1812   return indexHighestPtGamma;
1813
1814 }
1815
1816
1817 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
1818 {
1819   // Terminate analysis
1820   //
1821   AliDebug(1,"Do nothing in Terminate");
1822 }
1823
1824 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
1825 {
1826   //AOD
1827   if(fAODBranch==NULL){
1828     fAODBranch = new TClonesArray("AliGammaConversionAODObject", 0);
1829   }
1830   fAODBranch->Delete();
1831   fAODBranch->SetName(fAODBranchName);
1832   AddAODBranch("TClonesArray", &fAODBranch);
1833
1834   // Create the output container
1835   if(fOutputContainer != NULL){
1836     delete fOutputContainer;
1837     fOutputContainer = NULL;
1838   }
1839   if(fOutputContainer == NULL){
1840     fOutputContainer = new TList();
1841     fOutputContainer->SetOwner(kTRUE);
1842   }
1843
1844   //Adding the histograms to the output container
1845   fHistograms->GetOutputContainer(fOutputContainer);
1846
1847
1848   if(fWriteNtuple){
1849     if(fGammaNtuple == NULL){
1850       fGammaNtuple = new TNtuple("V0ntuple","V0ntuple","OnTheFly:HasVertex:NegPIDProb:PosPIDProb:X:Y:Z:R:MotherCandidateNDF:MotherCandidateChi2:MotherCandidateEnergy:MotherCandidateEta:MotherCandidatePt:MotherCandidateMass:MotherCandidateWidth:MCMotherCandidatePT:EPOpeningAngle:ElectronEnergy:ElectronPt:ElectronEta:ElectronPhi:PositronEnergy:PositronPt:PositronEta:PositronPhi:HasSameMCMother:MotherMCParticlePIDCode",50000);
1851     }
1852     if(fNeutralMesonNtuple == NULL){
1853       fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
1854     }
1855     TList * ntupleTList = new TList();
1856     ntupleTList->SetOwner(kTRUE);
1857     ntupleTList->SetName("Ntuple");
1858     ntupleTList->Add((TNtuple*)fGammaNtuple);
1859     fOutputContainer->Add(ntupleTList);
1860   }
1861
1862   fOutputContainer->SetName(GetName());
1863 }
1864
1865 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{
1866   //helper function
1867   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
1868   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
1869   return v3D0.Angle(v3D1);
1870 }
1871
1872 void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
1873   // see header file for documentation
1874
1875   vector<Int_t> indexOfGammaParticle;
1876
1877   fStack = fV0Reader->GetMCStack();
1878
1879   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
1880     return; // aborts if the primary vertex does not have contributors.
1881   }
1882
1883   for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
1884     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
1885     if(particle->GetPdgCode()==22){     //Gamma
1886       if(particle->GetNDaughters() >= 2){
1887         TParticle* electron=NULL;
1888         TParticle* positron=NULL;
1889         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
1890           TParticle *tmpDaughter = fStack->Particle(daughterIndex);
1891           if(tmpDaughter->GetUniqueID() == 5){
1892             if(tmpDaughter->GetPdgCode() == 11){
1893               electron = tmpDaughter;
1894             }
1895             else if(tmpDaughter->GetPdgCode() == -11){
1896               positron = tmpDaughter;
1897             }
1898           }
1899         }
1900         if(electron!=NULL && positron!=0){
1901           if(electron->R()<160){
1902             indexOfGammaParticle.push_back(iTracks);
1903           }
1904         }
1905       }
1906     }
1907   }
1908
1909   Int_t nFoundGammas=0;
1910   Int_t nNotFoundGammas=0;
1911
1912   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1913   for(Int_t i=0;i<numberOfV0s;i++){
1914     fV0Reader->GetV0(i);
1915
1916     if(fV0Reader->HasSameMCMother() == kFALSE){
1917       continue;
1918     }
1919
1920     TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1921     TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1922
1923     if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1924       continue;
1925     }
1926     if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1927       continue;
1928     }
1929
1930     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1931       //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
1932       for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
1933         if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
1934           nFoundGammas++;
1935         }
1936         else{
1937           nNotFoundGammas++;
1938         }
1939       }
1940     }
1941   }
1942 }
1943
1944
1945 void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
1946   // see header file for documantation
1947
1948   fESDEvent = fV0Reader->GetESDEvent();
1949
1950
1951   TClonesArray * vESDeNegTemp = new TClonesArray("AliESDtrack",0);
1952   TClonesArray * vESDePosTemp = new TClonesArray("AliESDtrack",0);
1953   TClonesArray * vESDxNegTemp = new TClonesArray("AliESDtrack",0);
1954   TClonesArray * vESDxPosTemp = new TClonesArray("AliESDtrack",0);
1955   TClonesArray * vESDeNegNoJPsi = new TClonesArray("AliESDtrack",0);
1956   TClonesArray * vESDePosNoJPsi = new TClonesArray("AliESDtrack",0);
1957
1958   /*
1959     vector <AliESDtrack*> vESDeNegTemp(0);
1960     vector <AliESDtrack*> vESDePosTemp(0);
1961     vector <AliESDtrack*> vESDxNegTemp(0);
1962     vector <AliESDtrack*> vESDxPosTemp(0);
1963     vector <AliESDtrack*> vESDeNegNoJPsi(0);
1964     vector <AliESDtrack*> vESDePosNoJPsi(0);
1965   */
1966
1967
1968   fHistograms->FillTable("Table_Electrons",0);//Count number of Events
1969
1970   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1971     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1972
1973     if(!curTrack){
1974       //print warning here
1975       continue;
1976     }
1977
1978     double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
1979     double r[3];curTrack->GetConstrainedXYZ(r);
1980
1981     TVector3 rXYZ(r);
1982
1983     fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
1984
1985     Bool_t flagKink       =  kTRUE;
1986     Bool_t flagTPCrefit   =  kTRUE;
1987     Bool_t flagTRDrefit   =  kTRUE;
1988     Bool_t flagITSrefit   =  kTRUE;
1989     Bool_t flagTRDout     =  kTRUE;
1990     Bool_t flagVertex     =  kTRUE;
1991
1992
1993     //Cuts ---------------------------------------------------------------
1994
1995     if(curTrack->GetKinkIndex(0) > 0){
1996       fHistograms->FillHistogram("Table_Electrons",5);//Count kink
1997       flagKink = kFALSE;
1998     }
1999
2000     ULong_t trkStatus = curTrack->GetStatus();
2001
2002     ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
2003
2004     if(!tpcRefit){
2005       fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
2006       flagTPCrefit = kFALSE;
2007     }
2008
2009     ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
2010     if(!itsRefit){
2011       fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
2012       flagITSrefit = kFALSE;
2013     }
2014
2015     ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
2016
2017     if(!trdRefit){
2018       fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
2019       flagTRDrefit = kFALSE;
2020     }
2021
2022     ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
2023
2024     if(!trdOut) {
2025       fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
2026       flagTRDout = kFALSE;
2027     }
2028
2029     double nsigmaToVxt = GetSigmaToVertex(curTrack);
2030
2031     if(nsigmaToVxt > 3){
2032       fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
2033       flagVertex = kFALSE;
2034     }
2035
2036     if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
2037     fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
2038
2039
2040     Stat_t pid, weight;
2041     GetPID(curTrack, pid, weight);
2042
2043     if(pid!=0){
2044       fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
2045     }
2046
2047     if(pid == 0){
2048       fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
2049     }
2050
2051
2052
2053
2054
2055
2056     TLorentzVector curElec;
2057     curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
2058
2059
2060     if(fDoMCTruth){
2061       Int_t labelMC = TMath::Abs(curTrack->GetLabel());
2062       TParticle* curParticle = fStack->Particle(labelMC);
2063       if(curTrack->GetSign() > 0){
2064         if( pid == 0){
2065           fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
2066           fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
2067         }
2068         else{
2069           fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
2070           fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
2071         }
2072       }
2073     }
2074
2075
2076     if(curTrack->GetSign() > 0){
2077
2078       //     vESDxPosTemp.push_back(curTrack);
2079       new((*vESDxPosTemp)[vESDxPosTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
2080
2081       if( pid == 0){
2082
2083         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
2084         fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
2085         //      fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
2086         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
2087         //      fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
2088         //      vESDePosTemp.push_back(curTrack);
2089         new((*vESDePosTemp)[vESDePosTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
2090
2091       }
2092
2093     }
2094     else {
2095
2096       new((*vESDxNegTemp)[vESDxNegTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
2097
2098       if( pid == 0){
2099
2100         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
2101         fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
2102         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
2103         new((*vESDeNegTemp)[vESDeNegTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
2104
2105       }
2106
2107     }
2108
2109   }
2110
2111
2112   Bool_t ePosJPsi = kFALSE;
2113   Bool_t eNegJPsi = kFALSE;
2114   Bool_t ePosPi0  = kFALSE;
2115   Bool_t eNegPi0  = kFALSE;
2116
2117   UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
2118
2119   for(Int_t iNeg=0; iNeg < vESDeNegTemp->GetEntriesFast(); iNeg++){
2120     if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetPdgCode() == 11)
2121       if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0) > -1){
2122         Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0);
2123         TParticle* partMother = fStack ->Particle(labelMother);
2124         if (partMother->GetPdgCode() == 111){
2125           ieNegPi0 = iNeg;
2126           eNegPi0 = kTRUE;
2127         }
2128         if(partMother->GetPdgCode() == 443){ //Mother JPsi
2129           fHistograms->FillTable("Table_Electrons",14);
2130           ieNegJPsi = iNeg;
2131           eNegJPsi = kTRUE;
2132         }
2133         else{
2134           //      vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
2135           new((*vESDeNegNoJPsi)[vESDeNegNoJPsi->GetEntriesFast()])  AliESDtrack(*(AliESDtrack*)(vESDeNegTemp->At(iNeg)));
2136           //            cout<<"ESD No Positivo JPsi "<<endl;
2137         }
2138
2139       }
2140   }
2141
2142   for(Int_t iPos=0; iPos < vESDePosTemp->GetEntriesFast(); iPos++){
2143     if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetPdgCode() == -11)
2144       if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0) > -1){
2145         Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0);
2146         TParticle* partMother = fStack ->Particle(labelMother);
2147         if (partMother->GetPdgCode() == 111){
2148           iePosPi0 = iPos;
2149           ePosPi0 = kTRUE;
2150         }
2151         if(partMother->GetPdgCode() == 443){ //Mother JPsi
2152           fHistograms->FillTable("Table_Electrons",15);
2153           iePosJPsi = iPos;
2154           ePosJPsi = kTRUE;
2155         }
2156         else{
2157           //      vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
2158           new((*vESDePosNoJPsi)[vESDePosNoJPsi->GetEntriesFast()])  AliESDtrack(*(AliESDtrack*)(vESDePosTemp->At(iPos)));
2159           //            cout<<"ESD No Negativo JPsi "<<endl;
2160         }
2161
2162       }
2163   }
2164
2165   if( eNegJPsi && ePosJPsi ){
2166     TVector3 tempeNegV,tempePosV;
2167     tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Pz());
2168     tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Pz());
2169     fHistograms->FillTable("Table_Electrons",16);
2170     fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));
2171     fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->GetLabel())),
2172                                                                               fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->GetLabel()))));
2173   }
2174
2175   if( eNegPi0 && ePosPi0 ){
2176     TVector3 tempeNegV,tempePosV;
2177     tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Pz());
2178     tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Pz());
2179     fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
2180     fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->GetLabel())),
2181                                                                              fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->GetLabel()))));
2182   }
2183
2184
2185   FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
2186
2187   CleanWithAngleCuts(*vESDeNegTemp,*vESDePosTemp,*fKFReconstructedGammasTClone);
2188
2189   //  vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
2190   //  vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
2191
2192   TClonesArray vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectronTClone);
2193   TClonesArray vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectronTClone);
2194
2195
2196   FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
2197
2198
2199
2200
2201   //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
2202
2203
2204   FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
2205   FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
2206
2207
2208
2209   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",*fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
2210
2211   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
2212                            *fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
2213
2214   //BackGround
2215
2216   //Like Sign e+e-
2217   ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
2218   ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
2219   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
2220   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
2221
2222   //        Like Sign e+e- no JPsi
2223   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
2224   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
2225
2226   //Mixed Event
2227
2228   if( fCurrentEventPosElectronTClone->GetEntriesFast() > 0 && fCurrentEventNegElectronTClone->GetEntriesFast() > 0 && fKFReconstructedGammasCutTClone->GetEntriesFast() > 0 ){
2229     FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
2230                              *fKFReconstructedGammasCutTClone,*fPreviousEventTLVNegElectronTClone,*fPreviousEventTLVPosElectronTClone);
2231     *fPreviousEventTLVNegElectronTClone = vCurrentTLVeNeg;
2232     *fPreviousEventTLVPosElectronTClone = vCurrentTLVePos;
2233
2234   }
2235
2236   /*
2237   //Photons P
2238   Double_t vtx[3];
2239   vtx[0]=0;vtx[1]=0;vtx[2]=0;
2240   for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
2241
2242   //      if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
2243
2244
2245
2246   Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
2247   //      cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
2248
2249   //      cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
2250
2251   if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
2252
2253   fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
2254
2255
2256   }
2257
2258
2259   */
2260
2261
2262   vESDeNegTemp->Delete();
2263   vESDePosTemp->Delete();
2264   vESDxNegTemp->Delete();
2265   vESDxPosTemp->Delete();
2266   vESDeNegNoJPsi->Delete();
2267   vESDePosNoJPsi->Delete();
2268
2269   delete vESDeNegTemp;
2270   delete vESDePosTemp;
2271   delete vESDxNegTemp;
2272   delete vESDxPosTemp;
2273   delete vESDeNegNoJPsi;
2274   delete vESDePosNoJPsi;
2275 }
2276
2277 /*
2278   void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
2279   //see header file for documentation
2280   for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
2281   for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
2282   fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
2283   }
2284   }
2285   }
2286 */
2287 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,TClonesArray const tlVeNeg, TClonesArray const tlVePos){
2288   //see header file for documentation
2289   for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++){
2290     for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
2291       fHistograms->FillHistogram(histoName.Data(),((TLorentzVector*)(tlVeNeg.At(iNeg)))->Vect().Angle(((TLorentzVector*)(tlVePos.At(iPos)))->Vect()));
2292     }
2293   }
2294 }
2295 void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, TClonesArray const eNeg, TClonesArray const ePos){
2296   //see header file for documentation
2297   for( Int_t n=0; n < eNeg.GetEntriesFast(); n++){
2298     TLorentzVector en = (*(TLorentzVector*)(eNeg.At(n)));
2299     for (Int_t p=0; p < ePos.GetEntriesFast(); p++){
2300       TLorentzVector ep = (*(TLorentzVector*)(ePos.At(p)));
2301       TLorentzVector np = ep + en;
2302       fHistograms->FillHistogram(histoName.Data(),np.M());
2303     }
2304   }
2305 }
2306
2307 void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,TClonesArray const fKFGammas,
2308                                                               TClonesArray const tlVeNeg,TClonesArray const tlVePos)
2309 {
2310   //see header file for documentation
2311
2312   for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++ ){
2313
2314     for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
2315
2316       TLorentzVector xy = *((TLorentzVector *)(tlVePos.At(iPos))) + *((TLorentzVector *)(tlVeNeg.At(iNeg)));
2317
2318       for (Int_t iGam=0; iGam < fKFGammas.GetEntriesFast(); iGam++){
2319
2320         //      AliKFParticle * gammaCandidate = &fKFGammas[iGam];
2321         AliKFParticle * gammaCandidate = (AliKFParticle *)(fKFGammas.At(iGam));
2322         TLorentzVector g;
2323
2324         g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
2325         TLorentzVector xyg = xy + g;
2326         fHistograms->FillHistogram(histoMass.Data(),xyg.M());
2327         fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
2328       }
2329     }
2330   }
2331
2332 }
2333 void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, TClonesArray e)
2334 {
2335   // see header file for documentation
2336   for(Int_t i=0; i < e.GetEntriesFast(); i++)
2337     {
2338       for (Int_t j=i+1; j < e.GetEntriesFast(); j++)
2339         {
2340           TLorentzVector ee = (*(TLorentzVector*)(e.At(i))) + (*(TLorentzVector*)(e.At(j)));
2341
2342           fHistograms->FillHistogram(hBg.Data(),ee.M());
2343         }
2344     }
2345 }
2346
2347
2348 void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negativeElectrons,
2349                                                         TClonesArray const positiveElectrons,
2350                                                         TClonesArray const gammas){
2351   // see header file for documentation
2352
2353   UInt_t  sizeN = negativeElectrons.GetEntriesFast();
2354   UInt_t  sizeP = positiveElectrons.GetEntriesFast();
2355   UInt_t  sizeG = gammas.GetEntriesFast();
2356
2357
2358
2359   vector <Bool_t> xNegBand(sizeN);
2360   vector <Bool_t> xPosBand(sizeP);
2361   vector <Bool_t> gammaBand(sizeG);
2362
2363
2364   for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
2365   for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
2366   for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
2367
2368
2369   for(UInt_t iPos=0; iPos < sizeP; iPos++){
2370
2371     Double_t aP[3];
2372     ((AliESDtrack*)(positiveElectrons.At(iPos)))->GetConstrainedPxPyPz(aP);
2373
2374     TVector3 ePosV(aP[0],aP[1],aP[2]);
2375
2376     for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
2377
2378       Double_t aN[3];
2379       ((AliESDtrack*)(negativeElectrons.At(iNeg)))->GetConstrainedPxPyPz(aN);
2380       TVector3 eNegV(aN[0],aN[1],aN[2]);
2381
2382       if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
2383         xPosBand[iPos]=kFALSE;
2384         xNegBand[iNeg]=kFALSE;
2385       }
2386
2387       for(UInt_t iGam=0; iGam < sizeG; iGam++){
2388         AliKFParticle* gammaCandidate = (AliKFParticle*)gammas.At(iGam);
2389         TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
2390         if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
2391           gammaBand[iGam]=kFALSE;
2392       }
2393     }
2394   }
2395
2396
2397
2398
2399   for(UInt_t iPos=0; iPos < sizeP; iPos++){
2400     if(xPosBand[iPos]){
2401       new((*fCurrentEventPosElectronTClone)[fCurrentEventPosElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(positiveElectrons.At(iPos))));
2402       //      fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
2403     }
2404   }
2405   for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
2406     if(xNegBand[iNeg]){
2407       new((*fCurrentEventNegElectronTClone)[fCurrentEventNegElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(negativeElectrons.At(iNeg))));
2408       //      fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
2409     }
2410   }
2411   for(UInt_t iGam=0; iGam < sizeG; iGam++){
2412     if(gammaBand[iGam]){
2413       new((*fKFReconstructedGammasCutTClone)[fKFReconstructedGammasCutTClone->GetEntriesFast()]) AliKFParticle((*(AliKFParticle*)(gammas.At(iGam))));
2414       //fKFReconstructedGammasCut.push_back(*(AliKFParticle*)gammas->At(iGam));
2415     }
2416   }
2417 }
2418
2419
2420 void  AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)
2421 {
2422   // see header file for documentation
2423   pid = -1;
2424   weight = -1;
2425
2426   double wpart[5];
2427   double wpartbayes[5];
2428
2429   //get probability of the diffenrent particle types
2430   track->GetESDpid(wpart);
2431
2432   // Tentative particle type "concentrations"
2433   double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
2434
2435   //Bayes' formula
2436   double rcc = 0.;
2437   for (int i = 0; i < 5; i++)
2438     {
2439       rcc+=(c[i] * wpart[i]);
2440     }
2441
2442
2443
2444   for (int i=0; i<5; i++) {
2445     if( rcc>0 || rcc<0){//Kenneth: not sure if the rcc<0 should be there, this is from fixing a coding violation where we are not allowed to say: rcc!=0 (RC19)
2446       wpartbayes[i] = c[i] * wpart[i] / rcc;
2447     }
2448   }
2449
2450
2451
2452   Float_t max=0.;
2453   int ipid=-1;
2454   //find most probable particle in ESD pid
2455   //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
2456   for (int i = 0; i < 5; i++)
2457     {
2458       if (wpartbayes[i] > max)
2459         {
2460           ipid = i;
2461           max = wpartbayes[i];
2462         }
2463     }
2464
2465   pid = ipid;
2466   weight = max;
2467 }
2468 double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)
2469 {
2470   // Calculates the number of sigma to the vertex.
2471
2472   Float_t b[2];
2473   Float_t bRes[2];
2474   Float_t bCov[3];
2475   t->GetImpactParameters(b,bCov);
2476   if (bCov[0]<=0 || bCov[2]<=0) {
2477     AliDebug(1, "Estimated b resolution lower or equal zero!");
2478     bCov[0]=0; bCov[2]=0;
2479   }
2480   bRes[0] = TMath::Sqrt(bCov[0]);
2481   bRes[1] = TMath::Sqrt(bCov[2]);
2482
2483   // -----------------------------------
2484   // How to get to a n-sigma cut?
2485   //
2486   // The accumulated statistics from 0 to d is
2487   //
2488   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
2489   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
2490   //
2491   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
2492   // Can this be expressed in a different way?
2493
2494   if (bRes[0] == 0 || bRes[1] ==0)
2495     return -1;
2496
2497   double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
2498
2499   // stupid rounding problem screws up everything:
2500   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
2501   if (TMath::Exp(-d * d / 2) < 1e-10)
2502     return 1000;
2503
2504
2505   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
2506   return d;
2507 }
2508
2509 //vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
2510 TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *const esdTrack){
2511   //Return TLoresntz vector of track?
2512   //  vector <TLorentzVector> tlVtrack(0);
2513   TClonesArray array("TLorentzVector",0);
2514
2515   for(Int_t itrack=0; itrack < esdTrack->GetEntriesFast(); itrack++){
2516     double p[3];
2517     //esdTrack[itrack]->GetConstrainedPxPyPz(p);
2518     ((AliESDtrack*)(esdTrack->At(itrack)))->GetConstrainedPxPyPz(p);
2519     TLorentzVector currentTrack;
2520     currentTrack.SetXYZM(p[0],p[1],p[2],fElectronMass);
2521     new((array)[array.GetEntriesFast()])  TLorentzVector(currentTrack);
2522     //    tlVtrack.push_back(currentTrack);
2523   }
2524
2525   return array;
2526
2527   //  return tlVtrack;
2528 }
2529