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