PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx

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