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