PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx

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