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