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