PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx

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