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