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