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