PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx

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