PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx

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