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 "AliESDInputHandler.h"
  34 #include "AliAnalysisManager.h"
  35 #include "AliGammaConversionAODObject.h"
  36 #include "AliGammaConversionBGHandler.h"
  37 #include "AliESDCaloCluster.h" // for combining PHOS and GammaConv
  38 #include "AliKFVertex.h"
  39 class AliESDTrackCuts;
  40 class AliCFContainer;
  41 class AliCFManager;
  42 class AliKFVertex;
  43 class AliAODHandler;
  44 class AliAODEvent;
  45 class ALiESDEvent;
  46 class AliMCEvent;
  47 class AliMCEventHandler;
  48 class AliESDInputHandler;
  49 class AliAnalysisManager;
  50 class Riostream;
  51 class TFile;
  52 class TInterpreter;
  53 class TSystem;
  54 class TROOT;
  55
  56 ClassImp(AliAnalysisTaskGammaConversion)
  57
  58
  59 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
  60 AliAnalysisTaskSE(),
  61   fV0Reader(NULL),
  62   fStack(NULL),
  63   fMCTruth(NULL),    // for CF
  64   fGCMCEvent(NULL),    // for CF
  65   fESDEvent(NULL),
  66   fOutputContainer(NULL),
  67   fCFManager(0x0),   // for CF
  68   fHistograms(NULL),
  69   fTriggerCINT1B(kFALSE),
  70   fDoMCTruth(kFALSE),
  71   fDoNeutralMeson(kFALSE),
  72   fDoOmegaMeson(kFALSE),
  73   fDoJet(kFALSE),
  74   fDoChic(kFALSE),
  75   fRecalculateV0ForGamma(kFALSE),
  76   fKFReconstructedGammasTClone(NULL),
  77   fKFReconstructedPi0sTClone(NULL),
  78   fKFRecalculatedGammasTClone(NULL),
  79   fCurrentEventPosElectronTClone(NULL),
  80   fCurrentEventNegElectronTClone(NULL),
  81   fKFReconstructedGammasCutTClone(NULL),
  82   fPreviousEventTLVNegElectronTClone(NULL),
  83   fPreviousEventTLVPosElectronTClone(NULL),
  84   fElectronv1(),
  85   fElectronv2(),
  86   fGammav1(),
  87   fGammav2(),
  88   fElectronRecalculatedv1(),
  89   fElectronRecalculatedv2(),
  90   fElectronMass(-1),
  91   fGammaMass(-1),
  92   fPi0Mass(-1),
  93   fEtaMass(-1),
  94   fGammaWidth(-1),
  95   fPi0Width(-1),
  96   fEtaWidth(-1),
  97   fMinOpeningAngleGhostCut(0.),
  98   fEsdTrackCuts(NULL),
  99   fCalculateBackground(kFALSE),
 100   fWriteNtuple(kFALSE),
 101   fGammaNtuple(NULL),
 102   fNeutralMesonNtuple(NULL),
 103   fTotalNumberOfAddedNtupleEntries(0),
 104   fChargedParticles(NULL),
 105   fChargedParticlesId(),
 106   fGammaPtHighest(0.),
 107   fMinPtForGammaJet(1.),
 108   fMinIsoConeSize(0.2),
 109   fMinPtIsoCone(0.7),
 110   fMinPtGamChargedCorr(0.5),
 111   fMinPtJetCone(0.5),
 112   fLeadingChargedIndex(-1),
 113   fLowPtMapping(1.),
 114   fHighPtMapping(3.),
 115   fDoCF(kFALSE),
 116   fAODGamma(NULL),
 117   fAODPi0(NULL),
 118   fAODOmega(NULL),
 119   fAODBranchName("GammaConv"),
 120   fKFForceAOD(kFALSE),
 121   fKFDeltaAODFileName(""),
 122   fDoNeutralMesonV0MCCheck(kFALSE),
 123   fKFReconstructedGammasV0Index()
 124 {
 125   // Default constructor
 126
 127   /*   Kenneth: the default constructor should not have any define input/output or the call to SetESDtrackCuts
 128   // Common I/O in slot 0
 129   DefineInput (0, TChain::Class());
 130   DefineOutput(0, TTree::Class());
 131
 132   // Your private output
 133   DefineOutput(1, TList::Class());
 134
 135   // Define standard ESD track cuts for Gamma-hadron correlation
 136   SetESDtrackCuts();
 137   */
 138 }
 139
 140 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
 141   AliAnalysisTaskSE(name),
 142   fV0Reader(NULL),
 143   fStack(NULL),
 144   fMCTruth(NULL),    // for CF
 145   fGCMCEvent(NULL),    // for CF
 146   fESDEvent(NULL),
 147   fOutputContainer(0x0),
 148   fCFManager(0x0),   // for CF
 149   fHistograms(NULL),
 150   fTriggerCINT1B(kFALSE),
 151   fDoMCTruth(kFALSE),
 152   fDoNeutralMeson(kFALSE),
 153   fDoOmegaMeson(kFALSE),
 154   fDoJet(kFALSE),
 155   fDoChic(kFALSE),
 156   fRecalculateV0ForGamma(kFALSE),
 157   fKFReconstructedGammasTClone(NULL),
 158   fKFReconstructedPi0sTClone(NULL),
 159   fKFRecalculatedGammasTClone(NULL),
 160   fCurrentEventPosElectronTClone(NULL),
 161   fCurrentEventNegElectronTClone(NULL),
 162   fKFReconstructedGammasCutTClone(NULL),
 163   fPreviousEventTLVNegElectronTClone(NULL),
 164   fPreviousEventTLVPosElectronTClone(NULL),
 165   fElectronv1(),
 166   fElectronv2(),
 167   fGammav1(),
 168   fGammav2(),
 169   fElectronRecalculatedv1(),
 170   fElectronRecalculatedv2(),
 171   fElectronMass(-1),
 172   fGammaMass(-1),
 173   fPi0Mass(-1),
 174   fEtaMass(-1),
 175   fGammaWidth(-1),
 176   fPi0Width(-1),
 177   fEtaWidth(-1),
 178   fMinOpeningAngleGhostCut(0.),
 179   fEsdTrackCuts(NULL),
 180   fCalculateBackground(kFALSE),
 181   fWriteNtuple(kFALSE),
 182   fGammaNtuple(NULL),
 183   fNeutralMesonNtuple(NULL),
 184   fTotalNumberOfAddedNtupleEntries(0),
 185   fChargedParticles(NULL),
 186   fChargedParticlesId(),
 187   fGammaPtHighest(0.),
 188   fMinPtForGammaJet(1.),
 189   fMinIsoConeSize(0.2),
 190   fMinPtIsoCone(0.7),
 191   fMinPtGamChargedCorr(0.5),
 192   fMinPtJetCone(0.5),
 193   fLeadingChargedIndex(-1),
 194   fLowPtMapping(1.),
 195   fHighPtMapping(3.),
 196    fDoCF(kFALSE),
 197    fAODGamma(NULL),
 198    fAODPi0(NULL),
 199    fAODOmega(NULL),
 200    fAODBranchName("GammaConv"),
 201    fKFForceAOD(kFALSE),
 202    fKFDeltaAODFileName(""),
 203    fDoNeutralMesonV0MCCheck(kFALSE),
 204   fKFReconstructedGammasV0Index()
 205 {
 206   // Common I/O in slot 0
 207   DefineInput (0, TChain::Class());
 208   DefineOutput(0, TTree::Class());
 209
 210   // Your private output
 211   DefineOutput(1, TList::Class());
 212   DefineOutput(2, AliCFContainer::Class());  // for CF
 213
 214
 215   // Define standard ESD track cuts for Gamma-hadron correlation
 216   SetESDtrackCuts();
 217
 218 }
 219
 220 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
 221 {
 222   // Remove all pointers
 223
 224   if(fOutputContainer){
 225     fOutputContainer->Clear() ;
 226     delete fOutputContainer ;
 227   }
 228   if(fHistograms){
 229     delete fHistograms;
 230   }
 231   if(fV0Reader){
 232     delete fV0Reader;
 233   }
 234
 235   // for CF
 236   if(fCFManager){
 237     delete fCFManager;
 238   }
 239
 240   if(fEsdTrackCuts){
 241     delete fEsdTrackCuts;
 242   }
 243
 244   //Delete AODs
 245   if (fAODGamma) {
 246     fAODGamma->Clear();
 247     delete fAODGamma;
 248   }
 249   fAODGamma = NULL;
 250
 251   if (fAODPi0) {
 252     fAODPi0->Clear();
 253     delete fAODPi0;
 254   }
 255   fAODPi0 = NULL;
 256
 257   if (fAODOmega) {
 258     fAODOmega->Clear();
 259     delete fAODOmega;
 260   }
 261   fAODOmega = NULL;
 262
 263
 264 }
 265
 266
 267 void AliAnalysisTaskGammaConversion::Init()
 268 {
 269   // Initialization
 270   // AliLog::SetGlobalLogLevel(AliLog::kError);
 271 }
 272 void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
 273 {
 274   // SetESDtrackCuts
 275   if (fEsdTrackCuts!=NULL){
 276     delete fEsdTrackCuts;
 277   }
 278   fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
 279   //standard cuts from:
 280   //http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
 281
 282   // Cuts used up to 3rd of March
 283
 284  //  fEsdTrackCuts->SetMinNClustersTPC(50);
 285 //   fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 286 //   fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
 287 //   fEsdTrackCuts->SetRequireITSRefit(kTRUE);
 288 //   fEsdTrackCuts->SetMaxNsigmaToVertex(3);
 289 //   fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
 290
 291   //------- To be tested-----------
 292   // Cuts used up to 26th of Agost
 293 //    Int_t minNClustersTPC = 70;
 294 //    Double_t maxChi2PerClusterTPC = 4.0;
 295 //    Double_t maxDCAtoVertexXY = 2.4; // cm
 296 //    Double_t maxDCAtoVertexZ  = 3.2; // cm
 297 //    fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 298 //    fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
 299 //    fEsdTrackCuts->SetRequireITSRefit(kTRUE);
 300 //    //   fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
 301 //    fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 302 //    fEsdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
 303 //    fEsdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
 304 //    fEsdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
 305 //    fEsdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
 306 //    fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
 307 //    fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
 308 //    fEsdTrackCuts->SetPtRange(0.15);
 309
 310 //    fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
 311
 312
 313 // Using standard function  for setting Cuts
 314   Bool_t selectPrimaries=kTRUE;
 315   fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selectPrimaries);
 316   fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
 317   fEsdTrackCuts->SetPtRange(0.15);
 318
 319   //-----  From Jacek 10.03.03 ------------------/
 320 //     minNClustersTPC = 70;
 321 //     maxChi2PerClusterTPC = 4.0;
 322 //     maxDCAtoVertexXY = 2.4; // cm
 323 //     maxDCAtoVertexZ  = 3.2; // cm
 324
 325 //     esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 326 //     esdTrackCuts->SetRequireTPCRefit(kFALSE);
 327 //     esdTrackCuts->SetRequireTPCStandAlone(kTRUE);
 328 //     esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 329 //     esdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
 330 //     esdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
 331 //     esdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
 332 //     esdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
 333 //     esdTrackCuts->SetDCAToVertex2D(kTRUE);
 334
 335
 336
 337   //  fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 338   //  fV0Reader->SetESDtrackCuts(fEsdTrackCuts);
 339 }
 340
 341 void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
 342 {
 343   // Execute analysis for current event
 344
 345   //  Load the esdpid from the esdhandler if exists (tender was applied) otherwise set the Bethe Bloch parameters
 346
 347   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 348   AliESDInputHandler *esdHandler=0x0;
 349   if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
 350     AliV0Reader::SetESDpid(esdHandler->GetESDpid());
 351   } else {
 352     //load esd pid bethe bloch parameters depending on the existance of the MC handler
 353     // yes: MC parameters
 354     // no:  data parameters
 355     if (!AliV0Reader::GetESDpid()){
 356       if (fMCEvent ) {
 357         AliV0Reader::InitESDpid();
 358       } else {
 359         AliV0Reader::InitESDpid(1);
 360       }
 361     }
 362   }
 363
 364   //Must set fForceAOD to true for the AOD to get filled. Should only be done when running independent chain / train.
 365   if(fKFForceAOD) {
 366     if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");
 367     AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
 368   }
 369
 370   if(fV0Reader == NULL){
 371     // Write warning here cuts and so on are default if this ever happens
 372   }
 373
 374   fV0Reader->SetInputAndMCEvent(InputEvent(), MCEvent());
 375
 376   fV0Reader->Initialize();
 377   fDoMCTruth = fV0Reader->GetDoMCTruth();
 378
 379
 380   //  if(fDoMCTruth){
 381   //ConnectInputData("");
 382     //}
 383   //Each event needs an empty branch
 384   //  fAODBranch->Clear();
 385   if(fAODGamma) fAODGamma->Delete();
 386   if(fAODPi0) fAODPi0->Delete();
 387   if(fAODOmega) fAODOmega->Delete();
 388
 389   if(fKFReconstructedGammasTClone == NULL){
 390     fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
 391   }
 392   if(fCurrentEventPosElectronTClone == NULL){
 393     fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
 394   }
 395   if(fCurrentEventNegElectronTClone == NULL){
 396     fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
 397   }
 398   if(fKFReconstructedGammasCutTClone == NULL){
 399     fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
 400   }
 401   if(fPreviousEventTLVNegElectronTClone == NULL){
 402     fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
 403   }
 404   if(fPreviousEventTLVPosElectronTClone == NULL){
 405     fPreviousEventTLVPosElectronTClone  = new TClonesArray("TLorentzVector",0);
 406   }
 407   if(fChargedParticles == NULL){
 408     fChargedParticles = new TClonesArray("AliESDtrack",0);
 409   }
 410
 411   if(fKFReconstructedPi0sTClone == NULL){
 412     fKFReconstructedPi0sTClone = new TClonesArray("AliKFParticle",0);
 413   }
 414
 415  if(fKFRecalculatedGammasTClone == NULL){
 416     fKFRecalculatedGammasTClone = new TClonesArray("AliKFParticle",0);
 417   }
 418
 419
 420   //clear TClones
 421   fKFReconstructedGammasTClone->Delete();
 422   fCurrentEventPosElectronTClone->Delete();
 423   fCurrentEventNegElectronTClone->Delete();
 424   fKFReconstructedGammasCutTClone->Delete();
 425   fPreviousEventTLVNegElectronTClone->Delete();
 426   fPreviousEventTLVPosElectronTClone->Delete();
 427   fKFReconstructedPi0sTClone->Delete();
 428   fKFRecalculatedGammasTClone->Delete();
 429
 430   //clear vectors
 431   fElectronv1.clear();
 432   fElectronv2.clear();
 433
 434   fGammav1.clear();
 435   fGammav2.clear();
 436
 437   fElectronRecalculatedv1.clear();
 438   fElectronRecalculatedv2.clear();
 439
 440
 441   fChargedParticles->Delete();
 442
 443   fChargedParticlesId.clear();
 444
 445   fKFReconstructedGammasV0Index.clear();
 446
 447   //Clear the data in the v0Reader
 448   //  fV0Reader->UpdateEventByEventData();
 449
 450   //Take Only events with proper trigger
 451   /*
 452   if(fTriggerCINT1B){
 453     if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
 454   }
 455   */
 456
 457   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
 458     //    cout<< "Event not taken"<< endl;
 459     return; // aborts if the primary vertex does not have contributors.
 460   }
 461
 462
 463   if(!fV0Reader->CheckForPrimaryVertexZ() ){
 464     return;
 465   }
 466
 467
 468
 469   // Process the MC information
 470   if(fDoMCTruth){
 471     ProcessMCData();
 472   }
 473
 474   //Process the v0 information with no cuts
 475   ProcessV0sNoCut();
 476
 477   // Process the v0 information
 478   ProcessV0s();
 479
 480
 481   //Fill Gamma AOD
 482   FillAODWithConversionGammas() ;
 483
 484   // Process reconstructed gammas
 485   if(fDoNeutralMeson == kTRUE){
 486     ProcessGammasForNeutralMesonAnalysis();
 487
 488   }
 489
 490   if(fDoMCTruth == kTRUE){
 491     CheckV0Efficiency();
 492   }
 493   //Process reconstructed gammas electrons for Chi_c Analysis
 494   if(fDoChic == kTRUE){
 495     ProcessGammaElectronsForChicAnalysis();
 496   }
 497   // Process reconstructed gammas for gamma Jet/hadron correlations
 498   if(fDoJet == kTRUE){
 499     ProcessGammasForGammaJetAnalysis();
 500   }
 501
 502   //calculate background if flag is set
 503   if(fCalculateBackground){
 504     CalculateBackground();
 505   }
 506
 507    if(fDoNeutralMeson == kTRUE){
 508 //     ProcessConvPHOSGammasForNeutralMesonAnalysis();
 509      if(fDoOmegaMeson == kTRUE){
 510        ProcessGammasForOmegaMesonAnalysis();
 511      }
 512    }
 513
 514   //Clear the data in the v0Reader
 515   fV0Reader->UpdateEventByEventData();
 516   if(fRecalculateV0ForGamma==kTRUE){
 517     RecalculateV0ForGamma();
 518   }
 519   PostData(1, fOutputContainer);
 520   PostData(2, fCFManager->GetParticleContainer());  // for CF
 521
 522 }
 523
 524 // void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t *option){
 525 //   // see header file for documentation
 526 //   //  printf("   ConnectInputData %s\n", GetName());
 527
 528 //   AliAnalysisTaskSE::ConnectInputData(option);
 529
 530 //   if(fV0Reader == NULL){
 531 //     // Write warning here cuts and so on are default if this ever happens
 532 //   }
 533 //   fV0Reader->Initialize();
 534 //   fDoMCTruth = fV0Reader->GetDoMCTruth();
 535 // }
 536
 537
 538
 539 void AliAnalysisTaskGammaConversion::ProcessMCData(){
 540   // see header file for documentation
 541   //InputEvent(), MCEvent());
 542   /* TestAnaMarin
 543   fStack = fV0Reader->GetMCStack();
 544   fMCTruth = fV0Reader->GetMCTruth();  // for CF
 545   fGCMCEvent = fV0Reader->GetMCEvent();  // for CF
 546   */
 547   fStack= MCEvent()->Stack();
 548   fGCMCEvent=MCEvent();
 549
 550   // for CF
 551   Double_t containerInput[3];
 552   if(fDoCF){
 553     if(!fGCMCEvent) cout << "NO MC INFO FOUND" << endl;
 554     fCFManager->SetEventInfo(fGCMCEvent);
 555   }
 556   // end for CF
 557
 558
 559   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
 560     return; // aborts if the primary vertex does not have contributors.
 561   }
 562
 563   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
 564     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
 565
 566     if (!particle) {
 567       //print warning here
 568       continue;
 569     }
 570
 571     ///////////////////////Begin Chic Analysis/////////////////////////////
 572
 573     if(particle->GetPdgCode() == 443){//Is JPsi
 574       if(particle->GetNDaughters()==2){
 575         if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
 576            TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
 577
 578           TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
 579           TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
 580           if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
 581             fHistograms->FillTable("Table_Electrons",3);//e+ e-  from J/Psi inside acceptance
 582
 583           if( TMath::Abs(daug0->Eta()) < 0.9){
 584             if(daug0->GetPdgCode() == -11)
 585               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
 586             else
 587               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
 588
 589           }
 590           if(TMath::Abs(daug1->Eta()) < 0.9){
 591             if(daug1->GetPdgCode() == -11)
 592               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
 593             else
 594               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
 595           }
 596         }
 597       }
 598     }
 599     //              const int CHI_C0   = 10441;
 600     //              const int CHI_C1   = 20443;
 601     //              const int CHI_C2   = 445
 602     if(particle->GetPdgCode() == 22){//gamma from JPsi
 603       if(particle->GetMother(0) > -1){
 604         if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
 605            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
 606            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
 607           if(TMath::Abs(particle->Eta()) < 1.2)
 608             fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
 609         }
 610       }
 611     }
 612     if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
 613       if( particle->GetNDaughters() == 2){
 614         TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
 615         TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
 616
 617         if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
 618           if( daug0->GetPdgCode() == 443){
 619             TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
 620             TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
 621             if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
 622               fHistograms->FillTable("Table_Electrons",18);
 623
 624           }//if
 625           else if (daug1->GetPdgCode() == 443){
 626             TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
 627             TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
 628             if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
 629               fHistograms->FillTable("Table_Electrons",18);
 630           }//else if
 631         }//gamma o Jpsi
 632       }//GetNDaughters
 633     }
 634
 635
 636     /////////////////////End Chic Analysis////////////////////////////
 637
 638
 639     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )    continue;
 640
 641     if(particle->R()>fV0Reader->GetMaxRCut())   continue; // cuts on distance from collision point
 642
 643     Double_t tmpPhi=particle->Phi();
 644
 645     if(particle->Phi()> TMath::Pi()){
 646       tmpPhi = particle->Phi()-(2*TMath::Pi());
 647     }
 648
 649     Double_t rapidity;
 650     if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
 651       rapidity=0;
 652     }
 653     else{
 654       rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
 655     }
 656
 657     //process the gammas
 658     if (particle->GetPdgCode() == 22){
 659
 660
 661       if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
 662         continue; // no photon as mothers!
 663       }
 664
 665       if(particle->GetMother(0) >= fStack->GetNprimary()){
 666         continue; // the gamma has a mother, and it is not a primary particle
 667       }
 668
 669       fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
 670       fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
 671       fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
 672       fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
 673       fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
 674
 675       // for CF
 676       if(fDoCF){
 677         containerInput[0] = particle->Pt();
 678         containerInput[1] = particle->Eta();
 679         if(particle->GetMother(0) >=0){
 680           containerInput[2] = fStack->Particle(particle->GetMother(0))->GetMass();
 681         }
 682         else{
 683           containerInput[2]=-1;
 684         }
 685
 686         fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);                                        // generated gamma
 687       }
 688       if(particle->GetMother(0) < 0){   // direct gamma
 689         fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
 690         fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
 691         fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
 692         fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
 693         fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);
 694       }
 695
 696       // looking for conversion (electron + positron from pairbuilding (= 5) )
 697       TParticle* ePos = NULL;
 698       TParticle* eNeg = NULL;
 699
 700       if(particle->GetNDaughters() >= 2){
 701         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
 702           TParticle *tmpDaughter = fStack->Particle(daughterIndex);
 703           if(tmpDaughter->GetUniqueID() == 5){
 704             if(tmpDaughter->GetPdgCode() == 11){
 705               eNeg = tmpDaughter;
 706             }
 707             else if(tmpDaughter->GetPdgCode() == -11){
 708               ePos = tmpDaughter;
 709             }
 710           }
 711         }
 712       }
 713
 714
 715       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
 716         continue;
 717       }
 718
 719
 720       Double_t ePosPhi = ePos->Phi();
 721       if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
 722
 723       Double_t eNegPhi = eNeg->Phi();
 724       if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
 725
 726
 727       if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
 728         continue; // no reconstruction below the Pt cut
 729       }
 730
 731       if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
 732         continue;
 733       }
 734
 735       if(ePos->R()>fV0Reader->GetMaxRCut()){
 736         continue; // cuts on distance from collision point
 737       }
 738
 739       if(TMath::Abs(ePos->Vz()) > fV0Reader->GetMaxZCut()){
 740         continue;   // outside material
 741       }
 742
 743
 744       if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  > ePos->R()){
 745         continue;               // line cut to exclude regions where we do not reconstruct
 746       }
 747
 748
 749       // for CF
 750       if(fDoCF){
 751         fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructable);  // reconstructable gamma
 752       }
 753       fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
 754       fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
 755       fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
 756       fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
 757       fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
 758       fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
 759
 760       fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
 761       fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
 762       fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
 763       fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
 764
 765       fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
 766       fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
 767       fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
 768       fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
 769
 770
 771       // begin Mapping
 772       Int_t rBin    = fHistograms->GetRBin(ePos->R());
 773       Int_t zBin    = fHistograms->GetZBin(ePos->Vz());
 774       Int_t phiBin  = fHistograms->GetPhiBin(particle->Phi());
 775       Double_t rFMD=30;
 776
 777       TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());
 778
 779       TString nameMCMappingPhiR="";
 780       nameMCMappingPhiR.Form("MC_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
 781       // fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
 782
 783       TString nameMCMappingPhi="";
 784       nameMCMappingPhi.Form("MC_Conversion_Mapping_Phi%02d",phiBin);
 785       //      fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
 786       //fHistograms->FillHistogram(nameMCMappingPhi, ePos->Vz(), particle->Eta());
 787
 788       TString nameMCMappingR="";
 789       nameMCMappingR.Form("MC_Conversion_Mapping_R%02d",rBin);
 790       //      fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
 791       //fHistograms->FillHistogram(nameMCMappingR,ePos->Vz(), particle->Eta());
 792
 793       TString nameMCMappingPhiInR="";
 794       nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_in_R_%02d",rBin);
 795       //      fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
 796       fHistograms->FillHistogram(nameMCMappingPhiInR, vtxPos.Phi());
 797
 798       TString nameMCMappingZInR="";
 799       nameMCMappingZInR.Form("MC_Conversion_Mapping_Z_in_R_%02d",rBin);
 800       fHistograms->FillHistogram(nameMCMappingZInR,ePos->Vz() );
 801
 802
 803       TString nameMCMappingPhiInZ="";
 804       nameMCMappingPhiInZ.Form("MC_Conversion_Mapping_Phi_in_Z_%02d",zBin);
 805       //      fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
 806       fHistograms->FillHistogram(nameMCMappingPhiInZ, vtxPos.Phi());
 807
 808
 809       if(ePos->R()<rFMD){
 810         TString nameMCMappingFMDPhiInZ="";
 811         nameMCMappingFMDPhiInZ.Form("MC_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
 812         fHistograms->FillHistogram(nameMCMappingFMDPhiInZ, vtxPos.Phi());
 813       }
 814
 815       TString nameMCMappingRInZ="";
 816       nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
 817       fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
 818
 819       if(particle->Pt() > fLowPtMapping && particle->Pt()< fHighPtMapping){
 820         TString nameMCMappingMidPtPhiInR="";
 821         nameMCMappingMidPtPhiInR.Form("MC_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
 822         fHistograms->FillHistogram(nameMCMappingMidPtPhiInR, vtxPos.Phi());
 823
 824         TString nameMCMappingMidPtZInR="";
 825         nameMCMappingMidPtZInR.Form("MC_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
 826         fHistograms->FillHistogram(nameMCMappingMidPtZInR,ePos->Vz() );
 827
 828
 829         TString nameMCMappingMidPtPhiInZ="";
 830         nameMCMappingMidPtPhiInZ.Form("MC_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
 831         fHistograms->FillHistogram(nameMCMappingMidPtPhiInZ, vtxPos.Phi());
 832
 833
 834         if(ePos->R()<rFMD){
 835           TString nameMCMappingMidPtFMDPhiInZ="";
 836           nameMCMappingMidPtFMDPhiInZ.Form("MC_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
 837           fHistograms->FillHistogram(nameMCMappingMidPtFMDPhiInZ, vtxPos.Phi());
 838         }
 839
 840
 841         TString nameMCMappingMidPtRInZ="";
 842         nameMCMappingMidPtRInZ.Form("MC_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
 843         fHistograms->FillHistogram(nameMCMappingMidPtRInZ,ePos->R() );
 844
 845       }
 846
 847       //end mapping
 848
 849       fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
 850       fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
 851       fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
 852       fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
 853       fHistograms->FillHistogram("MC_ConvGamma_E_AsymmetryP",particle->P(),eNeg->P()/particle->P());
 854       fHistograms->FillHistogram("MC_ConvGamma_P_AsymmetryP",particle->P(),ePos->P()/particle->P());
 855
 856
 857       if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
 858         fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
 859         fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
 860         fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
 861         fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
 862         fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
 863
 864       } // end direct gamma
 865       else{   // mother exits
 866         /*      if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0
 867                 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
 868                 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2
 869                 ){
 870                 fMCGammaChic.push_back(particle);
 871                 }
 872         */
 873       }  // end if mother exits
 874     } // end if particle is a photon
 875
 876
 877
 878     // process motherparticles (2 gammas as daughters)
 879     // the motherparticle had already to pass the R and the eta cut, but no line cut.
 880     // the line cut is just valid for the conversions!
 881
 882     if(particle->GetNDaughters() == 2){
 883
 884       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
 885       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
 886
 887       if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
 888
 889       // Check the acceptance for both gammas
 890       Bool_t gammaEtaCut = kTRUE;
 891       if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut()  ) gammaEtaCut = kFALSE;
 892       Bool_t gammaRCut = kTRUE;
 893       if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut()  ) gammaRCut = kFALSE;
 894
 895
 896
 897       // check for conversions now -> have to pass eta, R and line cut!
 898       Bool_t daughter0Electron = kFALSE;
 899       Bool_t daughter0Positron = kFALSE;
 900       Bool_t daughter1Electron = kFALSE;
 901       Bool_t daughter1Positron = kFALSE;
 902
 903       if(daughter0->GetNDaughters() >= 2){  // first gamma
 904         for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
 905           TParticle *tmpDaughter = fStack->Particle(TrackIndex);
 906           if(tmpDaughter->GetUniqueID() == 5){
 907             if(tmpDaughter->GetPdgCode() == 11){
 908               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
 909                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
 910                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
 911                     daughter0Electron = kTRUE;
 912                   }
 913                 }
 914
 915               }
 916             }
 917             else if(tmpDaughter->GetPdgCode() == -11){
 918               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
 919                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
 920                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
 921                     daughter0Positron = kTRUE;
 922                   }
 923                 }
 924               }
 925             }
 926           }
 927         }
 928       }
 929
 930
 931       if(daughter1->GetNDaughters() >= 2){   // second gamma
 932         for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
 933           TParticle *tmpDaughter = fStack->Particle(TrackIndex);
 934           if(tmpDaughter->GetUniqueID() == 5){
 935             if(tmpDaughter->GetPdgCode() == 11){
 936               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
 937                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
 938                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
 939                     daughter1Electron = kTRUE;
 940                   }
 941                 }
 942
 943               }
 944             }
 945             else if(tmpDaughter->GetPdgCode() == -11){
 946               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
 947                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){
 948                   if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
 949                     daughter1Positron = kTRUE;
 950                   }
 951                 }
 952
 953               }
 954
 955             }
 956           }
 957         }
 958       }
 959
 960
 961       if(particle->GetPdgCode()==111){     //Pi0
 962         if( iTracks >= fStack->GetNprimary()){
 963           fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
 964           fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
 965           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
 966           fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
 967           fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
 968           fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
 969           fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
 970           fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
 971           fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
 972
 973           if(gammaEtaCut && gammaRCut){
 974             //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
 975             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
 976             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
 977             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
 978               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
 979               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
 980             }
 981           }
 982         }
 983         else{
 984           fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());
 985           fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
 986           fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
 987           fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
 988           fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
 989           fHistograms->FillHistogram("MC_Pi0_R", particle->R());
 990           fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
 991           fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
 992           fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
 993           if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_Fiducial", particle->Pt());
 994
 995           if(gammaEtaCut && gammaRCut){
 996             //    if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
 997             fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
 998             fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
 999             if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_withinAcceptance_Fiducial", particle->Pt());
1000
1001             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1002               fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1003               fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1004               fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
1005               fHistograms->FillHistogram("MC_Pi0_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1006               fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1007               fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1008
1009               Double_t alfa=0.;
1010               if((daughter0->Energy()+daughter1->Energy())!= 0.){
1011                 alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
1012               }
1013               fHistograms->FillHistogram("MC_Pi0_alpha",alfa);
1014               if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_ConvGamma_withinAcceptance_Fiducial", particle->Pt());
1015
1016             }
1017           }
1018         }
1019       }
1020
1021       if(particle->GetPdgCode()==221){   //Eta
1022         fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
1023         fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
1024         fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
1025         fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
1026         fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
1027         fHistograms->FillHistogram("MC_Eta_R", particle->R());
1028         fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
1029         fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1030         fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
1031
1032         if(gammaEtaCut && gammaRCut){
1033           //    if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1034           fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1035           fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1036           if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1037             fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1038             fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1039             fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
1040             fHistograms->FillHistogram("MC_Eta_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1041             fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1042             fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1043
1044           }
1045
1046         }
1047
1048       }
1049
1050       // all motherparticles with 2 gammas as daughters
1051       fHistograms->FillHistogram("MC_Mother_R", particle->R());
1052       fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
1053       fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
1054       fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
1055       fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1056       fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
1057       fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
1058       fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
1059       fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
1060       fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
1061       fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());
1062
1063       if(gammaEtaCut && gammaRCut){
1064         //      if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1065         fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1066         fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1067         fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());
1068         if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1069           fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1070           fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1071           fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());
1072
1073         }
1074
1075
1076       }  // end passed R and eta cut
1077
1078     } // end if(particle->GetNDaughters() == 2)
1079
1080   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
1081
1082 } // end ProcessMCData
1083
1084
1085
1086 void AliAnalysisTaskGammaConversion::FillNtuple(){
1087   //Fills the ntuple with the different values
1088
1089   if(fGammaNtuple == NULL){
1090     return;
1091   }
1092   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1093   for(Int_t i=0;i<numberOfV0s;i++){
1094     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};
1095     AliESDv0 * cV0 = fV0Reader->GetV0(i);
1096     Double_t negPID=0;
1097     Double_t posPID=0;
1098     fV0Reader->GetPIDProbability(negPID,posPID);
1099     values[0]=cV0->GetOnFlyStatus();
1100     values[1]=fV0Reader->CheckForPrimaryVertex();
1101     values[2]=negPID;
1102     values[3]=posPID;
1103     values[4]=fV0Reader->GetX();
1104     values[5]=fV0Reader->GetY();
1105     values[6]=fV0Reader->GetZ();
1106     values[7]=fV0Reader->GetXYRadius();
1107     values[8]=fV0Reader->GetMotherCandidateNDF();
1108     values[9]=fV0Reader->GetMotherCandidateChi2();
1109     values[10]=fV0Reader->GetMotherCandidateEnergy();
1110     values[11]=fV0Reader->GetMotherCandidateEta();
1111     values[12]=fV0Reader->GetMotherCandidatePt();
1112     values[13]=fV0Reader->GetMotherCandidateMass();
1113     values[14]=fV0Reader->GetMotherCandidateWidth();
1114     //    values[15]=fV0Reader->GetMotherMCParticle()->Pt();   MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
1115     values[16]=fV0Reader->GetOpeningAngle();
1116     values[17]=fV0Reader->GetNegativeTrackEnergy();
1117     values[18]=fV0Reader->GetNegativeTrackPt();
1118     values[19]=fV0Reader->GetNegativeTrackEta();
1119     values[20]=fV0Reader->GetNegativeTrackPhi();
1120     values[21]=fV0Reader->GetPositiveTrackEnergy();
1121     values[22]=fV0Reader->GetPositiveTrackPt();
1122     values[23]=fV0Reader->GetPositiveTrackEta();
1123     values[24]=fV0Reader->GetPositiveTrackPhi();
1124     values[25]=fV0Reader->HasSameMCMother();
1125     if(values[25] != 0){
1126       values[26]=fV0Reader->GetMotherMCParticlePDGCode();
1127       values[15]=fV0Reader->GetMotherMCParticle()->Pt();
1128     }
1129     fTotalNumberOfAddedNtupleEntries++;
1130     fGammaNtuple->Fill(values);
1131   }
1132   fV0Reader->ResetV0IndexNumber();
1133
1134 }
1135
1136 void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
1137   // Process all the V0's without applying any cuts to it
1138
1139   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1140   for(Int_t i=0;i<numberOfV0s;i++){
1141     /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
1142
1143     if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
1144       continue;
1145     }
1146
1147     //    if( !fV0Reader->GetV0(i)->GetOnFlyStatus()){
1148     if( !fV0Reader->CheckV0FinderStatus(i)){
1149       continue;
1150     }
1151
1152
1153     if( !((fV0Reader->GetNegativeESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ||
1154         !((fV0Reader->GetPositiveESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ){
1155       continue;
1156     }
1157
1158     if( fV0Reader->GetNegativeESDTrack()->GetSign()== fV0Reader->GetPositiveESDTrack()->GetSign()){
1159       continue;
1160     }
1161
1162     if( fV0Reader->GetNegativeESDTrack()->GetKinkIndex(0) > 0 ||
1163         fV0Reader->GetPositiveESDTrack()->GetKinkIndex(0) > 0) {
1164       continue;
1165     }
1166     if(TMath::Abs(fV0Reader->GetMotherCandidateEta())> fV0Reader->GetEtaCut()){
1167       continue;
1168     }
1169     if(TMath::Abs(fV0Reader->GetPositiveTrackEta())> fV0Reader->GetEtaCut()){
1170       continue;
1171     }
1172     if(TMath::Abs(fV0Reader->GetNegativeTrackEta())> fV0Reader->GetEtaCut()){
1173       continue;
1174     }
1175     if((TMath::Abs(fV0Reader->GetZ())*fV0Reader->GetLineCutZRSlope())-fV0Reader->GetLineCutZValue() > fV0Reader->GetXYRadius() ){ // cuts out regions where we do not reconstruct
1176       continue;
1177     }
1178     if(fDoMCTruth){
1179
1180       if(fV0Reader->HasSameMCMother() == kFALSE){
1181         continue;
1182       }
1183
1184       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1185       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1186
1187       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1188         continue;
1189       }
1190       if(negativeMC->GetPdgCode() == positiveMC->GetPdgCode()){
1191         continue;
1192       }
1193
1194       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
1195         continue;
1196       }
1197
1198       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1199
1200         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1201         fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1202         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1203         fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1204         fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1205         fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1206         fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1207         fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1208         fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1209         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1210
1211         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1212         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1213
1214         fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1215         fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
1216         fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1217         fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1218         fHistograms->FillHistogram("ESD_NoCutConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1219         fHistograms->FillHistogram("ESD_NoCutConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1220         fHistograms->FillHistogram("ESD_NoCutConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1221         fHistograms->FillHistogram("ESD_NoCutConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1222
1223         fHistograms->FillHistogram("ESD_NoCutConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1224         fHistograms->FillHistogram("ESD_NoCutConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1225         fHistograms->FillHistogram("ESD_NoCutConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1226         fHistograms->FillHistogram("ESD_NoCutConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1227
1228         //store MCTruth properties
1229         fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1230         fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1231         fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1232       }
1233     }
1234   }
1235   fV0Reader->ResetV0IndexNumber();
1236 }
1237
1238 void AliAnalysisTaskGammaConversion::ProcessV0s(){
1239   // see header file for documentation
1240
1241
1242   if(fWriteNtuple == kTRUE){
1243     FillNtuple();
1244   }
1245
1246   Int_t nSurvivingV0s=0;
1247   while(fV0Reader->NextV0()){
1248     nSurvivingV0s++;
1249
1250
1251     //-------------------------- filling v0 information -------------------------------------
1252     fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
1253     fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1254     fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1255     fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1256
1257     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
1258     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
1259     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
1260     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
1261     fHistograms->FillHistogram("ESD_E_nTPCClusters", fV0Reader->GetNegativeTracknTPCClusters());
1262     fHistograms->FillHistogram("ESD_E_nITSClusters", fV0Reader->GetNegativeTracknITSClusters());
1263
1264     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
1265     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
1266     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
1267     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
1268     fHistograms->FillHistogram("ESD_P_nTPCClusters", fV0Reader->GetPositiveTracknTPCClusters());
1269     fHistograms->FillHistogram("ESD_P_nITSClusters", fV0Reader->GetPositiveTracknITSClusters());
1270
1271     fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1272     fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1273     fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1274     fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1275     fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1276     fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1277     fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1278     fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1279     fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1280     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1281
1282     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1283     fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1284
1285     fHistograms->FillHistogram("ESD_ConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1286     fHistograms->FillHistogram("ESD_ConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1287     fHistograms->FillHistogram("ESD_ConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1288     fHistograms->FillHistogram("ESD_ConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1289
1290     fHistograms->FillHistogram("ESD_ConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1291     fHistograms->FillHistogram("ESD_ConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1292     fHistograms->FillHistogram("ESD_ConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1293     fHistograms->FillHistogram("ESD_ConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1294
1295     Double_t negPID=0;
1296     Double_t posPID=0;
1297     fV0Reader->GetPIDProbability(negPID,posPID);
1298     fHistograms->FillHistogram("ESD_ConvGamma_E_EProbP",fV0Reader->GetNegativeTrackP(),negPID);
1299     fHistograms->FillHistogram("ESD_ConvGamma_P_EProbP",fV0Reader->GetPositiveTrackP(),posPID);
1300
1301     Double_t negPIDmupi=0;
1302     Double_t posPIDmupi=0;
1303     fV0Reader->GetPIDProbabilityMuonPion(negPIDmupi,posPIDmupi);
1304     fHistograms->FillHistogram("ESD_ConvGamma_E_mupiProbP",fV0Reader->GetNegativeTrackP(),negPIDmupi);
1305     fHistograms->FillHistogram("ESD_ConvGamma_P_mupiProbP",fV0Reader->GetPositiveTrackP(),posPIDmupi);
1306
1307     Double_t armenterosQtAlfa[2];
1308     fV0Reader->GetArmenterosQtAlfa(fV0Reader-> GetNegativeKFParticle(),
1309                                    fV0Reader-> GetPositiveKFParticle(),
1310                                    fV0Reader->GetMotherCandidateKFCombination(),
1311                                    armenterosQtAlfa);
1312
1313     fHistograms->FillHistogram("ESD_ConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1314
1315
1316     // begin mapping
1317     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());
1318     Int_t zBin    = fHistograms->GetZBin(fV0Reader->GetZ());
1319     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
1320     Double_t rFMD=30;
1321
1322     TVector3 vtxConv(fV0Reader->GetX(),fV0Reader->GetY(), fV0Reader->GetZ());
1323
1324     //    Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
1325
1326     TString nameESDMappingPhiR="";
1327     nameESDMappingPhiR.Form("ESD_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
1328     //fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
1329
1330     TString nameESDMappingPhi="";
1331     nameESDMappingPhi.Form("ESD_Conversion_Mapping_Phi%02d",phiBin);
1332     //fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
1333
1334     TString nameESDMappingR="";
1335     nameESDMappingR.Form("ESD_Conversion_Mapping_R%02d",rBin);
1336     //fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);
1337
1338     TString nameESDMappingPhiInR="";
1339     nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_in_R_%02d",rBin);
1340     //    fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1341     fHistograms->FillHistogram(nameESDMappingPhiInR, vtxConv.Phi());
1342
1343     TString nameESDMappingZInR="";
1344     nameESDMappingZInR.Form("ESD_Conversion_Mapping_Z_in_R_%02d",rBin);
1345     fHistograms->FillHistogram(nameESDMappingZInR, fV0Reader->GetZ());
1346
1347     TString nameESDMappingPhiInZ="";
1348     nameESDMappingPhiInZ.Form("ESD_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1349     //    fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1350     fHistograms->FillHistogram(nameESDMappingPhiInZ, vtxConv.Phi());
1351
1352     if(fV0Reader->GetXYRadius()<rFMD){
1353       TString nameESDMappingFMDPhiInZ="";
1354       nameESDMappingFMDPhiInZ.Form("ESD_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
1355       fHistograms->FillHistogram(nameESDMappingFMDPhiInZ, vtxConv.Phi());
1356     }
1357
1358
1359     TString nameESDMappingRInZ="";
1360     nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
1361     fHistograms->FillHistogram(nameESDMappingRInZ, fV0Reader->GetXYRadius());
1362
1363     if(fV0Reader->GetMotherCandidatePt() > fLowPtMapping && fV0Reader->GetMotherCandidatePt()< fHighPtMapping){
1364       TString nameESDMappingMidPtPhiInR="";
1365       nameESDMappingMidPtPhiInR.Form("ESD_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1366       fHistograms->FillHistogram(nameESDMappingMidPtPhiInR, vtxConv.Phi());
1367
1368       TString nameESDMappingMidPtZInR="";
1369       nameESDMappingMidPtZInR.Form("ESD_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1370       fHistograms->FillHistogram(nameESDMappingMidPtZInR, fV0Reader->GetZ());
1371
1372       TString nameESDMappingMidPtPhiInZ="";
1373       nameESDMappingMidPtPhiInZ.Form("ESD_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1374       fHistograms->FillHistogram(nameESDMappingMidPtPhiInZ, vtxConv.Phi());
1375       if(fV0Reader->GetXYRadius()<rFMD){
1376         TString nameESDMappingMidPtFMDPhiInZ="";
1377         nameESDMappingMidPtFMDPhiInZ.Form("ESD_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
1378         fHistograms->FillHistogram(nameESDMappingMidPtFMDPhiInZ, vtxConv.Phi());
1379       }
1380
1381
1382       TString nameESDMappingMidPtRInZ="";
1383       nameESDMappingMidPtRInZ.Form("ESD_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1384       fHistograms->FillHistogram(nameESDMappingMidPtRInZ, fV0Reader->GetXYRadius());
1385
1386     }
1387
1388
1389     // end mapping
1390
1391     new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()])  AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
1392     fKFReconstructedGammasV0Index.push_back(fV0Reader->GetCurrentV0IndexNumber()-1);
1393     //    fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
1394     fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
1395     fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
1396
1397     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
1398     if(fDoMCTruth){
1399
1400       if(fV0Reader->HasSameMCMother() == kFALSE){
1401         continue;
1402       }
1403       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1404       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1405
1406       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1407         continue;
1408       }
1409
1410       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1411         continue;
1412       }
1413       if(negativeMC->GetUniqueID() == 4 && positiveMC->GetUniqueID() ==4){// fill r distribution for Dalitz decays
1414         if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 111){ //pi0
1415           fHistograms->FillHistogram("ESD_TrueDalitzContamination_R", fV0Reader->GetXYRadius());
1416         }
1417       }
1418
1419       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
1420         continue;
1421       }
1422
1423       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1424
1425         if(fDoCF){
1426           Double_t containerInput[3];
1427           containerInput[0] = fV0Reader->GetMotherCandidatePt();
1428           containerInput[1] = fV0Reader->GetMotherCandidateEta();
1429           containerInput[2] = fV0Reader->GetMotherCandidateMass();
1430           fCFManager->GetParticleContainer()->Fill(containerInput,kStepTrueGamma); // for CF
1431         }
1432
1433         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1434         fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1435         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1436         fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1437         fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1438         fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1439         fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1440         fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1441         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1442         fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1443         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
1444         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
1445         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1446         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1447
1448         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1449         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1450
1451
1452         fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1453         fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
1454         fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1455         fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1456
1457         fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1458         fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1459         fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1460         fHistograms->FillHistogram("ESD_TrueConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1461         if (fV0Reader->GetMotherCandidateP() != 0) {
1462           fHistograms->FillHistogram("ESD_TrueConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1463           fHistograms->FillHistogram("ESD_TrueConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1464         } else { cout << "Error::fV0Reader->GetNegativeTrackP() == 0 !!!" << endl; }
1465
1466         fHistograms->FillHistogram("ESD_TrueConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1467         fHistograms->FillHistogram("ESD_TrueConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1468         fHistograms->FillHistogram("ESD_TrueConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1469
1470
1471
1472         //store MCTruth properties
1473         fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1474         fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1475         fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1476
1477         //resolution
1478         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();
1479         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();
1480         Double_t resdPt = 0.;
1481         if(mcpt > 0){
1482           resdPt = ((esdpt - mcpt)/mcpt)*100.;
1483         }
1484         else if(mcpt < 0){
1485           cout<<"Pt of MC particle is negative, this will cause wrong calculation of resPt"<<endl;
1486         }
1487
1488         fHistograms->FillHistogram("Resolution_Gamma_dPt_Pt", mcpt, resdPt);
1489         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
1490         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
1491         fHistograms->FillHistogram("Resolution_Gamma_dPt_Phi", fV0Reader->GetMotherCandidatePhi(), resdPt);
1492
1493         Double_t resdZ = 0.;
1494         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
1495           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100.;
1496         }
1497         Double_t resdZAbs = 0.;
1498         resdZAbs = (fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz());
1499
1500         fHistograms->FillHistogram("Resolution_dZAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdZAbs);
1501         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
1502         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
1503         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
1504
1505         // new for dPt_Pt-histograms for Electron and Positron
1506         Double_t mcEpt = fV0Reader->GetNegativeMCParticle()->Pt();
1507         Double_t resEdPt = 0.;
1508         if (mcEpt > 0){
1509                 resEdPt = ((fV0Reader->GetNegativeTrackPt()-mcEpt)/mcEpt)*100.;
1510         }
1511         UInt_t statusN = fV0Reader->GetNegativeESDTrack()->GetStatus();
1512  //    AliESDtrack * negTrk = fV0Reader->GetNegativeESDTrack();
1513         UInt_t kTRDoutN =  (statusN & AliESDtrack::kTRDout);
1514
1515         Int_t ITSclsE= fV0Reader->GetNegativeTracknITSClusters();
1516         // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
1517          switch(ITSclsE){
1518           case 0: // 0 ITS clusters
1519                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS0", mcEpt, resEdPt);
1520             break;
1521           case 1:  // 1 ITS cluster
1522                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS1", mcEpt, resEdPt);
1523             break;
1524           case 2:  // 2 ITS clusters
1525                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS2", mcEpt, resEdPt);
1526             break;
1527           case 3: // 3 ITS clusters
1528                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS3", mcEpt, resEdPt);
1529             break;
1530           case 4: // 4 ITS clusters
1531                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS4", mcEpt, resEdPt);
1532             break;
1533           case 5: // 5 ITS clusters
1534                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS5", mcEpt, resEdPt);
1535             break;
1536           case 6: // 6 ITS clusters
1537                 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS6", mcEpt, resEdPt);
1538             break;
1539           }
1540         //Filling histograms with respect to Electron resolution
1541         fHistograms->FillHistogram("Resolution_E_dPt_Pt", mcEpt, resEdPt);
1542         fHistograms->FillHistogram("Resolution_E_dPt_Phi", fV0Reader->GetNegativeTrackPhi(), resEdPt);
1543         if(kTRDoutN){
1544         fHistograms->FillHistogram("Resolution_E_nTRDtracklets_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
1545         fHistograms->FillHistogram("Resolution_E_nTRDtracklets_MCPt", mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
1546         fHistograms->FillHistogram("Resolution_E_nTRDclusters_ESDPt",fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1547         fHistograms->FillHistogram("Resolution_E_nTRDclusters_MCPt",mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1548         fHistograms->FillHistogram("Resolution_E_TRDsignal_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDsignal());
1549         }
1550
1551         Double_t mcPpt = fV0Reader->GetPositiveMCParticle()->Pt();
1552         Double_t resPdPt = 0;
1553         if (mcPpt > 0){
1554                 resPdPt = ((fV0Reader->GetPositiveTrackPt()-mcPpt)/mcPpt)*100.;
1555         }
1556
1557         UInt_t statusP = fV0Reader->GetPositiveESDTrack()->GetStatus();
1558 //     AliESDtrack * posTr= fV0Reader->GetPositiveESDTrack();
1559         UInt_t kTRDoutP =  (statusP & AliESDtrack::kTRDout);
1560
1561         Int_t ITSclsP = fV0Reader->GetPositiveTracknITSClusters();
1562         // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
1563          switch(ITSclsP){
1564           case 0: // 0 ITS clusters
1565                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS0", mcPpt, resPdPt);
1566             break;
1567           case 1:  // 1 ITS cluster
1568                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS1", mcPpt, resPdPt);
1569             break;
1570           case 2:  // 2 ITS clusters
1571                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS2", mcPpt, resPdPt);
1572             break;
1573           case 3: // 3 ITS clusters
1574                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS3", mcPpt, resPdPt);
1575             break;
1576           case 4: // 4 ITS clusters
1577                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS4", mcPpt, resPdPt);
1578             break;
1579           case 5: // 5 ITS clusters
1580                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS5", mcPpt, resPdPt);
1581             break;
1582           case 6: // 6 ITS clusters
1583                 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS6", mcPpt, resPdPt);
1584             break;
1585           }
1586         //Filling histograms with respect to Positron resolution
1587         fHistograms->FillHistogram("Resolution_P_dPt_Pt", mcPpt, resPdPt);
1588         fHistograms->FillHistogram("Resolution_P_dPt_Phi", fV0Reader->GetPositiveTrackPhi(), resPdPt);
1589         if(kTRDoutP){
1590         fHistograms->FillHistogram("Resolution_P_nTRDtracklets_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1591         fHistograms->FillHistogram("Resolution_P_nTRDtracklets_MCPt", mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1592         fHistograms->FillHistogram("Resolution_P_nTRDclusters_ESDPt",fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1593         fHistograms->FillHistogram("Resolution_P_nTRDclusters_MCPt",mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1594         fHistograms->FillHistogram("Resolution_P_TRDsignal_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDsignal());
1595         }
1596
1597
1598         Double_t resdR = 0.;
1599         if(fV0Reader->GetNegativeMCParticle()->R() != 0){
1600           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100.;
1601         }
1602         Double_t resdRAbs = 0.;
1603         resdRAbs = (fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R());
1604
1605         fHistograms->FillHistogram("Resolution_dRAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdRAbs);
1606         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
1607         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
1608         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
1609         fHistograms->FillHistogram("Resolution_R_dPt", fV0Reader->GetNegativeMCParticle()->R(), resdPt);
1610
1611         Double_t resdPhiAbs=0.;
1612         resdPhiAbs=0.;
1613         resdPhiAbs= (fV0Reader->GetMotherCandidatePhi()-fV0Reader->GetNegativeMCParticle()->Phi());
1614         fHistograms->FillHistogram("Resolution_dPhiAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdPhiAbs);
1615
1616       }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
1617     }//if(fDoMCTruth)
1618   }//while(fV0Reader->NextV0)
1619   fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
1620   fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
1621   fHistograms->FillHistogram("ESD_NumberOfContributorsVtx", fV0Reader->GetNumberOfContributorsVtx());
1622
1623   fV0Reader->ResetV0IndexNumber();
1624 }
1625
1626 void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
1627   // Fill AOD with reconstructed Gamma
1628
1629   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
1630     //  for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
1631     //Create AOD particle object from AliKFParticle
1632     //You could add directly AliKFParticle objects to the AOD, avoiding dependences with PartCorr
1633     //but this means that I have to work a little bit more in my side.
1634     //AODPWG4Particle objects are simpler and lighter, I think
1635     /*
1636     AliKFParticle * gammakf = dynamic_cast<AliKFParticle*>(fKFReconstructedGammasTClone->At(gammaIndex));
1637     AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
1638     //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
1639     gamma.SetTrackLabel( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
1640     gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
1641     gamma.SetPdg(AliPID::kEleCon); //photon id
1642     gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
1643     gamma.SetChi2(gammakf->Chi2());
1644     Int_t i = fAODBranch->GetEntriesFast();
1645     new((*fAODBranch)[i])  AliAODPWG4Particle(gamma);
1646     */
1647
1648     AliKFParticle * gammakf = (AliKFParticle *)fKFReconstructedGammasTClone->At(gammaIndex);
1649     AliGammaConversionAODObject aodObject;
1650     aodObject.SetPx(gammakf->GetPx());
1651     aodObject.SetPy(gammakf->GetPy());
1652     aodObject.SetPz(gammakf->GetPz());
1653     aodObject.SetLabel1(fElectronv1[gammaIndex]);
1654     aodObject.SetLabel2(fElectronv2[gammaIndex]);
1655     aodObject.SetChi2(gammakf->Chi2());
1656     aodObject.SetE(gammakf->E());
1657     Int_t i = fAODGamma->GetEntriesFast();
1658     new((*fAODGamma)[i])  AliGammaConversionAODObject(aodObject);
1659   }
1660
1661 }
1662
1663 void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
1664   // omega meson analysis pi0+gamma decay
1665   for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
1666     AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
1667     for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1668
1669       AliKFParticle * omegaCandidateGammaDaughter = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1670       if(fGammav1[firstPi0Index]==firstGammaIndex || fGammav2[firstPi0Index]==firstGammaIndex){
1671         continue;
1672       }
1673
1674       AliKFParticle omegaCandidate(*omegaCandidatePi0Daughter,*omegaCandidateGammaDaughter);
1675       Double_t massOmegaCandidate = 0.;
1676       Double_t widthOmegaCandidate = 0.;
1677
1678       omegaCandidate.GetMass(massOmegaCandidate,widthOmegaCandidate);
1679
1680       AddOmegaToAOD(&omegaCandidate, massOmegaCandidate, firstPi0Index, firstGammaIndex);
1681
1682       fHistograms->FillHistogram("ESD_Omega_InvMass_vs_Pt",massOmegaCandidate ,omegaCandidate.GetPt());
1683       fHistograms->FillHistogram("ESD_Omega_InvMass",massOmegaCandidate);
1684
1685       //delete omegaCandidate;
1686
1687     }// end of omega reconstruction in pi0+gamma channel
1688
1689     if(fDoJet == kTRUE){
1690       AliKFParticle* negPiKF=NULL;
1691       AliKFParticle* posPiKF=NULL;
1692
1693       // look at the pi+pi+pi0 channel
1694       for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
1695         AliESDtrack* posTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
1696         if (posTrack->GetSign()<0) continue;
1697         if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion))>2.) continue;
1698         if (posPiKF) delete posPiKF; posPiKF=NULL;
1699         posPiKF = new AliKFParticle( *(posTrack) ,211);
1700
1701         for(Int_t jCh=0;jCh<fChargedParticles->GetEntriesFast();jCh++){
1702           AliESDtrack* negTrack = (AliESDtrack*)(fChargedParticles->At(jCh));
1703           if( negTrack->GetSign()>0) continue;
1704           if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion))>2.) continue;
1705           if (negPiKF) delete negPiKF; negPiKF=NULL;
1706           negPiKF = new AliKFParticle( *(negTrack) ,-211);
1707           AliKFParticle omegaCandidatePipPinPi0(*omegaCandidatePi0Daughter,*posPiKF,*negPiKF);
1708           Double_t massOmegaCandidatePipPinPi0 = 0.;
1709           Double_t widthOmegaCandidatePipPinPi0 = 0.;
1710
1711           omegaCandidatePipPinPi0.GetMass(massOmegaCandidatePipPinPi0,widthOmegaCandidatePipPinPi0);
1712
1713           fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass_vs_Pt",massOmegaCandidatePipPinPi0 ,omegaCandidatePipPinPi0.GetPt());
1714           fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass",massOmegaCandidatePipPinPi0);
1715
1716           //  delete omegaCandidatePipPinPi0;
1717         }
1718       }
1719     } // checking ig gammajet because in that case the chargedparticle list is created
1720
1721
1722
1723   }
1724
1725   if(fCalculateBackground){
1726     // Background calculation for the omega
1727     for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
1728       AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
1729       for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1730         AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
1731         for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
1732           AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
1733           AliKFParticle * omegaBckCandidate = new AliKFParticle(*omegaCandidatePi0Daughter,previousGoodV0);
1734           Double_t massOmegaBckCandidate = 0.;
1735           Double_t widthOmegaBckCandidate = 0.;
1736
1737           omegaBckCandidate->GetMass(massOmegaBckCandidate,widthOmegaBckCandidate);
1738           fHistograms->FillHistogram("ESD_Omega_Bck_InvMass_vs_Pt",massOmegaBckCandidate ,omegaBckCandidate->GetPt());
1739           fHistograms->FillHistogram("ESD_Omega_Bck_InvMass",massOmegaBckCandidate);
1740
1741           delete omegaBckCandidate;
1742         }
1743       }
1744     }
1745   } // end of checking if background calculation is available
1746 }
1747
1748
1749 void AliAnalysisTaskGammaConversion::AddOmegaToAOD(AliKFParticle * omegakf, Double_t mass, Int_t omegaDaughter, Int_t gammaDaughter) {
1750   //See header file for documentation
1751   AliGammaConversionAODObject omega;
1752   omega.SetPx(omegakf->Px());
1753   omega.SetPy(omegakf->Py());
1754   omega.SetPz(omegakf->Pz());
1755   omega.SetChi2(omegakf->GetChi2());
1756   omega.SetE(omegakf->E());
1757   omega.SetIMass(mass);
1758   omega.SetLabel1(omegaDaughter);
1759   //dynamic_cast<AliGammaConversionAODObject*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
1760   omega.SetLabel2(gammaDaughter);
1761   new((*fAODOmega)[fAODOmega->GetEntriesFast()])  AliGammaConversionAODObject(omega);
1762 }
1763
1764
1765
1766 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
1767   // see header file for documentation
1768
1769   //  for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
1770   //    for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
1771
1772   if(fKFReconstructedGammasTClone->GetEntriesFast()>fV0Reader->GetNumberOfV0s()){
1773     cout<<"Warning, number of entries in the tclone is bigger than number of v0s"<<endl;
1774   }
1775
1776   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1777     for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
1778
1779       //      AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
1780       //      AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
1781
1782       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1783       AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
1784
1785       if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1786         continue;
1787       }
1788       if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1789         continue;
1790       }
1791
1792       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
1793
1794       Double_t massTwoGammaCandidate = 0.;
1795       Double_t widthTwoGammaCandidate = 0.;
1796       Double_t chi2TwoGammaCandidate =10000.;
1797       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
1798       //      if(twoGammaCandidate->GetNDF()>0){
1799       //        chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
1800       chi2TwoGammaCandidate = twoGammaCandidate->GetChi2();
1801
1802         fHistograms->FillHistogram("ESD_Mother_Chi2",chi2TwoGammaCandidate);
1803         //if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){
1804
1805           TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
1806           TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
1807
1808           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
1809           Double_t rapidity;
1810           if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() <= 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() <= 0){
1811             cout << "Error: |Pz| > E !!!! " << endl;
1812             rapidity=0;
1813           }
1814           else{
1815             rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
1816           }
1817           Double_t alfa=0.0;
1818           if( (twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()) != 0){
1819             alfa=TMath::Abs((twoGammaDecayCandidateDaughter0->GetE()-twoGammaDecayCandidateDaughter1->GetE())
1820               /(twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()));
1821           }
1822
1823           if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut){
1824             delete twoGammaCandidate;
1825             continue;   // minimum opening angle to avoid using ghosttracks
1826           }
1827
1828           fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
1829           fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
1830           fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
1831           fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
1832           fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);
1833           fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
1834           fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
1835           fHistograms->FillHistogram("ESD_Mother_alfa", alfa);
1836           fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!
1837           fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
1838           fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
1839           fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1840           if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
1841             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1842           }
1843           fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
1844
1845           /* Kenneth, just for testing*/
1846           AliGammaConversionBGHandler * bgHandlerTest = fV0Reader->GetBGHandler();
1847
1848           Int_t zbin= bgHandlerTest->GetZBinIndex(fV0Reader->GetVertexZ());
1849           Int_t mbintest= bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1850
1851           fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass",zbin,mbintest),massTwoGammaCandidate);
1852
1853           /* end Kenneth, just for testing*/
1854
1855           if(fCalculateBackground){
1856             AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
1857             Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1858             fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1859           }
1860           //      if(fDoNeutralMesonV0MCCheck){
1861           if(fDoMCTruth){
1862             //Kenneth: Checking the eta of the gamma to check the difference between 0.9 and 1.2
1863             Int_t indexKF1 = fKFReconstructedGammasV0Index.at(firstGammaIndex);
1864             if(indexKF1<fV0Reader->GetNumberOfV0s()){
1865               fV0Reader->GetV0(indexKF1);//updates to the correct v0
1866               Double_t eta1 = fV0Reader->GetMotherCandidateEta();
1867               Bool_t isRealPi0=kFALSE;
1868               Bool_t isRealEta=kFALSE;
1869               Int_t gamma1MotherLabel=-1;
1870               if(fV0Reader->HasSameMCMother() == kTRUE){
1871                 //cout<<"This v0 is a real v0!!!!"<<endl;
1872                 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1873                 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1874                 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1875                   if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
1876                     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1877                       gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
1878                     }
1879                   }
1880                 }
1881               }
1882               Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
1883               if(indexKF1 == indexKF2){
1884                 cout<<"index of the two KF particles are the same.... should not happen"<<endl;
1885               }
1886               if(indexKF2<fV0Reader->GetNumberOfV0s()){
1887                 fV0Reader->GetV0(indexKF2);
1888                 Double_t eta2 = fV0Reader->GetMotherCandidateEta();
1889                 Int_t gamma2MotherLabel=-1;
1890                 if(fV0Reader->HasSameMCMother() == kTRUE){
1891                   TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1892                   TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1893                   if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1894                     if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
1895                       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1896                         gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
1897                       }
1898                     }
1899                   }
1900                 }
1901                 if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
1902                   if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
1903                     isRealPi0=kTRUE;
1904                   }
1905                   if(fV0Reader->CheckIfEtaIsMother(gamma1MotherLabel)){
1906                     isRealEta=kTRUE;
1907                   }
1908
1909                 }
1910
1911                 if(TMath::Abs(eta1)>0.9 && TMath::Abs(eta2)>0.9){
1912                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
1913                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1914                   if(isRealPi0 || isRealEta){
1915                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_1212",massTwoGammaCandidate);
1916                     fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_1212",openingAngleTwoGammaCandidate);
1917                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1918                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1919                     fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
1920                     if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
1921                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1922                     }
1923                   }
1924                 }
1925                 else if(TMath::Abs(eta1)>0.9 || TMath::Abs(eta2)>0.9){
1926                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
1927                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1928                   if(isRealPi0 || isRealEta){
1929                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_0912",massTwoGammaCandidate);
1930                     fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0912",openingAngleTwoGammaCandidate);
1931                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1932                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1933                     fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
1934                     if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
1935                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1936                     }
1937                   }
1938                 }
1939                 else{
1940                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
1941                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1942                   if(isRealPi0 || isRealEta){
1943                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_0909",massTwoGammaCandidate);
1944                     fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0909",openingAngleTwoGammaCandidate);
1945                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1946                     fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1947                     fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
1948                     if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
1949                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1950                     }
1951                   }
1952                 }
1953               }
1954             }
1955           }
1956           if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 &&  TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
1957             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1958             fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
1959           }
1960
1961           if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
1962             fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
1963             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1964           }
1965           else if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 || TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
1966             fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
1967             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1968           }
1969           else{
1970             fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
1971             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1972           }
1973           Double_t lowMassPi0=0.1;
1974           Double_t highMassPi0=0.15;
1975           if (massTwoGammaCandidate > lowMassPi0 && massTwoGammaCandidate < highMassPi0 ){
1976             new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()])  AliKFParticle(*twoGammaCandidate);
1977             fGammav1.push_back(firstGammaIndex);
1978             fGammav2.push_back(secondGammaIndex);
1979         AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
1980           }
1981
1982           //        }
1983           //}
1984           delete twoGammaCandidate;
1985     }
1986   }
1987 }
1988
1989 void AliAnalysisTaskGammaConversion::AddPionToAOD(AliKFParticle * pionkf, Double_t mass, Int_t daughter1, Int_t daughter2) {
1990   //See header file for documentation
1991   AliGammaConversionAODObject pion;
1992   pion.SetPx(pionkf->Px());
1993   pion.SetPy(pionkf->Py());
1994   pion.SetPz(pionkf->Pz());
1995   pion.SetChi2(pionkf->GetChi2());
1996   pion.SetE(pionkf->E());
1997   pion.SetIMass(mass);
1998   pion.SetLabel1(daughter1);
1999   //dynamic_cast<AliGammaConversionAODObject*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
2000   pion.SetLabel2(daughter2);
2001   new((*fAODPi0)[fAODPi0->GetEntriesFast()])  AliGammaConversionAODObject(pion);
2002
2003 }
2004
2005
2006
2007 void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
2008 /*
2009   // see header file for documentation
2010   // Analyse Pi0 with one photon from Phos and 1 photon from conversions
2011
2012
2013
2014   Double_t vtx[3];
2015   vtx[0] = fV0Reader->GetPrimaryVertex()->GetX();
2016   vtx[1] = fV0Reader->GetPrimaryVertex()->GetY();
2017   vtx[2] = fV0Reader->GetPrimaryVertex()->GetZ();
2018
2019
2020   // Loop over all CaloClusters and consider only the PHOS ones:
2021   AliESDCaloCluster *clu;
2022   TLorentzVector pPHOS;
2023   TLorentzVector gammaPHOS;
2024   TLorentzVector gammaGammaConv;
2025   TLorentzVector pi0GammaConvPHOS;
2026   TLorentzVector gammaGammaConvBck;
2027   TLorentzVector pi0GammaConvPHOSBck;
2028
2029
2030   for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2031     clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2032     if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
2033     clu ->GetMomentum(pPHOS ,vtx);
2034     for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2035       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2036       gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
2037       gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
2038       pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
2039       fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
2040       fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
2041
2042       TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
2043       TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
2044       Double_t opanConvPHOS= v3D0.Angle(v3D1);
2045       if ( opanConvPHOS < 0.35){
2046         fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
2047       }else{
2048         fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
2049       }
2050
2051     }
2052
2053     //  Now the LorentVector pPHOS is obtained and can be paired with the converted proton
2054   }
2055   //==== End of the PHOS cluster selection ============
2056   TLorentzVector pEMCAL;
2057   TLorentzVector gammaEMCAL;
2058   TLorentzVector pi0GammaConvEMCAL;
2059   TLorentzVector pi0GammaConvEMCALBck;
2060
2061    for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2062     clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2063     if ( !clu->IsEMCAL()  || clu->E()<0.1 ) continue;
2064     if (clu->GetNCells() <= 1) continue;
2065     if ( clu->GetTOF()*1e9 < 550  || clu->GetTOF()*1e9 > 750) continue;
2066
2067     clu ->GetMomentum(pEMCAL ,vtx);
2068     for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2069       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2070       gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
2071                              twoGammaDecayCandidateDaughter0->Py(),
2072                              twoGammaDecayCandidateDaughter0->Pz(),0.);
2073       gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
2074       pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
2075       fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
2076       fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
2077       TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
2078                     twoGammaDecayCandidateDaughter0->Py(),
2079                     twoGammaDecayCandidateDaughter0->Pz());
2080       TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
2081
2082
2083       Double_t opanConvEMCAL= v3D0.Angle(v3D1);
2084       if ( opanConvEMCAL < 0.35){
2085         fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
2086       }else{
2087         fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
2088       }
2089
2090     }
2091     if(fCalculateBackground){
2092       for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2093         AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2094         for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2095           AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2096           gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
2097                                     previousGoodV0.Py(),
2098                                     previousGoodV0.Pz(),0.);
2099           pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
2100           fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
2101           fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
2102                                      pi0GammaConvEMCALBck.Pt());
2103         }
2104       }
2105
2106       //  Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
2107     } // end of checking if background photons are available
2108    }
2109   //==== End of the PHOS cluster selection ============
2110 */
2111 }
2112
2113 void AliAnalysisTaskGammaConversion::CalculateBackground(){
2114   // see header file for documentation
2115
2116
2117   TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
2118
2119   for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2120     AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2121     for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2122       AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2123       for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2124         AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2125
2126         AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2127
2128         Double_t massBG =0.;
2129         Double_t widthBG = 0.;
2130         Double_t chi2BG =10000.;
2131         backgroundCandidate->GetMass(massBG,widthBG);
2132         //      if(backgroundCandidate->GetNDF()>0){
2133         //        chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2134           chi2BG = backgroundCandidate->GetChi2();
2135           //      if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){
2136
2137             TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2138             TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2139
2140             Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2141
2142             Double_t rapidity;
2143
2144             if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() <= 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() <= 0){
2145                     cout << "Error: |Pz| > E !!!! " << endl;
2146                     rapidity=0;
2147             } else {
2148               rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2149             }
2150
2151             Double_t alfa=0.0;
2152             if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2153               alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2154                               /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2155             }
2156
2157
2158             if(openingAngleBG < fMinOpeningAngleGhostCut ){
2159               delete backgroundCandidate;
2160               continue;   // minimum opening angle to avoid using ghosttracks
2161             }
2162
2163             // original
2164             fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2165             fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2166             fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
2167             fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2168             fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2169             fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2170             fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2171             fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2172             fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2173             fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2174             fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2175             fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2176
2177             if(alfa>fV0Reader->GetAlphaMinCutMeson() &&   alfa<fV0Reader->GetAlphaCutMeson()){
2178               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2179             }
2180             if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2181               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2182               fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2183             }
2184
2185
2186             // test
2187             AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2188
2189             Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2190             Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2191
2192             fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2193             fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2194             fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
2195             fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2196             fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2197             fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2198             fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2199             fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2200             fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2201             fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2202             fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2203             fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2204
2205             if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2206               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2207               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2208             }
2209             //    }
2210             //  }
2211         delete backgroundCandidate;
2212       }
2213     }
2214   }
2215 }
2216
2217
2218 void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
2219   //ProcessGammasForGammaJetAnalysis
2220
2221   Double_t distIsoMin;
2222
2223   CreateListOfChargedParticles();
2224
2225
2226   //  for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
2227   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
2228     AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
2229     TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
2230
2231     if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
2232       distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
2233
2234       if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
2235         CalculateJetCone(gammaIndex);
2236       }
2237     }
2238   }
2239 }
2240
2241 //____________________________________________________________________
2242 Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(AliESDtrack *const track)
2243 {
2244 //
2245 // check whether particle has good DCAr(Pt) impact
2246 // parameter. Only for TPC+ITS tracks (7*sigma cut)
2247 // Origin: Andrea Dainese
2248 //
2249
2250 Float_t d0z0[2],covd0z0[3];
2251 track->GetImpactParameters(d0z0,covd0z0);
2252 Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
2253 Float_t d0max = 7.*sigma;
2254 if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
2255
2256 return kFALSE;
2257 }
2258
2259
2260 void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
2261   // CreateListOfChargedParticles
2262
2263   fESDEvent = fV0Reader->GetESDEvent();
2264   Int_t numberOfESDTracks=0;
2265   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2266     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
2267
2268     if(!curTrack){
2269       continue;
2270     }
2271     // Not needed if Standard function used.
2272 //     if(!IsGoodImpPar(curTrack)){
2273 //       continue;
2274 //     }
2275
2276     if(fEsdTrackCuts->AcceptTrack(curTrack) ){
2277       new((*fChargedParticles)[fChargedParticles->GetEntriesFast()])  AliESDtrack(*curTrack);
2278       //      fChargedParticles.push_back(curTrack);
2279       fChargedParticlesId.push_back(iTracks);
2280       numberOfESDTracks++;
2281     }
2282   }
2283   fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
2284
2285   if (fV0Reader->GetNumberOfContributorsVtx()>=1){
2286     fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
2287   }
2288 }
2289 void AliAnalysisTaskGammaConversion::RecalculateV0ForGamma(){
2290
2291  Double_t massE=0.00051099892;
2292  TLorentzVector curElecPos;
2293  TLorentzVector curElecNeg;
2294  TLorentzVector curGamma;
2295
2296  TLorentzVector curGammaAt;
2297  TLorentzVector curElecPosAt;
2298  TLorentzVector curElecNegAt;
2299  AliKFVertex primVtxGamma(*(fESDEvent->GetPrimaryVertex()));
2300  AliKFVertex primVtxImprovedGamma = primVtxGamma;
2301
2302  const AliESDVertex *vtxT3D=fESDEvent->GetPrimaryVertex();
2303
2304  Double_t xPrimaryVertex=vtxT3D->GetXv();
2305  Double_t yPrimaryVertex=vtxT3D->GetYv();
2306  Double_t zPrimaryVertex=vtxT3D->GetZv();
2307  // Float_t primvertex[3]={xPrimaryVertex,yPrimaryVertex,zPrimaryVertex};
2308
2309  Float_t nsigmaTPCtrackPos;
2310  Float_t nsigmaTPCtrackNeg;
2311  Float_t nsigmaTPCtrackPosToPion;
2312  Float_t nsigmaTPCtrackNegToPion;
2313  AliKFParticle* negKF=NULL;
2314  AliKFParticle* posKF=NULL;
2315
2316  for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2317    AliESDtrack* posTrack = fESDEvent->GetTrack(iTracks);
2318    if(!posTrack){
2319      continue;
2320    }
2321    if (posKF) delete posKF; posKF=NULL;
2322    if(posTrack->GetSign()<0) continue;
2323    if(!(posTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
2324    if(posTrack->GetKinkIndex(0)>0 ) continue;
2325    if(posTrack->GetNcls(1)<50)continue;
2326    Double_t momPos[3];
2327    //    posTrack->GetConstrainedPxPyPz(momPos);
2328    posTrack->GetPxPyPz(momPos);
2329    AliESDtrack *ptrk=fESDEvent->GetTrack(iTracks);
2330    curElecPos.SetXYZM(momPos[0],momPos[1],momPos[2],massE);
2331    if(TMath::Abs(curElecPos.Eta())<0.9) continue;
2332    posKF = new AliKFParticle( *(posTrack),-11);
2333
2334    nsigmaTPCtrackPos = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kElectron);
2335    nsigmaTPCtrackPosToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion);
2336
2337    if ( nsigmaTPCtrackPos>5.|| nsigmaTPCtrackPos<-2.){
2338      continue;
2339    }
2340
2341    if(pow((momPos[0]*momPos[0]+momPos[1]*momPos[1]+momPos[2]*momPos[2]),0.5)>0.5 && nsigmaTPCtrackPosToPion<1){
2342      continue;
2343    }
2344
2345
2346
2347    for(Int_t jTracks = 0; jTracks < fESDEvent->GetNumberOfTracks(); jTracks++){
2348      AliESDtrack* negTrack = fESDEvent->GetTrack(jTracks);
2349      if(!negTrack){
2350        continue;
2351      }
2352      if (negKF) delete negKF; negKF=NULL;
2353      if(negTrack->GetSign()>0) continue;
2354      if(!(negTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
2355      if(negTrack->GetKinkIndex(0)>0 ) continue;
2356      if(negTrack->GetNcls(1)<50)continue;
2357      Double_t momNeg[3];
2358      //    negTrack->GetConstrainedPxPyPz(momNeg);
2359      negTrack->GetPxPyPz(momNeg);
2360
2361      nsigmaTPCtrackNeg = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kElectron);
2362      nsigmaTPCtrackNegToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion);
2363      if ( nsigmaTPCtrackNeg>5. || nsigmaTPCtrackNeg<-2.){
2364        continue;
2365      }
2366      if(pow((momNeg[0]*momNeg[0]+momNeg[1]*momNeg[1]+momNeg[2]*momNeg[2]),0.5)>0.5 && nsigmaTPCtrackNegToPion<1){
2367        continue;
2368      }
2369      AliESDtrack *ntrk=fESDEvent->GetTrack(jTracks);
2370      curElecNeg.SetXYZM(momNeg[0],momNeg[1],momNeg[2],massE);
2371      if(TMath::Abs(curElecNeg.Eta())<0.9) continue;
2372      negKF = new AliKFParticle( *(negTrack) ,11);
2373
2374      Double_t b=fESDEvent->GetMagneticField();
2375      Double_t xn, xp, dca=ntrk->GetDCA(ptrk,b,xn,xp);
2376      AliExternalTrackParam nt(*ntrk), pt(*ptrk);
2377      nt.PropagateTo(xn,b); pt.PropagateTo(xp,b);
2378
2379
2380      //--- Like in ITSV0Finder
2381      AliExternalTrackParam ntAt0(*ntrk), ptAt0(*ptrk);
2382      Double_t xxP,yyP,alphaP;
2383      Double_t rP[3];
2384
2385      //     if (!ptAt0.GetGlobalXYZat(ptAt0->GetX(),xxP,yyP,zzP)) continue;
2386      if (!ptAt0.GetXYZAt(ptAt0.GetX(),b,rP)) continue;
2387      xxP=rP[0];
2388      yyP=rP[1];
2389      alphaP = TMath::ATan2(yyP,xxP);
2390
2391
2392      ptAt0.Propagate(alphaP,0,b);
2393      Float_t ptfacP  = (1.+100.*TMath::Abs(ptAt0.GetC(b)));
2394
2395      //     Double_t distP      = ptAt0.GetY();
2396      Double_t normP      = ptfacP*TMath::Sqrt(ptAt0.GetSigmaY2());
2397      Double_t normdist0P = TMath::Abs(ptAt0.GetY()/normP);
2398      Double_t normdist1P = TMath::Abs((ptAt0.GetZ()-zPrimaryVertex)/(ptfacP*TMath::Sqrt(ptAt0.GetSigmaZ2())));
2399      Double_t normdistP  = TMath::Sqrt(normdist0P*normdist0P+normdist1P*normdist1P);
2400
2401
2402      Double_t xxN,yyN,alphaN;
2403      Double_t rN[3];
2404      //     if (!ntAt0.GetGlobalXYZat(ntAt0->GetX(),xxN,yyN,zzN)) continue;
2405      if (!ntAt0.GetXYZAt(ntAt0.GetX(),b,rN)) continue;
2406      xxN=rN[0];
2407      yyN=rN[1];
2408
2409      alphaN = TMath::ATan2(yyN,xxN);
2410
2411      ntAt0.Propagate(alphaN,0,b);
2412
2413      Float_t ptfacN  = (1.+100.*TMath::Abs(ntAt0.GetC(b)));
2414      //     Double_t distN      = ntAt0.GetY();
2415      Double_t normN      = ptfacN*TMath::Sqrt(ntAt0.GetSigmaY2());
2416      Double_t normdist0N = TMath::Abs(ntAt0.GetY()/normN);
2417      Double_t normdist1N = TMath::Abs((ntAt0.GetZ()-zPrimaryVertex)/(ptfacN*TMath::Sqrt(ntAt0.GetSigmaZ2())));
2418      Double_t normdistN  = TMath::Sqrt(normdist0N*normdist0N+normdist1N*normdist1N);
2419
2420      //-----------------------------
2421
2422      Double_t momNegAt[3];
2423      nt.GetPxPyPz(momNegAt);
2424      curElecNegAt.SetXYZM(momNegAt[0],momNegAt[1],momNegAt[2],massE);
2425
2426      Double_t momPosAt[3];
2427      pt.GetPxPyPz(momPosAt);
2428      curElecPosAt.SetXYZM(momPosAt[0],momPosAt[1],momPosAt[2],massE);
2429      if(dca>1){
2430        continue;
2431      }
2432
2433      //     Double_t dneg= negTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
2434      //     Double_t dpos= posTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
2435      AliESDv0 vertex(nt,jTracks,pt,iTracks);
2436
2437
2438      Float_t cpa=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
2439
2440
2441
2442      //  cout<< "v0Rr::"<< v0Rr<<endl;
2443      // if (pvertex.GetRr()<0.5){
2444      // continue;
2445      //}
2446      //     cout<<"vertex.GetChi2V0()"<<vertex.GetChi2V0()<<endl;
2447      if(cpa<0.9)continue;
2448      //     if (vertex.GetChi2V0() > 30) continue;
2449      //     cout<<"xp+xn::"<<xp<<" "<<xn<<endl;
2450      if ((xn+xp) < 0.4) continue;
2451      if (TMath::Abs(ntrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05)
2452        if (TMath::Abs(ptrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05) continue;
2453
2454      //cout<<"pass"<<endl;
2455
2456      AliKFParticle v0GammaC;
2457      v0GammaC+=(*negKF);
2458      v0GammaC+=(*posKF);
2459      v0GammaC.SetMassConstraint(0,0.001);
2460      primVtxImprovedGamma+=v0GammaC;
2461      v0GammaC.SetProductionVertex(primVtxImprovedGamma);
2462
2463
2464      curGamma=curElecNeg+curElecPos;
2465      curGammaAt=curElecNegAt+curElecPosAt;
2466
2467      // invariant mass versus pt of K0short
2468
2469      Double_t chi2V0GammaC=100000.;
2470      if( v0GammaC.GetNDF() != 0) {
2471        chi2V0GammaC = v0GammaC.GetChi2()/v0GammaC.GetNDF();
2472      }else{
2473        cout<< "ERROR::v0K0C.GetNDF()" << endl;
2474      }
2475
2476      if(chi2V0GammaC<200 &&chi2V0GammaC>0 ){
2477        if(fHistograms != NULL){
2478          fHistograms->FillHistogram("ESD_RecalculateV0_InvMass",v0GammaC.GetMass());
2479          fHistograms->FillHistogram("ESD_RecalculateV0_Pt",v0GammaC.GetPt());
2480          fHistograms->FillHistogram("ESD_RecalculateV0_E_dEdxP",curElecNegAt.P(),negTrack->GetTPCsignal());
2481          fHistograms->FillHistogram("ESD_RecalculateV0_P_dEdxP",curElecPosAt.P(),posTrack->GetTPCsignal());
2482          fHistograms->FillHistogram("ESD_RecalculateV0_cpa",cpa);
2483          fHistograms->FillHistogram("ESD_RecalculateV0_dca",dca);
2484          fHistograms->FillHistogram("ESD_RecalculateV0_normdistP",normdistP);
2485          fHistograms->FillHistogram("ESD_RecalculateV0_normdistN",normdistN);
2486
2487          new((*fKFRecalculatedGammasTClone)[fKFRecalculatedGammasTClone->GetEntriesFast()])  AliKFParticle(v0GammaC);
2488          fElectronRecalculatedv1.push_back(iTracks);
2489          fElectronRecalculatedv2.push_back(jTracks);
2490        }
2491      }
2492    }
2493  }
2494
2495  for(Int_t firstGammaIndex=0;firstGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();firstGammaIndex++){
2496    for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();secondGammaIndex++){
2497       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(firstGammaIndex);
2498       AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(secondGammaIndex);
2499
2500       if(fElectronRecalculatedv1[firstGammaIndex]==fElectronRecalculatedv1[secondGammaIndex]){
2501         continue;
2502       }
2503       if( fElectronRecalculatedv2[firstGammaIndex]==fElectronRecalculatedv2[secondGammaIndex]){
2504         continue;
2505       }
2506
2507       AliKFParticle twoGammaCandidate(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
2508       fHistograms->FillHistogram("ESD_RecalculateGG_InvMass",twoGammaCandidate.GetMass());
2509       fHistograms->FillHistogram("ESD_RecalculateGG_InvMass_vs_Pt",twoGammaCandidate.GetMass(),twoGammaCandidate.GetPt());
2510    }
2511  }
2512 }
2513 void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
2514   // CaculateJetCone
2515
2516   Double_t cone;
2517   Double_t coneSize=0.3;
2518   Double_t ptJet=0;
2519
2520   //  AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
2521   AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
2522
2523   TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
2524
2525   AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
2526
2527   Double_t momLeadingCharged[3];
2528   leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
2529
2530   TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
2531
2532   Double_t phi1=momentumVectorLeadingCharged.Phi();
2533   Double_t eta1=momentumVectorLeadingCharged.Eta();
2534   Double_t phi3=momentumVectorCurrentGamma.Phi();
2535
2536   for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
2537     AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
2538     Int_t chId = fChargedParticlesId[iCh];
2539     if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
2540     Double_t mom[3];
2541     curTrack->GetConstrainedPxPyPz(mom);
2542     TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
2543     Double_t phi2=momentumVectorChargedParticle.Phi();
2544     Double_t eta2=momentumVectorChargedParticle.Eta();
2545
2546
2547     cone=100.;
2548     if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
2549       cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
2550     }else{
2551       if( (phi2 - phi1)> TMath::TwoPi()-coneSize ){
2552         cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-TMath::TwoPi()-phi1),2) );
2553       }
2554       if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
2555         cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
2556       }
2557     }
2558
2559     if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
2560       ptJet+= momentumVectorChargedParticle.Pt();
2561       Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
2562       Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
2563       fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
2564       fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
2565
2566     }
2567
2568     Double_t dphiHdrGam=phi3-phi2;
2569     if ( dphiHdrGam < (-TMath::PiOver2())){
2570       dphiHdrGam+=(TMath::TwoPi());
2571     }
2572
2573     if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
2574       dphiHdrGam-=(TMath::TwoPi());
2575     }
2576
2577     if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
2578       fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
2579     }
2580   }//track loop
2581
2582
2583 }
2584
2585 Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
2586   // GetMinimumDistanceToCharge
2587
2588   Double_t fIsoMin=100.;
2589   Double_t ptLeadingCharged=-1.;
2590
2591   fLeadingChargedIndex=-1;
2592
2593   AliKFParticle * gammaHighestPt = (AliKFParticle*)fKFReconstructedGammasTClone->At(indexHighestPtGamma);
2594   TVector3 momentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
2595
2596   Double_t phi1=momentumVectorgammaHighestPt.Phi();
2597   Double_t eta1=momentumVectorgammaHighestPt.Eta();
2598
2599   for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
2600     AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
2601     Int_t chId = fChargedParticlesId[iCh];
2602     if(fElectronv1[indexHighestPtGamma]==chId || fElectronv2[indexHighestPtGamma]==chId) continue;
2603     Double_t mom[3];
2604     curTrack->GetConstrainedPxPyPz(mom);
2605     TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
2606     Double_t phi2=momentumVectorChargedParticle.Phi();
2607     Double_t eta2=momentumVectorChargedParticle.Eta();
2608     Double_t iso=pow(  (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
2609
2610     if(momentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
2611       if (iso<fIsoMin){
2612         fIsoMin=iso;
2613       }
2614     }
2615
2616     Double_t dphiHdrGam=phi1-phi2;
2617     if ( dphiHdrGam < (-TMath::PiOver2())){
2618       dphiHdrGam+=(TMath::TwoPi());
2619     }
2620
2621     if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
2622       dphiHdrGam-=(TMath::TwoPi());
2623     }
2624     if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
2625       fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
2626     }
2627
2628     if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
2629       if (momentumVectorChargedParticle.Pt()> ptLeadingCharged && momentumVectorChargedParticle.Pt()>0.1*momentumVectorgammaHighestPt.Pt()){
2630         ptLeadingCharged=momentumVectorChargedParticle.Pt();
2631         fLeadingChargedIndex=iCh;
2632       }
2633     }
2634
2635   }//track loop
2636   fHistograms->FillHistogram("ESD_MinimumIsoDistance",fIsoMin);
2637   return fIsoMin;
2638
2639 }
2640
2641 Int_t  AliAnalysisTaskGammaConversion::GetIndexHighestPtGamma(){
2642   //GetIndexHighestPtGamma
2643
2644   Int_t indexHighestPtGamma=-1;
2645   //Double_t
2646   fGammaPtHighest = -100.;
2647
2648   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2649     AliKFParticle * gammaHighestPtCandidate = (AliKFParticle*)fKFReconstructedGammasTClone->At(firstGammaIndex);
2650     TVector3 momentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
2651     if (momentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
2652       fGammaPtHighest=momentumVectorgammaHighestPtCandidate.Pt();
2653       //gammaHighestPt = gammaHighestPtCandidate;
2654       indexHighestPtGamma=firstGammaIndex;
2655     }
2656   }
2657
2658   return indexHighestPtGamma;
2659
2660 }
2661
2662
2663 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
2664 {
2665   // Terminate analysis
2666   //
2667   AliDebug(1,"Do nothing in Terminate");
2668 }
2669
2670 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
2671 {
2672   //AOD
2673   if(!fAODGamma) fAODGamma = new TClonesArray("AliGammaConversionAODObject", 0);
2674   else fAODGamma->Delete();
2675   fAODGamma->SetName(Form("%s_gamma", fAODBranchName.Data()));
2676
2677   if(!fAODPi0) fAODPi0 = new TClonesArray("AliGammaConversionAODObject", 0);
2678   else fAODPi0->Delete();
2679   fAODPi0->SetName(Form("%s_Pi0", fAODBranchName.Data()));
2680
2681   if(!fAODOmega) fAODOmega = new TClonesArray("AliGammaConversionAODObject", 0);
2682   else fAODOmega->Delete();
2683   fAODOmega->SetName(Form("%s_Omega", fAODBranchName.Data()));
2684
2685   //If delta AOD file name set, add in separate file. Else add in standard aod file.
2686   if(fKFDeltaAODFileName.Length() > 0) {
2687     AddAODBranch("TClonesArray", &fAODGamma, fKFDeltaAODFileName.Data());
2688     AddAODBranch("TClonesArray", &fAODPi0, fKFDeltaAODFileName.Data());
2689     AddAODBranch("TClonesArray", &fAODOmega, fKFDeltaAODFileName.Data());
2690     AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fKFDeltaAODFileName.Data());
2691   } else  {
2692     AddAODBranch("TClonesArray", &fAODGamma);
2693     AddAODBranch("TClonesArray", &fAODPi0);
2694     AddAODBranch("TClonesArray", &fAODOmega);
2695   }
2696
2697   // Create the output container
2698   if(fOutputContainer != NULL){
2699     delete fOutputContainer;
2700     fOutputContainer = NULL;
2701   }
2702   if(fOutputContainer == NULL){
2703     fOutputContainer = new TList();
2704     fOutputContainer->SetOwner(kTRUE);
2705   }
2706
2707   //Adding the histograms to the output container
2708   fHistograms->GetOutputContainer(fOutputContainer);
2709
2710
2711   if(fWriteNtuple){
2712     if(fGammaNtuple == NULL){
2713       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);
2714     }
2715     if(fNeutralMesonNtuple == NULL){
2716       fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
2717     }
2718     TList * ntupleTList = new TList();
2719     ntupleTList->SetOwner(kTRUE);
2720     ntupleTList->SetName("Ntuple");
2721     ntupleTList->Add((TNtuple*)fGammaNtuple);
2722     fOutputContainer->Add(ntupleTList);
2723   }
2724
2725   fOutputContainer->SetName(GetName());
2726 }
2727
2728 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{
2729   //helper function
2730   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
2731   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
2732   return v3D0.Angle(v3D1);
2733 }
2734
2735 void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
2736   // see header file for documentation
2737
2738   vector<Int_t> indexOfGammaParticle;
2739
2740   fStack = fV0Reader->GetMCStack();
2741
2742   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
2743     return; // aborts if the primary vertex does not have contributors.
2744   }
2745
2746   for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
2747     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
2748     if(particle->GetPdgCode()==22){     //Gamma
2749       if(particle->GetNDaughters() >= 2){
2750         TParticle* electron=NULL;
2751         TParticle* positron=NULL;
2752         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
2753           TParticle *tmpDaughter = fStack->Particle(daughterIndex);
2754           if(tmpDaughter->GetUniqueID() == 5){
2755             if(tmpDaughter->GetPdgCode() == 11){
2756               electron = tmpDaughter;
2757             }
2758             else if(tmpDaughter->GetPdgCode() == -11){
2759               positron = tmpDaughter;
2760             }
2761           }
2762         }
2763         if(electron!=NULL && positron!=0){
2764           if(electron->R()<160){
2765             indexOfGammaParticle.push_back(iTracks);
2766           }
2767         }
2768       }
2769     }
2770   }
2771
2772   Int_t nFoundGammas=0;
2773   Int_t nNotFoundGammas=0;
2774
2775   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
2776   for(Int_t i=0;i<numberOfV0s;i++){
2777     fV0Reader->GetV0(i);
2778
2779     if(fV0Reader->HasSameMCMother() == kFALSE){
2780       continue;
2781     }
2782
2783     TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2784     TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2785
2786     if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
2787       continue;
2788     }
2789     if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
2790       continue;
2791     }
2792
2793     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2794       //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
2795       for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
2796         if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
2797           nFoundGammas++;
2798         }
2799         else{
2800           nNotFoundGammas++;
2801         }
2802       }
2803     }
2804   }
2805 }
2806
2807
2808 void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
2809   // see header file for documantation
2810
2811   fESDEvent = fV0Reader->GetESDEvent();
2812
2813
2814   TClonesArray * vESDeNegTemp = new TClonesArray("AliESDtrack",0);
2815   TClonesArray * vESDePosTemp = new TClonesArray("AliESDtrack",0);
2816   TClonesArray * vESDxNegTemp = new TClonesArray("AliESDtrack",0);
2817   TClonesArray * vESDxPosTemp = new TClonesArray("AliESDtrack",0);
2818   TClonesArray * vESDeNegNoJPsi = new TClonesArray("AliESDtrack",0);
2819   TClonesArray * vESDePosNoJPsi = new TClonesArray("AliESDtrack",0);
2820
2821   /*
2822     vector <AliESDtrack*> vESDeNegTemp(0);
2823     vector <AliESDtrack*> vESDePosTemp(0);
2824     vector <AliESDtrack*> vESDxNegTemp(0);
2825     vector <AliESDtrack*> vESDxPosTemp(0);
2826     vector <AliESDtrack*> vESDeNegNoJPsi(0);
2827     vector <AliESDtrack*> vESDePosNoJPsi(0);
2828   */
2829
2830
2831   fHistograms->FillTable("Table_Electrons",0);//Count number of Events
2832
2833   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2834     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
2835
2836     if(!curTrack){
2837       //print warning here
2838       continue;
2839     }
2840
2841     double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
2842     double r[3];curTrack->GetConstrainedXYZ(r);
2843
2844     TVector3 rXYZ(r);
2845
2846     fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
2847
2848     Bool_t flagKink       =  kTRUE;
2849     Bool_t flagTPCrefit   =  kTRUE;
2850     Bool_t flagTRDrefit   =  kTRUE;
2851     Bool_t flagITSrefit   =  kTRUE;
2852     Bool_t flagTRDout     =  kTRUE;
2853     Bool_t flagVertex     =  kTRUE;
2854
2855
2856     //Cuts ---------------------------------------------------------------
2857
2858     if(curTrack->GetKinkIndex(0) > 0){
2859       fHistograms->FillHistogram("Table_Electrons",5);//Count kink
2860       flagKink = kFALSE;
2861     }
2862
2863     ULong_t trkStatus = curTrack->GetStatus();
2864
2865     ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
2866
2867     if(!tpcRefit){
2868       fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
2869       flagTPCrefit = kFALSE;
2870     }
2871
2872     ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
2873     if(!itsRefit){
2874       fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
2875       flagITSrefit = kFALSE;
2876     }
2877
2878     ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
2879
2880     if(!trdRefit){
2881       fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
2882       flagTRDrefit = kFALSE;
2883     }
2884
2885     ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
2886
2887     if(!trdOut) {
2888       fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
2889       flagTRDout = kFALSE;
2890     }
2891
2892     double nsigmaToVxt = GetSigmaToVertex(curTrack);
2893
2894     if(nsigmaToVxt > 3){
2895       fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
2896       flagVertex = kFALSE;
2897     }
2898
2899     if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
2900     fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
2901
2902
2903     Stat_t pid, weight;
2904     GetPID(curTrack, pid, weight);
2905
2906     if(pid!=0){
2907       fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
2908     }
2909
2910     if(pid == 0){
2911       fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
2912     }
2913
2914
2915
2916
2917
2918
2919     TLorentzVector curElec;
2920     curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
2921
2922
2923     if(fDoMCTruth){
2924       Int_t labelMC = TMath::Abs(curTrack->GetLabel());
2925       TParticle* curParticle = fStack->Particle(labelMC);
2926       if(curTrack->GetSign() > 0){
2927         if( pid == 0){
2928           fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
2929           fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
2930         }
2931         else{
2932           fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
2933           fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
2934         }
2935       }
2936     }
2937
2938
2939     if(curTrack->GetSign() > 0){
2940
2941       //     vESDxPosTemp.push_back(curTrack);
2942       new((*vESDxPosTemp)[vESDxPosTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
2943
2944       if( pid == 0){
2945
2946         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
2947         fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
2948         //      fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
2949         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
2950         //      fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
2951         //      vESDePosTemp.push_back(curTrack);
2952         new((*vESDePosTemp)[vESDePosTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
2953
2954       }
2955
2956     }
2957     else {
2958
2959       new((*vESDxNegTemp)[vESDxNegTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
2960
2961       if( pid == 0){
2962
2963         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
2964         fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
2965         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
2966         new((*vESDeNegTemp)[vESDeNegTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
2967
2968       }
2969
2970     }
2971
2972   }
2973
2974
2975   Bool_t ePosJPsi = kFALSE;
2976   Bool_t eNegJPsi = kFALSE;
2977   Bool_t ePosPi0  = kFALSE;
2978   Bool_t eNegPi0  = kFALSE;
2979
2980   UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
2981
2982   for(Int_t iNeg=0; iNeg < vESDeNegTemp->GetEntriesFast(); iNeg++){
2983     if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetPdgCode() == 11)
2984       if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0) > -1){
2985         Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0);
2986         TParticle* partMother = fStack ->Particle(labelMother);
2987         if (partMother->GetPdgCode() == 111){
2988           ieNegPi0 = iNeg;
2989           eNegPi0 = kTRUE;
2990         }
2991         if(partMother->GetPdgCode() == 443){ //Mother JPsi
2992           fHistograms->FillTable("Table_Electrons",14);
2993           ieNegJPsi = iNeg;
2994           eNegJPsi = kTRUE;
2995         }
2996         else{
2997           //      vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
2998           new((*vESDeNegNoJPsi)[vESDeNegNoJPsi->GetEntriesFast()])  AliESDtrack(*(AliESDtrack*)(vESDeNegTemp->At(iNeg)));
2999           //            cout<<"ESD No Positivo JPsi "<<endl;
3000         }
3001
3002       }
3003   }
3004
3005   for(Int_t iPos=0; iPos < vESDePosTemp->GetEntriesFast(); iPos++){
3006     if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetPdgCode() == -11)
3007       if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0) > -1){
3008         Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0);
3009         TParticle* partMother = fStack ->Particle(labelMother);
3010         if (partMother->GetPdgCode() == 111){
3011           iePosPi0 = iPos;
3012           ePosPi0 = kTRUE;
3013         }
3014         if(partMother->GetPdgCode() == 443){ //Mother JPsi
3015           fHistograms->FillTable("Table_Electrons",15);
3016           iePosJPsi = iPos;
3017           ePosJPsi = kTRUE;
3018         }
3019         else{
3020           //      vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
3021           new((*vESDePosNoJPsi)[vESDePosNoJPsi->GetEntriesFast()])  AliESDtrack(*(AliESDtrack*)(vESDePosTemp->At(iPos)));
3022           //            cout<<"ESD No Negativo JPsi "<<endl;
3023         }
3024
3025       }
3026   }
3027
3028   if( eNegJPsi && ePosJPsi ){
3029     TVector3 tempeNegV,tempePosV;
3030     tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Pz());
3031     tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Pz());
3032     fHistograms->FillTable("Table_Electrons",16);
3033     fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));
3034     fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->GetLabel())),
3035                                                                               fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->GetLabel()))));
3036   }
3037
3038   if( eNegPi0 && ePosPi0 ){
3039     TVector3 tempeNegV,tempePosV;
3040     tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Pz());
3041     tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Pz());
3042     fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
3043     fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->GetLabel())),
3044                                                                              fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->GetLabel()))));
3045   }
3046
3047
3048   FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
3049
3050   CleanWithAngleCuts(*vESDeNegTemp,*vESDePosTemp,*fKFReconstructedGammasTClone);
3051
3052   //  vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
3053   //  vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
3054
3055   TClonesArray vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectronTClone);
3056   TClonesArray vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectronTClone);
3057
3058
3059   FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
3060
3061
3062
3063
3064   //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
3065
3066
3067   FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
3068   FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
3069
3070
3071
3072   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",*fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
3073
3074   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
3075                            *fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
3076
3077   //BackGround
3078
3079   //Like Sign e+e-
3080   ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
3081   ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
3082   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
3083   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
3084
3085   //        Like Sign e+e- no JPsi
3086   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
3087   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
3088
3089   //Mixed Event
3090
3091   if( fCurrentEventPosElectronTClone->GetEntriesFast() > 0 && fCurrentEventNegElectronTClone->GetEntriesFast() > 0 && fKFReconstructedGammasCutTClone->GetEntriesFast() > 0 ){
3092     FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
3093                              *fKFReconstructedGammasCutTClone,*fPreviousEventTLVNegElectronTClone,*fPreviousEventTLVPosElectronTClone);
3094     *fPreviousEventTLVNegElectronTClone = vCurrentTLVeNeg;
3095     *fPreviousEventTLVPosElectronTClone = vCurrentTLVePos;
3096
3097   }
3098
3099   /*
3100   //Photons P
3101   Double_t vtx[3];
3102   vtx[0]=0;vtx[1]=0;vtx[2]=0;
3103   for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
3104
3105   //      if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
3106
3107
3108
3109   Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
3110   //      cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
3111
3112   //      cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
3113
3114   if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
3115
3116   fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
3117
3118
3119   }
3120
3121
3122   */
3123
3124
3125   vESDeNegTemp->Delete();
3126   vESDePosTemp->Delete();
3127   vESDxNegTemp->Delete();
3128   vESDxPosTemp->Delete();
3129   vESDeNegNoJPsi->Delete();
3130   vESDePosNoJPsi->Delete();
3131
3132   delete vESDeNegTemp;
3133   delete vESDePosTemp;
3134   delete vESDxNegTemp;
3135   delete vESDxPosTemp;
3136   delete vESDeNegNoJPsi;
3137   delete vESDePosNoJPsi;
3138 }
3139
3140 /*
3141   void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
3142   //see header file for documentation
3143   for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
3144   for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
3145   fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
3146   }
3147   }
3148   }
3149 */
3150 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,TClonesArray const tlVeNeg, TClonesArray const tlVePos){
3151   //see header file for documentation
3152   for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++){
3153     for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
3154       fHistograms->FillHistogram(histoName.Data(),((TLorentzVector*)(tlVeNeg.At(iNeg)))->Vect().Angle(((TLorentzVector*)(tlVePos.At(iPos)))->Vect()));
3155     }
3156   }
3157 }
3158 void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, TClonesArray const eNeg, TClonesArray const ePos){
3159   //see header file for documentation
3160   for( Int_t n=0; n < eNeg.GetEntriesFast(); n++){
3161     TLorentzVector en = (*(TLorentzVector*)(eNeg.At(n)));
3162     for (Int_t p=0; p < ePos.GetEntriesFast(); p++){
3163       TLorentzVector ep = (*(TLorentzVector*)(ePos.At(p)));
3164       TLorentzVector np = ep + en;
3165       fHistograms->FillHistogram(histoName.Data(),np.M());
3166     }
3167   }
3168 }
3169
3170 void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,TClonesArray const fKFGammas,
3171                                                               TClonesArray const tlVeNeg,TClonesArray const tlVePos)
3172 {
3173   //see header file for documentation
3174
3175   for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++ ){
3176
3177     for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
3178
3179       TLorentzVector xy = *((TLorentzVector *)(tlVePos.At(iPos))) + *((TLorentzVector *)(tlVeNeg.At(iNeg)));
3180
3181       for (Int_t iGam=0; iGam < fKFGammas.GetEntriesFast(); iGam++){
3182
3183         //      AliKFParticle * gammaCandidate = &fKFGammas[iGam];
3184         AliKFParticle * gammaCandidate = (AliKFParticle *)(fKFGammas.At(iGam));
3185         TLorentzVector g;
3186
3187         g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
3188         TLorentzVector xyg = xy + g;
3189         fHistograms->FillHistogram(histoMass.Data(),xyg.M());
3190         fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
3191       }
3192     }
3193   }
3194
3195 }
3196 void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, TClonesArray e)
3197 {
3198   // see header file for documentation
3199   for(Int_t i=0; i < e.GetEntriesFast(); i++)
3200     {
3201       for (Int_t j=i+1; j < e.GetEntriesFast(); j++)
3202         {
3203           TLorentzVector ee = (*(TLorentzVector*)(e.At(i))) + (*(TLorentzVector*)(e.At(j)));
3204
3205           fHistograms->FillHistogram(hBg.Data(),ee.M());
3206         }
3207     }
3208 }
3209
3210
3211 void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negativeElectrons,
3212                                                         TClonesArray const positiveElectrons,
3213                                                         TClonesArray const gammas){
3214   // see header file for documentation
3215
3216   UInt_t  sizeN = negativeElectrons.GetEntriesFast();
3217   UInt_t  sizeP = positiveElectrons.GetEntriesFast();
3218   UInt_t  sizeG = gammas.GetEntriesFast();
3219
3220
3221
3222   vector <Bool_t> xNegBand(sizeN);
3223   vector <Bool_t> xPosBand(sizeP);
3224   vector <Bool_t> gammaBand(sizeG);
3225
3226
3227   for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
3228   for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
3229   for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
3230
3231
3232   for(UInt_t iPos=0; iPos < sizeP; iPos++){
3233
3234     Double_t aP[3];
3235     ((AliESDtrack*)(positiveElectrons.At(iPos)))->GetConstrainedPxPyPz(aP);
3236
3237     TVector3 ePosV(aP[0],aP[1],aP[2]);
3238
3239     for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
3240
3241       Double_t aN[3];
3242       ((AliESDtrack*)(negativeElectrons.At(iNeg)))->GetConstrainedPxPyPz(aN);
3243       TVector3 eNegV(aN[0],aN[1],aN[2]);
3244
3245       if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
3246         xPosBand[iPos]=kFALSE;
3247         xNegBand[iNeg]=kFALSE;
3248       }
3249
3250       for(UInt_t iGam=0; iGam < sizeG; iGam++){
3251         AliKFParticle* gammaCandidate = (AliKFParticle*)gammas.At(iGam);
3252         TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
3253         if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
3254           gammaBand[iGam]=kFALSE;
3255       }
3256     }
3257   }
3258
3259
3260
3261
3262   for(UInt_t iPos=0; iPos < sizeP; iPos++){
3263     if(xPosBand[iPos]){
3264       new((*fCurrentEventPosElectronTClone)[fCurrentEventPosElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(positiveElectrons.At(iPos))));
3265       //      fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
3266     }
3267   }
3268   for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
3269     if(xNegBand[iNeg]){
3270       new((*fCurrentEventNegElectronTClone)[fCurrentEventNegElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(negativeElectrons.At(iNeg))));
3271       //      fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
3272     }
3273   }
3274   for(UInt_t iGam=0; iGam < sizeG; iGam++){
3275     if(gammaBand[iGam]){
3276       new((*fKFReconstructedGammasCutTClone)[fKFReconstructedGammasCutTClone->GetEntriesFast()]) AliKFParticle((*(AliKFParticle*)(gammas.At(iGam))));
3277       //fKFReconstructedGammasCut.push_back(*(AliKFParticle*)gammas->At(iGam));
3278     }
3279   }
3280 }
3281
3282
3283 void  AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)
3284 {
3285   // see header file for documentation
3286   pid = -1;
3287   weight = -1;
3288
3289   double wpart[5];
3290   double wpartbayes[5];
3291
3292   //get probability of the diffenrent particle types
3293   track->GetESDpid(wpart);
3294
3295   // Tentative particle type "concentrations"
3296   double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
3297
3298   //Bayes' formula
3299   double rcc = 0.;
3300   for (int i = 0; i < 5; i++)
3301     {
3302       rcc+=(c[i] * wpart[i]);
3303     }
3304
3305
3306
3307   for (int i=0; i<5; i++) {
3308     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)
3309       wpartbayes[i] = c[i] * wpart[i] / rcc;
3310     }
3311   }
3312
3313
3314
3315   Float_t max=0.;
3316   int ipid=-1;
3317   //find most probable particle in ESD pid
3318   //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
3319   for (int i = 0; i < 5; i++)
3320     {
3321       if (wpartbayes[i] > max)
3322         {
3323           ipid = i;
3324           max = wpartbayes[i];
3325         }
3326     }
3327
3328   pid = ipid;
3329   weight = max;
3330 }
3331 double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)
3332 {
3333   // Calculates the number of sigma to the vertex.
3334
3335   Float_t b[2];
3336   Float_t bRes[2];
3337   Float_t bCov[3];
3338   t->GetImpactParameters(b,bCov);
3339   if (bCov[0]<=0 || bCov[2]<=0) {
3340     AliDebug(1, "Estimated b resolution lower or equal zero!");
3341     bCov[0]=0; bCov[2]=0;
3342   }
3343   bRes[0] = TMath::Sqrt(bCov[0]);
3344   bRes[1] = TMath::Sqrt(bCov[2]);
3345
3346   // -----------------------------------
3347   // How to get to a n-sigma cut?
3348   //
3349   // The accumulated statistics from 0 to d is
3350   //
3351   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
3352   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
3353   //
3354   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
3355   // Can this be expressed in a different way?
3356
3357   if (bRes[0] == 0 || bRes[1] ==0)
3358     return -1;
3359
3360   double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
3361
3362   // stupid rounding problem screws up everything:
3363   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
3364   if (TMath::Exp(-d * d / 2) < 1e-10)
3365     return 1000;
3366
3367
3368   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
3369   return d;
3370 }
3371
3372 //vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
3373 TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *const esdTrack){
3374   //Return TLoresntz vector of track?
3375   //  vector <TLorentzVector> tlVtrack(0);
3376   TClonesArray array("TLorentzVector",0);
3377
3378   for(Int_t itrack=0; itrack < esdTrack->GetEntriesFast(); itrack++){
3379     double p[3];
3380     //esdTrack[itrack]->GetConstrainedPxPyPz(p);
3381     ((AliESDtrack*)(esdTrack->At(itrack)))->GetConstrainedPxPyPz(p);
3382     TLorentzVector currentTrack;
3383     currentTrack.SetXYZM(p[0],p[1],p[2],fElectronMass);
3384     new((array)[array.GetEntriesFast()])  TLorentzVector(currentTrack);
3385     //    tlVtrack.push_back(currentTrack);
3386   }
3387
3388   return array;
3389
3390   //  return tlVtrack;
3391 }
3392