PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx

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