PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx

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