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