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