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::GetStandardITSTPCTrackCuts2010(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             if(massTwoGammaCandidate>0.1 && massTwoGammaCandidate<0.15){
2057               fHistograms->FillHistogram("ESD_Mother_alfa_Pi0", alfa);
2058             }
2059             fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!
2060             fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
2061             fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
2062             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2063             fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
2064             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2065           }
2066           if(alfa<0.1){
2067             fHistograms->FillHistogram("ESD_Mother_InvMass_vs_E_alpha",massTwoGammaCandidate ,twoGammaCandidate->GetE());
2068           }
2069
2070
2071           if(fCalculateBackground){
2072             /* Kenneth, just for testing*/
2073             AliGammaConversionBGHandler * bgHandlerTest = fV0Reader->GetBGHandler();
2074
2075             Int_t zbin= bgHandlerTest->GetZBinIndex(fV0Reader->GetVertexZ());
2076             Int_t mbin=0;
2077             Int_t multKAA=0;
2078             if(fUseTrackMultiplicityForBG == kTRUE){
2079               multKAA=fV0Reader->CountESDTracks();
2080               mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2081             }
2082             else{// means we use #v0s for multiplicity
2083               multKAA=fV0Reader->GetNGoodV0s();
2084               mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2085             }
2086             //      cout<<"Filling bin number "<<zbin<<" and "<<mbin<<endl;
2087             //      cout<<"Corresponding to z = "<<fV0Reader->GetVertexZ()<<" and m = "<<multKAA<<endl;
2088             if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2089               fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass",zbin,mbin),massTwoGammaCandidate);
2090               fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass_vs_Pt",zbin,mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2091               /* end Kenneth, just for testing*/
2092               fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2093             }
2094           }
2095           /*      if(fCalculateBackground){
2096             AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2097             Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2098             fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2099             }*/
2100           //      if(fDoNeutralMesonV0MCCheck){
2101           if(fDoMCTruth){
2102             //Kenneth: Checking the eta of the gamma to check the difference between 0.9 and 1.2
2103             Int_t indexKF1 = fKFReconstructedGammasV0Index.at(firstGammaIndex);
2104             if(indexKF1<fV0Reader->GetNumberOfV0s()){
2105               fV0Reader->GetV0(indexKF1);//updates to the correct v0
2106               Double_t eta1 = fV0Reader->GetMotherCandidateEta();
2107               Bool_t isRealPi0=kFALSE;
2108               Bool_t isRealEta=kFALSE;
2109               Int_t gamma1MotherLabel=-1;
2110               if(fV0Reader->HasSameMCMother() == kTRUE){
2111                 //cout<<"This v0 is a real v0!!!!"<<endl;
2112                 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2113                 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2114                 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2115                   if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2116                     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2117                       gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2118                     }
2119                   }
2120                 }
2121               }
2122               Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
2123               if(indexKF1 == indexKF2){
2124                 cout<<"index of the two KF particles are the same.... should not happen"<<endl;
2125               }
2126               if(indexKF2<fV0Reader->GetNumberOfV0s()){
2127                 fV0Reader->GetV0(indexKF2);
2128                 Double_t eta2 = fV0Reader->GetMotherCandidateEta();
2129                 Int_t gamma2MotherLabel=-1;
2130                 if(fV0Reader->HasSameMCMother() == kTRUE){
2131                   TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2132                   TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2133                   if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2134                     if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2135                       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2136                         gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2137                       }
2138                     }
2139                   }
2140                 }
2141                 if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
2142                   if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
2143                     isRealPi0=kTRUE;
2144                   }
2145                   if(fV0Reader->CheckIfEtaIsMother(gamma1MotherLabel)){
2146                     isRealEta=kTRUE;
2147                   }
2148
2149                 }
2150                 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2151                   if(TMath::Abs(eta1)>0.9 && TMath::Abs(eta2)>0.9){
2152                     //            fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2153                     //            fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2154                     if(isRealPi0 || isRealEta){
2155                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_1212",massTwoGammaCandidate);
2156                       fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_1212",openingAngleTwoGammaCandidate);
2157                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2158                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2159                       fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
2160                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2161                     }
2162
2163                     if(!isRealPi0 && !isRealEta){
2164                       if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2165                         fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2166                       }else{
2167                         fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2168                       }
2169                     }
2170
2171                   }
2172                   else if(TMath::Abs(eta1)>0.9 || TMath::Abs(eta2)>0.9){
2173                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2174                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2175
2176                     if(isRealPi0 || isRealEta){
2177                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_0912",massTwoGammaCandidate);
2178                       fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0912",openingAngleTwoGammaCandidate);
2179                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2180                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2181                       fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
2182                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2183                     }
2184                     if(!isRealPi0 && !isRealEta){
2185                       if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2186                         fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2187                       }else{
2188                         fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2189                       }
2190                     }
2191                   }
2192                   else{
2193                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2194                   //              fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2195                     if(isRealPi0 || isRealEta){
2196                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_0909",massTwoGammaCandidate);
2197                       fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0909",openingAngleTwoGammaCandidate);
2198                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2199                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2200                       fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
2201                       fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2202                       if(gamma1MotherLabel > fV0Reader->GetMCStack()->GetNprimary()){
2203                         fHistograms->FillHistogram("ESD_TruePi0Sec_InvMass",massTwoGammaCandidate);
2204                       }
2205                     }
2206                     if(!isRealPi0 && !isRealEta){
2207                       if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2208                         fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2209                       }else{
2210                         fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2211                       }
2212                     }
2213                   }
2214                 }
2215               }
2216             }
2217           }
2218           if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2219             if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 &&  TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
2220               fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2221               fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
2222             }
2223
2224             if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2225               fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2226               fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2227             }
2228             else if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 || TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2229               fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2230               fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2231             }
2232             else{
2233               fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2234               fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2235             }
2236
2237             Double_t lowMassPi0=0.1;
2238             Double_t highMassPi0=0.15;
2239             if (massTwoGammaCandidate > lowMassPi0 && massTwoGammaCandidate < highMassPi0 ){
2240               new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()])  AliKFParticle(*twoGammaCandidate);
2241               fGammav1.push_back(firstGammaIndex);
2242               fGammav2.push_back(secondGammaIndex);
2243               AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
2244             }
2245           }
2246
2247         }
2248           //}
2249         delete twoGammaCandidate;
2250     }
2251   }
2252 }
2253
2254 void AliAnalysisTaskGammaConversion::AddPionToAOD(AliKFParticle * pionkf, Double_t mass, Int_t daughter1, Int_t daughter2) {
2255   //See header file for documentation
2256   AliGammaConversionAODObject pion;
2257   pion.SetPx(pionkf->Px());
2258   pion.SetPy(pionkf->Py());
2259   pion.SetPz(pionkf->Pz());
2260   pion.SetChi2(pionkf->GetChi2());
2261   pion.SetE(pionkf->E());
2262   pion.SetIMass(mass);
2263   pion.SetLabel1(daughter1);
2264   //dynamic_cast<AliGammaConversionAODObject*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
2265   pion.SetLabel2(daughter2);
2266   new((*fAODPi0)[fAODPi0->GetEntriesFast()])  AliGammaConversionAODObject(pion);
2267
2268 }
2269
2270
2271
2272 void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
2273 /*
2274   // see header file for documentation
2275   // Analyse Pi0 with one photon from Phos and 1 photon from conversions
2276
2277
2278
2279   Double_t vtx[3];
2280   vtx[0] = fV0Reader->GetPrimaryVertex()->GetX();
2281   vtx[1] = fV0Reader->GetPrimaryVertex()->GetY();
2282   vtx[2] = fV0Reader->GetPrimaryVertex()->GetZ();
2283
2284
2285   // Loop over all CaloClusters and consider only the PHOS ones:
2286   AliESDCaloCluster *clu;
2287   TLorentzVector pPHOS;
2288   TLorentzVector gammaPHOS;
2289   TLorentzVector gammaGammaConv;
2290   TLorentzVector pi0GammaConvPHOS;
2291   TLorentzVector gammaGammaConvBck;
2292   TLorentzVector pi0GammaConvPHOSBck;
2293
2294
2295   for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2296     clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2297     if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
2298     clu ->GetMomentum(pPHOS ,vtx);
2299     for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2300       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2301       gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
2302       gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
2303       pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
2304       fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
2305       fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
2306
2307       TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
2308       TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
2309       Double_t opanConvPHOS= v3D0.Angle(v3D1);
2310       if ( opanConvPHOS < 0.35){
2311         fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
2312       }else{
2313         fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
2314       }
2315
2316     }
2317
2318     //  Now the LorentVector pPHOS is obtained and can be paired with the converted proton
2319   }
2320   //==== End of the PHOS cluster selection ============
2321   TLorentzVector pEMCAL;
2322   TLorentzVector gammaEMCAL;
2323   TLorentzVector pi0GammaConvEMCAL;
2324   TLorentzVector pi0GammaConvEMCALBck;
2325
2326    for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2327     clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2328     if ( !clu->IsEMCAL()  || clu->E()<0.1 ) continue;
2329     if (clu->GetNCells() <= 1) continue;
2330     if ( clu->GetTOF()*1e9 < 550  || clu->GetTOF()*1e9 > 750) continue;
2331
2332     clu ->GetMomentum(pEMCAL ,vtx);
2333     for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2334       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2335       gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
2336                              twoGammaDecayCandidateDaughter0->Py(),
2337                              twoGammaDecayCandidateDaughter0->Pz(),0.);
2338       gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
2339       pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
2340       fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
2341       fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
2342       TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
2343                     twoGammaDecayCandidateDaughter0->Py(),
2344                     twoGammaDecayCandidateDaughter0->Pz());
2345       TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
2346
2347
2348       Double_t opanConvEMCAL= v3D0.Angle(v3D1);
2349       if ( opanConvEMCAL < 0.35){
2350         fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
2351       }else{
2352         fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
2353       }
2354
2355     }
2356     if(fCalculateBackground){
2357       for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2358         AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2359         for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2360           AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2361           gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
2362                                     previousGoodV0.Py(),
2363                                     previousGoodV0.Pz(),0.);
2364           pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
2365           fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
2366           fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
2367                                      pi0GammaConvEMCALBck.Pt());
2368         }
2369       }
2370
2371       //  Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
2372     } // end of checking if background photons are available
2373    }
2374   //==== End of the PHOS cluster selection ============
2375 */
2376 }
2377
2378 void AliAnalysisTaskGammaConversion::MoveParticleAccordingToVertex(AliKFParticle *particle,AliGammaConversionBGHandler::GammaConversionVertex *vertex){
2379   //see header file for documentation
2380
2381   Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
2382   Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
2383   Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
2384
2385   //  cout<<"dx, dy, dz: ["<<dx<<","<<dy<<","<<dz<<"]"<<endl;
2386   particle->X() = particle->GetX() - dx;
2387   particle->Y() = particle->GetY() - dy;
2388   particle->Z() = particle->GetZ() - dz;
2389 }
2390
2391 void AliAnalysisTaskGammaConversion::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle){
2392
2393   Double_t c = cos(angle);
2394   Double_t s = sin(angle);
2395
2396   Double_t A[7][ 7];
2397   for( Int_t i=0; i<7; i++ ){
2398     for( Int_t j=0; j<7; j++){
2399       A[i][j] = 0;
2400     }
2401   }
2402   for( int i=0; i<7; i++ ){
2403     A[i][i] = 1;
2404   }
2405   A[0][0] =  c;   A[0][1] = s;
2406   A[1][0] = -s;  A[1][1] = c;
2407   A[3][3] =  c;  A[3][4] = s;
2408   A[4][3] = -s;  A[4][4] = c;
2409
2410   Double_t AC[7][7];
2411   Double_t Ap[7];
2412
2413   for( Int_t i=0; i<7; i++ ){
2414     Ap[i] = 0;
2415     for( Int_t k=0; k<7; k++){
2416       Ap[i]+=A[i][k] * kfParticle->GetParameter(k);
2417     }
2418   }
2419
2420   for( Int_t i=0; i<7; i++){
2421     kfParticle->Parameter(i) = Ap[i];
2422   }
2423
2424   for( Int_t i=0; i<7; i++ ){
2425     for( Int_t j=0; j<7; j++ ){
2426       AC[i][j] = 0;
2427       for( Int_t k=0; k<7; k++ ){
2428         AC[i][j]+= A[i][k] * kfParticle->GetCovariance(k,j);
2429       }
2430     }
2431   }
2432
2433   for( Int_t i=0; i<7; i++ ){
2434     for( Int_t j=0; j<=i; j++ ){
2435       Double_t xx = 0;
2436       for( Int_t k=0; k<7; k++){
2437         xx+= AC[i][k]*A[j][k];
2438       }
2439       kfParticle->Covariance(i,j) = xx;
2440     }
2441   }
2442 }
2443
2444
2445 void AliAnalysisTaskGammaConversion::CalculateBackground(){
2446   // see header file for documentation
2447
2448
2449   TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
2450
2451   AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2452
2453   Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2454   Int_t mbin = 0;
2455   if(fUseTrackMultiplicityForBG == kTRUE){
2456     mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2457   }
2458   else{
2459     mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2460   }
2461
2462   if(fDoRotation == kTRUE){
2463     TRandom3 *random = new TRandom3();
2464     for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2465       AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2466       for(Int_t iCurrent2=iCurrent+1;iCurrent2<currentEventV0s->GetEntriesFast();iCurrent2++){
2467         for(Int_t nRandom=0;nRandom<nRandomEventsForBG;nRandom++){
2468
2469           AliKFParticle currentEventGoodV02 = *(AliKFParticle *)(currentEventV0s->At(iCurrent2));
2470
2471           if(fCheckBGProbability == kTRUE){
2472             Double_t massBGprob =0.;
2473             Double_t widthBGprob = 0.;
2474             AliKFParticle *backgroundCandidateProb = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
2475             backgroundCandidateProb->GetMass(massBGprob,widthBGprob);
2476             if(massBGprob>0.1 && massBGprob<0.14){
2477               if(random->Rndm()>bgHandler->GetBGProb(zbin,mbin)){
2478                 delete backgroundCandidateProb;
2479                 continue;
2480               }
2481             }
2482             delete backgroundCandidateProb;
2483           }
2484
2485           Double_t nRadiansPM = nDegreesPMBackground*TMath::Pi()/180;
2486
2487           Double_t rotationValue = random->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
2488
2489           RotateKFParticle(&currentEventGoodV02,rotationValue);
2490
2491           AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
2492
2493           Double_t massBG =0.;
2494           Double_t widthBG = 0.;
2495           Double_t chi2BG =10000.;
2496           backgroundCandidate->GetMass(massBG,widthBG);
2497
2498           //      if(backgroundCandidate->GetNDF()>0){
2499           chi2BG = backgroundCandidate->GetChi2();
2500           if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson())  || fApplyChi2Cut == kFALSE){
2501
2502             TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2503             TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2504
2505             Double_t openingAngleBG = currentEventGoodV0.GetAngle(currentEventGoodV02);
2506
2507             Double_t rapidity;
2508             if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2509             else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2510
2511             if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2512               delete backgroundCandidate;
2513               continue;   // rapidity cut
2514             }
2515
2516
2517             Double_t alfa=0.0;
2518             if( (currentEventGoodV0.GetE()+currentEventGoodV02.GetE()) != 0){
2519               alfa=TMath::Abs((currentEventGoodV0.GetE()-currentEventGoodV02.GetE())
2520                               /(currentEventGoodV0.GetE()+currentEventGoodV02.GetE()));
2521             }
2522
2523
2524             if(openingAngleBG < fMinOpeningAngleGhostCut ){
2525               delete backgroundCandidate;
2526               continue;   // minimum opening angle to avoid using ghosttracks
2527             }
2528
2529             // original
2530             if(alfa>fV0Reader->GetAlphaMinCutMeson() &&   alfa<fV0Reader->GetAlphaCutMeson()){
2531               fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2532               fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2533               fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
2534               fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2535               fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2536               fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2537               fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2538               fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2539               fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2540               fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2541               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2542               fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2543               fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2544
2545
2546               if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2547                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2548                 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2549               }
2550
2551               fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2552               fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2553               fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
2554               fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2555               fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2556               fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2557               fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2558               fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2559               fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2560               fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2561               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2562               fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2563
2564               if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2565                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2566                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2567               }
2568             }
2569             if(alfa<0.1){
2570               fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2571             }
2572
2573           }
2574           //}
2575           delete backgroundCandidate;
2576         }
2577       }
2578     }
2579     delete random;
2580   }
2581   else{ // means no rotation
2582     AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2583
2584     if(fUseTrackMultiplicityForBG){
2585       //    cout<<"Using charged track multiplicity for background calculation"<<endl;
2586       for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2587
2588         AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);//fV0Reader->GetBGGoodV0s(nEventsInBG);
2589
2590         if(fMoveParticleAccordingToVertex == kTRUE){
2591           bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2592         }
2593
2594         for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2595           AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2596           for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2597             AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2598             AliKFParticle previousGoodV0test = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2599
2600             //cout<<"Primary Vertex event: ["<<fESDEvent->GetPrimaryVertex()->GetX()<<","<<fESDEvent->GetPrimaryVertex()->GetY()<<","<<fESDEvent->GetPrimaryVertex()->GetZ()<<"]"<<endl;
2601             //cout<<"BG prim Vertex event: ["<<bgEventVertex->fX<<","<<bgEventVertex->fY<<","<<bgEventVertex->fZ<<"]"<<endl;
2602
2603             //cout<<"XYZ of particle before transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
2604             if(fMoveParticleAccordingToVertex == kTRUE){
2605               MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2606             }
2607             //cout<<"XYZ of particle after transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
2608
2609             AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2610
2611             Double_t massBG =0.;
2612             Double_t widthBG = 0.;
2613             Double_t chi2BG =10000.;
2614             backgroundCandidate->GetMass(massBG,widthBG);
2615             //  if(backgroundCandidate->GetNDF()>0){
2616             //    chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2617             chi2BG = backgroundCandidate->GetChi2();
2618             if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2619
2620               TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2621               TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2622
2623               Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2624
2625               Double_t rapidity;
2626
2627               if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() <= 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() <= 0){
2628                 cout << "Error: |Pz| > E !!!! " << endl;
2629                 rapidity=0;
2630               } else {
2631                 rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2632               }
2633               if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2634                 delete backgroundCandidate;
2635                 continue;   // rapidity cut
2636               }
2637
2638
2639               Double_t alfa=0.0;
2640               if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2641                 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2642                                 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2643               }
2644
2645
2646               if(openingAngleBG < fMinOpeningAngleGhostCut ){
2647                 delete backgroundCandidate;
2648                 continue;   // minimum opening angle to avoid using ghosttracks
2649               }
2650
2651               // original
2652               if(alfa>fV0Reader->GetAlphaMinCutMeson() &&   alfa<fV0Reader->GetAlphaCutMeson()){
2653                 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2654                 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2655                 fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
2656                 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2657                 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2658                 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2659                 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2660                 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2661                 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2662                 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2663                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2664                 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2665                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2666
2667
2668                 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2669                   fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2670                   fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2671                 }
2672
2673                 // test
2674                 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2675                 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2676                 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
2677                 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2678                 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2679                 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2680                 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2681                 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2682                 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2683                 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2684                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2685                 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2686
2687                 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2688                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2689                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2690                 }
2691                 //        }
2692               }
2693               if(alfa<0.1){
2694                 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2695               }
2696
2697             }
2698             delete backgroundCandidate;
2699           }
2700         }
2701       }
2702     }
2703     else{ // means using #V0s for multiplicity
2704
2705       //    cout<<"Using the v0 multiplicity to calculate background"<<endl;
2706
2707       fHistograms->FillHistogram("ESD_Background_z_m",zbin,mbin);
2708       fHistograms->FillHistogram("ESD_Mother_multpilicityVSv0s",fV0Reader->CountESDTracks(),fV0Reader->GetNumberOfV0s());
2709
2710       for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2711         AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);// fV0Reader->GetBGGoodV0s(nEventsInBG);
2712         if(previousEventV0s){
2713
2714           if(fMoveParticleAccordingToVertex == kTRUE){
2715             bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2716           }
2717
2718           for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2719             AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2720             for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2721               AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2722
2723               if(fMoveParticleAccordingToVertex == kTRUE){
2724                 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2725               }
2726
2727               AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2728               Double_t massBG =0.;
2729               Double_t widthBG = 0.;
2730               Double_t chi2BG =10000.;
2731               backgroundCandidate->GetMass(massBG,widthBG);
2732               /*            if(backgroundCandidate->GetNDF()>0){
2733                             chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2734                             {//remember to remove
2735                             TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2736                             TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2737
2738                             Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2739                             fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle_nochi2", openingAngleBG);
2740                             }
2741               */
2742               chi2BG = backgroundCandidate->GetChi2();
2743               if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2744                 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2745                 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2746
2747                 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2748
2749                 Double_t rapidity;
2750                 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2751                 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2752
2753                 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2754                   delete backgroundCandidate;
2755                   continue;   // rapidity cut
2756                 }
2757
2758
2759                 Double_t alfa=0.0;
2760                 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2761                   alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2762                                   /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2763                 }
2764
2765
2766                 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2767                   delete backgroundCandidate;
2768                   continue;   // minimum opening angle to avoid using ghosttracks
2769                 }
2770
2771                 if(alfa>fV0Reader->GetAlphaMinCutMeson() &&   alfa<fV0Reader->GetAlphaCutMeson()){
2772                   fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2773                   fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2774                   fHistograms->FillHistogram("ESD_Background_Pt",  momentumVectorbackgroundCandidate.Pt());
2775                   fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2776                   fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2777                   fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2778                   fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2779                   fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2780                   fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2781                   fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2782                   fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2783                   fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2784
2785
2786                   fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2787
2788
2789                   if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2790                     fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2791                     fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2792                   }
2793
2794                   if(massBG>0.5 && massBG<0.6){
2795                     fHistograms->FillHistogram("ESD_Background_alfa_pt0506",momentumVectorbackgroundCandidate.Pt(),alfa);
2796                   }
2797                   if(massBG>0.3 && massBG<0.4){
2798                     fHistograms->FillHistogram("ESD_Background_alfa_pt0304",momentumVectorbackgroundCandidate.Pt(),alfa);
2799                   }
2800
2801                   // test
2802                   fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2803                   fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2804                   fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin),  momentumVectorbackgroundCandidate.Pt());
2805                   fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2806                   fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2807                   fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2808                   fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2809                   fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!
2810                   fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2811                   fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2812                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2813                   fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2814
2815                   if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2816                     fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2817                     fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2818                   }
2819                 }
2820
2821                 if(alfa<0.1){
2822                   fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2823                 }
2824                 //  }
2825               }
2826               delete backgroundCandidate;
2827             }
2828           }
2829         }
2830       }
2831     } // end else (means use #v0s as multiplicity)
2832   } // end no rotation
2833 }
2834
2835
2836 void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
2837   //ProcessGammasForGammaJetAnalysis
2838
2839   Double_t distIsoMin;
2840
2841   CreateListOfChargedParticles();
2842
2843
2844   //  for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
2845   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
2846     AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
2847     TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
2848
2849     if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
2850       distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
2851
2852       if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
2853         CalculateJetCone(gammaIndex);
2854       }
2855     }
2856   }
2857 }
2858
2859 //____________________________________________________________________
2860 Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(AliESDtrack *const track)
2861 {
2862 //
2863 // check whether particle has good DCAr(Pt) impact
2864 // parameter. Only for TPC+ITS tracks (7*sigma cut)
2865 // Origin: Andrea Dainese
2866 //
2867
2868 Float_t d0z0[2],covd0z0[3];
2869 track->GetImpactParameters(d0z0,covd0z0);
2870 Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
2871 Float_t d0max = 7.*sigma;
2872 if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
2873
2874 return kFALSE;
2875 }
2876
2877
2878 void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
2879   // CreateListOfChargedParticles
2880
2881   fESDEvent = fV0Reader->GetESDEvent();
2882   Int_t numberOfESDTracks=0;
2883   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2884     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
2885
2886     if(!curTrack){
2887       continue;
2888     }
2889     // Not needed if Standard function used.
2890 //     if(!IsGoodImpPar(curTrack)){
2891 //       continue;
2892 //     }
2893
2894     if(fEsdTrackCuts->AcceptTrack(curTrack) ){
2895       new((*fChargedParticles)[fChargedParticles->GetEntriesFast()])  AliESDtrack(*curTrack);
2896       //      fChargedParticles.push_back(curTrack);
2897       fChargedParticlesId.push_back(iTracks);
2898       numberOfESDTracks++;
2899     }
2900   }
2901   fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
2902
2903   if (fV0Reader->GetNumberOfContributorsVtx()>=1){
2904     fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
2905   }
2906 }
2907 void AliAnalysisTaskGammaConversion::RecalculateV0ForGamma(){
2908
2909  Double_t massE=0.00051099892;
2910  TLorentzVector curElecPos;
2911  TLorentzVector curElecNeg;
2912  TLorentzVector curGamma;
2913
2914  TLorentzVector curGammaAt;
2915  TLorentzVector curElecPosAt;
2916  TLorentzVector curElecNegAt;
2917  AliKFVertex primVtxGamma(*(fESDEvent->GetPrimaryVertex()));
2918  AliKFVertex primVtxImprovedGamma = primVtxGamma;
2919
2920  const AliESDVertex *vtxT3D=fESDEvent->GetPrimaryVertex();
2921
2922  Double_t xPrimaryVertex=vtxT3D->GetXv();
2923  Double_t yPrimaryVertex=vtxT3D->GetYv();
2924  Double_t zPrimaryVertex=vtxT3D->GetZv();
2925  // Float_t primvertex[3]={xPrimaryVertex,yPrimaryVertex,zPrimaryVertex};
2926
2927  Float_t nsigmaTPCtrackPos;
2928  Float_t nsigmaTPCtrackNeg;
2929  Float_t nsigmaTPCtrackPosToPion;
2930  Float_t nsigmaTPCtrackNegToPion;
2931  AliKFParticle* negKF=NULL;
2932  AliKFParticle* posKF=NULL;
2933
2934  for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2935    AliESDtrack* posTrack = fESDEvent->GetTrack(iTracks);
2936    if(!posTrack){
2937      continue;
2938    }
2939    if (posKF) delete posKF; posKF=NULL;
2940    if(posTrack->GetSign()<0) continue;
2941    if(!(posTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
2942    if(posTrack->GetKinkIndex(0)>0 ) continue;
2943    if(posTrack->GetNcls(1)<50)continue;
2944    Double_t momPos[3];
2945    //    posTrack->GetConstrainedPxPyPz(momPos);
2946    posTrack->GetPxPyPz(momPos);
2947    AliESDtrack *ptrk=fESDEvent->GetTrack(iTracks);
2948    curElecPos.SetXYZM(momPos[0],momPos[1],momPos[2],massE);
2949    if(TMath::Abs(curElecPos.Eta())<0.9) continue;
2950    posKF = new AliKFParticle( *(posTrack),-11);
2951
2952    nsigmaTPCtrackPos = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kElectron);
2953    nsigmaTPCtrackPosToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion);
2954
2955    if ( nsigmaTPCtrackPos>5.|| nsigmaTPCtrackPos<-2.){
2956      continue;
2957    }
2958
2959    if(pow((momPos[0]*momPos[0]+momPos[1]*momPos[1]+momPos[2]*momPos[2]),0.5)>0.5 && nsigmaTPCtrackPosToPion<1){
2960      continue;
2961    }
2962
2963
2964
2965    for(Int_t jTracks = 0; jTracks < fESDEvent->GetNumberOfTracks(); jTracks++){
2966      AliESDtrack* negTrack = fESDEvent->GetTrack(jTracks);
2967      if(!negTrack){
2968        continue;
2969      }
2970      if (negKF) delete negKF; negKF=NULL;
2971      if(negTrack->GetSign()>0) continue;
2972      if(!(negTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
2973      if(negTrack->GetKinkIndex(0)>0 ) continue;
2974      if(negTrack->GetNcls(1)<50)continue;
2975      Double_t momNeg[3];
2976      //    negTrack->GetConstrainedPxPyPz(momNeg);
2977      negTrack->GetPxPyPz(momNeg);
2978
2979      nsigmaTPCtrackNeg = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kElectron);
2980      nsigmaTPCtrackNegToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion);
2981      if ( nsigmaTPCtrackNeg>5. || nsigmaTPCtrackNeg<-2.){
2982        continue;
2983      }
2984      if(pow((momNeg[0]*momNeg[0]+momNeg[1]*momNeg[1]+momNeg[2]*momNeg[2]),0.5)>0.5 && nsigmaTPCtrackNegToPion<1){
2985        continue;
2986      }
2987      AliESDtrack *ntrk=fESDEvent->GetTrack(jTracks);
2988      curElecNeg.SetXYZM(momNeg[0],momNeg[1],momNeg[2],massE);
2989      if(TMath::Abs(curElecNeg.Eta())<0.9) continue;
2990      negKF = new AliKFParticle( *(negTrack) ,11);
2991
2992      Double_t b=fESDEvent->GetMagneticField();
2993      Double_t xn, xp, dca=ntrk->GetDCA(ptrk,b,xn,xp);
2994      AliExternalTrackParam nt(*ntrk), pt(*ptrk);
2995      nt.PropagateTo(xn,b); pt.PropagateTo(xp,b);
2996
2997
2998      //--- Like in ITSV0Finder
2999      AliExternalTrackParam ntAt0(*ntrk), ptAt0(*ptrk);
3000      Double_t xxP,yyP,alphaP;
3001      Double_t rP[3];
3002
3003      //     if (!ptAt0.GetGlobalXYZat(ptAt0->GetX(),xxP,yyP,zzP)) continue;
3004      if (!ptAt0.GetXYZAt(ptAt0.GetX(),b,rP)) continue;
3005      xxP=rP[0];
3006      yyP=rP[1];
3007      alphaP = TMath::ATan2(yyP,xxP);
3008
3009
3010      ptAt0.Propagate(alphaP,0,b);
3011      Float_t ptfacP  = (1.+100.*TMath::Abs(ptAt0.GetC(b)));
3012
3013      //     Double_t distP      = ptAt0.GetY();
3014      Double_t normP      = ptfacP*TMath::Sqrt(ptAt0.GetSigmaY2());
3015      Double_t normdist0P = TMath::Abs(ptAt0.GetY()/normP);
3016      Double_t normdist1P = TMath::Abs((ptAt0.GetZ()-zPrimaryVertex)/(ptfacP*TMath::Sqrt(ptAt0.GetSigmaZ2())));
3017      Double_t normdistP  = TMath::Sqrt(normdist0P*normdist0P+normdist1P*normdist1P);
3018
3019
3020      Double_t xxN,yyN,alphaN;
3021      Double_t rN[3];
3022      //     if (!ntAt0.GetGlobalXYZat(ntAt0->GetX(),xxN,yyN,zzN)) continue;
3023      if (!ntAt0.GetXYZAt(ntAt0.GetX(),b,rN)) continue;
3024      xxN=rN[0];
3025      yyN=rN[1];
3026
3027      alphaN = TMath::ATan2(yyN,xxN);
3028
3029      ntAt0.Propagate(alphaN,0,b);
3030
3031      Float_t ptfacN  = (1.+100.*TMath::Abs(ntAt0.GetC(b)));
3032      //     Double_t distN      = ntAt0.GetY();
3033      Double_t normN      = ptfacN*TMath::Sqrt(ntAt0.GetSigmaY2());
3034      Double_t normdist0N = TMath::Abs(ntAt0.GetY()/normN);
3035      Double_t normdist1N = TMath::Abs((ntAt0.GetZ()-zPrimaryVertex)/(ptfacN*TMath::Sqrt(ntAt0.GetSigmaZ2())));
3036      Double_t normdistN  = TMath::Sqrt(normdist0N*normdist0N+normdist1N*normdist1N);
3037
3038      //-----------------------------
3039
3040      Double_t momNegAt[3];
3041      nt.GetPxPyPz(momNegAt);
3042      curElecNegAt.SetXYZM(momNegAt[0],momNegAt[1],momNegAt[2],massE);
3043
3044      Double_t momPosAt[3];
3045      pt.GetPxPyPz(momPosAt);
3046      curElecPosAt.SetXYZM(momPosAt[0],momPosAt[1],momPosAt[2],massE);
3047      if(dca>1){
3048        continue;
3049      }
3050
3051      //     Double_t dneg= negTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
3052      //     Double_t dpos= posTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
3053      AliESDv0 vertex(nt,jTracks,pt,iTracks);
3054
3055
3056      Float_t cpa=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
3057
3058
3059
3060      //  cout<< "v0Rr::"<< v0Rr<<endl;
3061      // if (pvertex.GetRr()<0.5){
3062      // continue;
3063      //}
3064      //     cout<<"vertex.GetChi2V0()"<<vertex.GetChi2V0()<<endl;
3065      if(cpa<0.9)continue;
3066      //     if (vertex.GetChi2V0() > 30) continue;
3067      //     cout<<"xp+xn::"<<xp<<" "<<xn<<endl;
3068      if ((xn+xp) < 0.4) continue;
3069      if (TMath::Abs(ntrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05)
3070        if (TMath::Abs(ptrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05) continue;
3071
3072      //cout<<"pass"<<endl;
3073
3074      AliKFParticle v0GammaC;
3075      v0GammaC+=(*negKF);
3076      v0GammaC+=(*posKF);
3077      v0GammaC.SetMassConstraint(0,0.001);
3078      primVtxImprovedGamma+=v0GammaC;
3079      v0GammaC.SetProductionVertex(primVtxImprovedGamma);
3080
3081
3082      curGamma=curElecNeg+curElecPos;
3083      curGammaAt=curElecNegAt+curElecPosAt;
3084
3085      // invariant mass versus pt of K0short
3086
3087      Double_t chi2V0GammaC=100000.;
3088      if( v0GammaC.GetNDF() != 0) {
3089        chi2V0GammaC = v0GammaC.GetChi2()/v0GammaC.GetNDF();
3090      }else{
3091        cout<< "ERROR::v0K0C.GetNDF()" << endl;
3092      }
3093
3094      if(chi2V0GammaC<200 &&chi2V0GammaC>0 ){
3095        if(fHistograms != NULL){
3096          fHistograms->FillHistogram("ESD_RecalculateV0_InvMass",v0GammaC.GetMass());
3097          fHistograms->FillHistogram("ESD_RecalculateV0_Pt",v0GammaC.GetPt());
3098          fHistograms->FillHistogram("ESD_RecalculateV0_E_dEdxP",curElecNegAt.P(),negTrack->GetTPCsignal());
3099          fHistograms->FillHistogram("ESD_RecalculateV0_P_dEdxP",curElecPosAt.P(),posTrack->GetTPCsignal());
3100          fHistograms->FillHistogram("ESD_RecalculateV0_cpa",cpa);
3101          fHistograms->FillHistogram("ESD_RecalculateV0_dca",dca);
3102          fHistograms->FillHistogram("ESD_RecalculateV0_normdistP",normdistP);
3103          fHistograms->FillHistogram("ESD_RecalculateV0_normdistN",normdistN);
3104
3105          new((*fKFRecalculatedGammasTClone)[fKFRecalculatedGammasTClone->GetEntriesFast()])  AliKFParticle(v0GammaC);
3106          fElectronRecalculatedv1.push_back(iTracks);
3107          fElectronRecalculatedv2.push_back(jTracks);
3108        }
3109      }
3110    }
3111  }
3112
3113  for(Int_t firstGammaIndex=0;firstGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();firstGammaIndex++){
3114    for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();secondGammaIndex++){
3115       AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(firstGammaIndex);
3116       AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(secondGammaIndex);
3117
3118       if(fElectronRecalculatedv1[firstGammaIndex]==fElectronRecalculatedv1[secondGammaIndex]){
3119         continue;
3120       }
3121       if( fElectronRecalculatedv2[firstGammaIndex]==fElectronRecalculatedv2[secondGammaIndex]){
3122         continue;
3123       }
3124
3125       AliKFParticle twoGammaCandidate(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
3126       fHistograms->FillHistogram("ESD_RecalculateGG_InvMass",twoGammaCandidate.GetMass());
3127       fHistograms->FillHistogram("ESD_RecalculateGG_InvMass_vs_Pt",twoGammaCandidate.GetMass(),twoGammaCandidate.GetPt());
3128    }
3129  }
3130 }
3131 void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
3132   // CaculateJetCone
3133
3134   Double_t cone;
3135   Double_t coneSize=0.3;
3136   Double_t ptJet=0;
3137
3138   //  AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
3139   AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
3140
3141   TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
3142
3143   AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
3144
3145   Double_t momLeadingCharged[3];
3146   leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
3147
3148   TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
3149
3150   Double_t phi1=momentumVectorLeadingCharged.Phi();
3151   Double_t eta1=momentumVectorLeadingCharged.Eta();
3152   Double_t phi3=momentumVectorCurrentGamma.Phi();
3153
3154   for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
3155     AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
3156     Int_t chId = fChargedParticlesId[iCh];
3157     if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
3158     Double_t mom[3];
3159     curTrack->GetConstrainedPxPyPz(mom);
3160     TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
3161     Double_t phi2=momentumVectorChargedParticle.Phi();
3162     Double_t eta2=momentumVectorChargedParticle.Eta();
3163
3164
3165     cone=100.;
3166     if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
3167       cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
3168     }else{
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       if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
3173         cone = TMath::Sqrt(  TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
3174       }
3175     }
3176
3177     if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
3178       ptJet+= momentumVectorChargedParticle.Pt();
3179       Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
3180       Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
3181       fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
3182       fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
3183
3184     }
3185
3186     Double_t dphiHdrGam=phi3-phi2;
3187     if ( dphiHdrGam < (-TMath::PiOver2())){
3188       dphiHdrGam+=(TMath::TwoPi());
3189     }
3190
3191     if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
3192       dphiHdrGam-=(TMath::TwoPi());
3193     }
3194
3195     if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
3196       fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
3197     }
3198   }//track loop
3199
3200
3201 }
3202
3203 Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
3204   // GetMinimumDistanceToCharge
3205
3206   Double_t fIsoMin=100.;
3207   Double_t ptLeadingCharged=-1.;
3208
3209   fLeadingChargedIndex=-1;
3210
3211   AliKFParticle * gammaHighestPt = (AliKFParticle*)fKFReconstructedGammasTClone->At(indexHighestPtGamma);
3212   TVector3 momentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
3213
3214   Double_t phi1=momentumVectorgammaHighestPt.Phi();
3215   Double_t eta1=momentumVectorgammaHighestPt.Eta();
3216
3217   for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
3218     AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
3219     Int_t chId = fChargedParticlesId[iCh];
3220     if(fElectronv1[indexHighestPtGamma]==chId || fElectronv2[indexHighestPtGamma]==chId) continue;
3221     Double_t mom[3];
3222     curTrack->GetConstrainedPxPyPz(mom);
3223     TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
3224     Double_t phi2=momentumVectorChargedParticle.Phi();
3225     Double_t eta2=momentumVectorChargedParticle.Eta();
3226     Double_t iso=pow(  (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
3227
3228     if(momentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
3229       if (iso<fIsoMin){
3230         fIsoMin=iso;
3231       }
3232     }
3233
3234     Double_t dphiHdrGam=phi1-phi2;
3235     if ( dphiHdrGam < (-TMath::PiOver2())){
3236       dphiHdrGam+=(TMath::TwoPi());
3237     }
3238
3239     if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
3240       dphiHdrGam-=(TMath::TwoPi());
3241     }
3242     if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
3243       fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
3244     }
3245
3246     if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
3247       if (momentumVectorChargedParticle.Pt()> ptLeadingCharged && momentumVectorChargedParticle.Pt()>0.1*momentumVectorgammaHighestPt.Pt()){
3248         ptLeadingCharged=momentumVectorChargedParticle.Pt();
3249         fLeadingChargedIndex=iCh;
3250       }
3251     }
3252
3253   }//track loop
3254   fHistograms->FillHistogram("ESD_MinimumIsoDistance",fIsoMin);
3255   return fIsoMin;
3256
3257 }
3258
3259 Int_t  AliAnalysisTaskGammaConversion::GetIndexHighestPtGamma(){
3260   //GetIndexHighestPtGamma
3261
3262   Int_t indexHighestPtGamma=-1;
3263   //Double_t
3264   fGammaPtHighest = -100.;
3265
3266   for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
3267     AliKFParticle * gammaHighestPtCandidate = (AliKFParticle*)fKFReconstructedGammasTClone->At(firstGammaIndex);
3268     TVector3 momentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
3269     if (momentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
3270       fGammaPtHighest=momentumVectorgammaHighestPtCandidate.Pt();
3271       //gammaHighestPt = gammaHighestPtCandidate;
3272       indexHighestPtGamma=firstGammaIndex;
3273     }
3274   }
3275
3276   return indexHighestPtGamma;
3277
3278 }
3279
3280
3281 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
3282 {
3283   // Terminate analysis
3284   //
3285   AliDebug(1,"Do nothing in Terminate");
3286 }
3287
3288 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
3289 {
3290   //AOD
3291   if(!fAODGamma) fAODGamma = new TClonesArray("AliGammaConversionAODObject", 0);
3292   else fAODGamma->Delete();
3293   fAODGamma->SetName(Form("%s_gamma", fAODBranchName.Data()));
3294
3295   if(!fAODPi0) fAODPi0 = new TClonesArray("AliGammaConversionAODObject", 0);
3296   else fAODPi0->Delete();
3297   fAODPi0->SetName(Form("%s_Pi0", fAODBranchName.Data()));
3298
3299   if(!fAODOmega) fAODOmega = new TClonesArray("AliGammaConversionAODObject", 0);
3300   else fAODOmega->Delete();
3301   fAODOmega->SetName(Form("%s_Omega", fAODBranchName.Data()));
3302
3303   //If delta AOD file name set, add in separate file. Else add in standard aod file.
3304   if(fKFDeltaAODFileName.Length() > 0) {
3305     AddAODBranch("TClonesArray", &fAODGamma, fKFDeltaAODFileName.Data());
3306     AddAODBranch("TClonesArray", &fAODPi0, fKFDeltaAODFileName.Data());
3307     AddAODBranch("TClonesArray", &fAODOmega, fKFDeltaAODFileName.Data());
3308     AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fKFDeltaAODFileName.Data());
3309   } else  {
3310     AddAODBranch("TClonesArray", &fAODGamma);
3311     AddAODBranch("TClonesArray", &fAODPi0);
3312     AddAODBranch("TClonesArray", &fAODOmega);
3313   }
3314
3315   // Create the output container
3316   if(fOutputContainer != NULL){
3317     delete fOutputContainer;
3318     fOutputContainer = NULL;
3319   }
3320   if(fOutputContainer == NULL){
3321     fOutputContainer = new TList();
3322     fOutputContainer->SetOwner(kTRUE);
3323   }
3324
3325   //Adding the histograms to the output container
3326   fHistograms->GetOutputContainer(fOutputContainer);
3327
3328
3329   if(fWriteNtuple){
3330     if(fGammaNtuple == NULL){
3331       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);
3332     }
3333     if(fNeutralMesonNtuple == NULL){
3334       fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
3335     }
3336     TList * ntupleTList = new TList();
3337     ntupleTList->SetOwner(kTRUE);
3338     ntupleTList->SetName("Ntuple");
3339     ntupleTList->Add((TNtuple*)fGammaNtuple);
3340     fOutputContainer->Add(ntupleTList);
3341   }
3342
3343   fOutputContainer->SetName(GetName());
3344 }
3345
3346 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{
3347   //helper function
3348   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
3349   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
3350   return v3D0.Angle(v3D1);
3351 }
3352
3353 void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
3354   // see header file for documentation
3355
3356   vector<Int_t> indexOfGammaParticle;
3357
3358   fStack = fV0Reader->GetMCStack();
3359
3360   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
3361     return; // aborts if the primary vertex does not have contributors.
3362   }
3363
3364   for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
3365     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
3366     if(particle->GetPdgCode()==22){     //Gamma
3367       if(particle->GetNDaughters() >= 2){
3368         TParticle* electron=NULL;
3369         TParticle* positron=NULL;
3370         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
3371           TParticle *tmpDaughter = fStack->Particle(daughterIndex);
3372           if(tmpDaughter->GetUniqueID() == 5){
3373             if(tmpDaughter->GetPdgCode() == 11){
3374               electron = tmpDaughter;
3375             }
3376             else if(tmpDaughter->GetPdgCode() == -11){
3377               positron = tmpDaughter;
3378             }
3379           }
3380         }
3381         if(electron!=NULL && positron!=0){
3382           if(electron->R()<160){
3383             indexOfGammaParticle.push_back(iTracks);
3384           }
3385         }
3386       }
3387     }
3388   }
3389
3390   Int_t nFoundGammas=0;
3391   Int_t nNotFoundGammas=0;
3392
3393   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
3394   for(Int_t i=0;i<numberOfV0s;i++){
3395     fV0Reader->GetV0(i);
3396
3397     if(fV0Reader->HasSameMCMother() == kFALSE){
3398       continue;
3399     }
3400
3401     TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
3402     TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
3403
3404     if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
3405       continue;
3406     }
3407     if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
3408       continue;
3409     }
3410
3411     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
3412       //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
3413       for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
3414         if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
3415           nFoundGammas++;
3416         }
3417         else{
3418           nNotFoundGammas++;
3419         }
3420       }
3421     }
3422   }
3423 }
3424
3425
3426 void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
3427   // see header file for documantation
3428
3429   fESDEvent = fV0Reader->GetESDEvent();
3430
3431
3432   TClonesArray * vESDeNegTemp = new TClonesArray("AliESDtrack",0);
3433   TClonesArray * vESDePosTemp = new TClonesArray("AliESDtrack",0);
3434   TClonesArray * vESDxNegTemp = new TClonesArray("AliESDtrack",0);
3435   TClonesArray * vESDxPosTemp = new TClonesArray("AliESDtrack",0);
3436   TClonesArray * vESDeNegNoJPsi = new TClonesArray("AliESDtrack",0);
3437   TClonesArray * vESDePosNoJPsi = new TClonesArray("AliESDtrack",0);
3438
3439   /*
3440     vector <AliESDtrack*> vESDeNegTemp(0);
3441     vector <AliESDtrack*> vESDePosTemp(0);
3442     vector <AliESDtrack*> vESDxNegTemp(0);
3443     vector <AliESDtrack*> vESDxPosTemp(0);
3444     vector <AliESDtrack*> vESDeNegNoJPsi(0);
3445     vector <AliESDtrack*> vESDePosNoJPsi(0);
3446   */
3447
3448
3449   fHistograms->FillTable("Table_Electrons",0);//Count number of Events
3450
3451   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
3452     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
3453
3454     if(!curTrack){
3455       //print warning here
3456       continue;
3457     }
3458
3459     double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
3460     double r[3];curTrack->GetConstrainedXYZ(r);
3461
3462     TVector3 rXYZ(r);
3463
3464     fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
3465
3466     Bool_t flagKink       =  kTRUE;
3467     Bool_t flagTPCrefit   =  kTRUE;
3468     Bool_t flagTRDrefit   =  kTRUE;
3469     Bool_t flagITSrefit   =  kTRUE;
3470     Bool_t flagTRDout     =  kTRUE;
3471     Bool_t flagVertex     =  kTRUE;
3472
3473
3474     //Cuts ---------------------------------------------------------------
3475
3476     if(curTrack->GetKinkIndex(0) > 0){
3477       fHistograms->FillHistogram("Table_Electrons",5);//Count kink
3478       flagKink = kFALSE;
3479     }
3480
3481     ULong_t trkStatus = curTrack->GetStatus();
3482
3483     ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
3484
3485     if(!tpcRefit){
3486       fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
3487       flagTPCrefit = kFALSE;
3488     }
3489
3490     ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
3491     if(!itsRefit){
3492       fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
3493       flagITSrefit = kFALSE;
3494     }
3495
3496     ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
3497
3498     if(!trdRefit){
3499       fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
3500       flagTRDrefit = kFALSE;
3501     }
3502
3503     ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
3504
3505     if(!trdOut) {
3506       fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
3507       flagTRDout = kFALSE;
3508     }
3509
3510     double nsigmaToVxt = GetSigmaToVertex(curTrack);
3511
3512     if(nsigmaToVxt > 3){
3513       fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
3514       flagVertex = kFALSE;
3515     }
3516
3517     if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
3518     fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
3519
3520
3521     Stat_t pid, weight;
3522     GetPID(curTrack, pid, weight);
3523
3524     if(pid!=0){
3525       fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
3526     }
3527
3528     if(pid == 0){
3529       fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
3530     }
3531
3532
3533
3534
3535
3536
3537     TLorentzVector curElec;
3538     curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
3539
3540
3541     if(fDoMCTruth){
3542       Int_t labelMC = TMath::Abs(curTrack->GetLabel());
3543       TParticle* curParticle = fStack->Particle(labelMC);
3544       if(curTrack->GetSign() > 0){
3545         if( pid == 0){
3546           fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
3547           fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
3548         }
3549         else{
3550           fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
3551           fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
3552         }
3553       }
3554     }
3555
3556
3557     if(curTrack->GetSign() > 0){
3558
3559       //     vESDxPosTemp.push_back(curTrack);
3560       new((*vESDxPosTemp)[vESDxPosTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
3561
3562       if( pid == 0){
3563
3564         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
3565         fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
3566         //      fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
3567         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
3568         //      fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
3569         //      vESDePosTemp.push_back(curTrack);
3570         new((*vESDePosTemp)[vESDePosTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
3571
3572       }
3573
3574     }
3575     else {
3576
3577       new((*vESDxNegTemp)[vESDxNegTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
3578
3579       if( pid == 0){
3580
3581         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
3582         fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
3583         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
3584         new((*vESDeNegTemp)[vESDeNegTemp->GetEntriesFast()])  AliESDtrack(*curTrack);
3585
3586       }
3587
3588     }
3589
3590   }
3591
3592
3593   Bool_t ePosJPsi = kFALSE;
3594   Bool_t eNegJPsi = kFALSE;
3595   Bool_t ePosPi0  = kFALSE;
3596   Bool_t eNegPi0  = kFALSE;
3597
3598   UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
3599
3600   for(Int_t iNeg=0; iNeg < vESDeNegTemp->GetEntriesFast(); iNeg++){
3601     if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetPdgCode() == 11)
3602       if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0) > -1){
3603         Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0);
3604         TParticle* partMother = fStack ->Particle(labelMother);
3605         if (partMother->GetPdgCode() == 111){
3606           ieNegPi0 = iNeg;
3607           eNegPi0 = kTRUE;
3608         }
3609         if(partMother->GetPdgCode() == 443){ //Mother JPsi
3610           fHistograms->FillTable("Table_Electrons",14);
3611           ieNegJPsi = iNeg;
3612           eNegJPsi = kTRUE;
3613         }
3614         else{
3615           //      vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
3616           new((*vESDeNegNoJPsi)[vESDeNegNoJPsi->GetEntriesFast()])  AliESDtrack(*(AliESDtrack*)(vESDeNegTemp->At(iNeg)));
3617           //            cout<<"ESD No Positivo JPsi "<<endl;
3618         }
3619
3620       }
3621   }
3622
3623   for(Int_t iPos=0; iPos < vESDePosTemp->GetEntriesFast(); iPos++){
3624     if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetPdgCode() == -11)
3625       if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0) > -1){
3626         Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0);
3627         TParticle* partMother = fStack ->Particle(labelMother);
3628         if (partMother->GetPdgCode() == 111){
3629           iePosPi0 = iPos;
3630           ePosPi0 = kTRUE;
3631         }
3632         if(partMother->GetPdgCode() == 443){ //Mother JPsi
3633           fHistograms->FillTable("Table_Electrons",15);
3634           iePosJPsi = iPos;
3635           ePosJPsi = kTRUE;
3636         }
3637         else{
3638           //      vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
3639           new((*vESDePosNoJPsi)[vESDePosNoJPsi->GetEntriesFast()])  AliESDtrack(*(AliESDtrack*)(vESDePosTemp->At(iPos)));
3640           //            cout<<"ESD No Negativo JPsi "<<endl;
3641         }
3642
3643       }
3644   }
3645
3646   if( eNegJPsi && ePosJPsi ){
3647     TVector3 tempeNegV,tempePosV;
3648     tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Pz());
3649     tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Pz());
3650     fHistograms->FillTable("Table_Electrons",16);
3651     fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));
3652     fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->GetLabel())),
3653                                                                               fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->GetLabel()))));
3654   }
3655
3656   if( eNegPi0 && ePosPi0 ){
3657     TVector3 tempeNegV,tempePosV;
3658     tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Pz());
3659     tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Pz());
3660     fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
3661     fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->GetLabel())),
3662                                                                              fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->GetLabel()))));
3663   }
3664
3665
3666   FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
3667
3668   CleanWithAngleCuts(*vESDeNegTemp,*vESDePosTemp,*fKFReconstructedGammasTClone);
3669
3670   //  vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
3671   //  vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
3672
3673   TClonesArray vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectronTClone);
3674   TClonesArray vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectronTClone);
3675
3676
3677   FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
3678
3679
3680
3681
3682   //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
3683
3684
3685   FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
3686   FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
3687
3688
3689
3690   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",*fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
3691
3692   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
3693                            *fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
3694
3695   //BackGround
3696
3697   //Like Sign e+e-
3698   ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
3699   ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
3700   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
3701   ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
3702
3703   //        Like Sign e+e- no JPsi
3704   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
3705   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
3706
3707   //Mixed Event
3708
3709   if( fCurrentEventPosElectronTClone->GetEntriesFast() > 0 && fCurrentEventNegElectronTClone->GetEntriesFast() > 0 && fKFReconstructedGammasCutTClone->GetEntriesFast() > 0 ){
3710     FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
3711                              *fKFReconstructedGammasCutTClone,*fPreviousEventTLVNegElectronTClone,*fPreviousEventTLVPosElectronTClone);
3712     *fPreviousEventTLVNegElectronTClone = vCurrentTLVeNeg;
3713     *fPreviousEventTLVPosElectronTClone = vCurrentTLVePos;
3714
3715   }
3716
3717   /*
3718   //Photons P
3719   Double_t vtx[3];
3720   vtx[0]=0;vtx[1]=0;vtx[2]=0;
3721   for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
3722
3723   //      if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
3724
3725
3726
3727   Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
3728   //      cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
3729
3730   //      cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
3731
3732   if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
3733
3734   fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
3735
3736
3737   }
3738
3739
3740   */
3741
3742
3743   vESDeNegTemp->Delete();
3744   vESDePosTemp->Delete();
3745   vESDxNegTemp->Delete();
3746   vESDxPosTemp->Delete();
3747   vESDeNegNoJPsi->Delete();
3748   vESDePosNoJPsi->Delete();
3749
3750   delete vESDeNegTemp;
3751   delete vESDePosTemp;
3752   delete vESDxNegTemp;
3753   delete vESDxPosTemp;
3754   delete vESDeNegNoJPsi;
3755   delete vESDePosNoJPsi;
3756 }
3757
3758 /*
3759   void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
3760   //see header file for documentation
3761   for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
3762   for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
3763   fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
3764   }
3765   }
3766   }
3767 */
3768 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,TClonesArray const tlVeNeg, TClonesArray const tlVePos){
3769   //see header file for documentation
3770   for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++){
3771     for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
3772       fHistograms->FillHistogram(histoName.Data(),((TLorentzVector*)(tlVeNeg.At(iNeg)))->Vect().Angle(((TLorentzVector*)(tlVePos.At(iPos)))->Vect()));
3773     }
3774   }
3775 }
3776 void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, TClonesArray const eNeg, TClonesArray const ePos){
3777   //see header file for documentation
3778   for( Int_t n=0; n < eNeg.GetEntriesFast(); n++){
3779     TLorentzVector en = (*(TLorentzVector*)(eNeg.At(n)));
3780     for (Int_t p=0; p < ePos.GetEntriesFast(); p++){
3781       TLorentzVector ep = (*(TLorentzVector*)(ePos.At(p)));
3782       TLorentzVector np = ep + en;
3783       fHistograms->FillHistogram(histoName.Data(),np.M());
3784     }
3785   }
3786 }
3787
3788 void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,TClonesArray const fKFGammas,
3789                                                               TClonesArray const tlVeNeg,TClonesArray const tlVePos)
3790 {
3791   //see header file for documentation
3792
3793   for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++ ){
3794
3795     for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
3796
3797       TLorentzVector xy = *((TLorentzVector *)(tlVePos.At(iPos))) + *((TLorentzVector *)(tlVeNeg.At(iNeg)));
3798
3799       for (Int_t iGam=0; iGam < fKFGammas.GetEntriesFast(); iGam++){
3800
3801         //      AliKFParticle * gammaCandidate = &fKFGammas[iGam];
3802         AliKFParticle * gammaCandidate = (AliKFParticle *)(fKFGammas.At(iGam));
3803         TLorentzVector g;
3804
3805         g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
3806         TLorentzVector xyg = xy + g;
3807         fHistograms->FillHistogram(histoMass.Data(),xyg.M());
3808         fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
3809       }
3810     }
3811   }
3812
3813 }
3814 void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, TClonesArray e)
3815 {
3816   // see header file for documentation
3817   for(Int_t i=0; i < e.GetEntriesFast(); i++)
3818     {
3819       for (Int_t j=i+1; j < e.GetEntriesFast(); j++)
3820         {
3821           TLorentzVector ee = (*(TLorentzVector*)(e.At(i))) + (*(TLorentzVector*)(e.At(j)));
3822
3823           fHistograms->FillHistogram(hBg.Data(),ee.M());
3824         }
3825     }
3826 }
3827
3828
3829 void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negativeElectrons,
3830                                                         TClonesArray const positiveElectrons,
3831                                                         TClonesArray const gammas){
3832   // see header file for documentation
3833
3834   UInt_t  sizeN = negativeElectrons.GetEntriesFast();
3835   UInt_t  sizeP = positiveElectrons.GetEntriesFast();
3836   UInt_t  sizeG = gammas.GetEntriesFast();
3837
3838
3839
3840   vector <Bool_t> xNegBand(sizeN);
3841   vector <Bool_t> xPosBand(sizeP);
3842   vector <Bool_t> gammaBand(sizeG);
3843
3844
3845   for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
3846   for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
3847   for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
3848
3849
3850   for(UInt_t iPos=0; iPos < sizeP; iPos++){
3851
3852     Double_t aP[3];
3853     ((AliESDtrack*)(positiveElectrons.At(iPos)))->GetConstrainedPxPyPz(aP);
3854
3855     TVector3 ePosV(aP[0],aP[1],aP[2]);
3856
3857     for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
3858
3859       Double_t aN[3];
3860       ((AliESDtrack*)(negativeElectrons.At(iNeg)))->GetConstrainedPxPyPz(aN);
3861       TVector3 eNegV(aN[0],aN[1],aN[2]);
3862
3863       if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
3864         xPosBand[iPos]=kFALSE;
3865         xNegBand[iNeg]=kFALSE;
3866       }
3867
3868       for(UInt_t iGam=0; iGam < sizeG; iGam++){
3869         AliKFParticle* gammaCandidate = (AliKFParticle*)gammas.At(iGam);
3870         TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
3871         if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
3872           gammaBand[iGam]=kFALSE;
3873       }
3874     }
3875   }
3876
3877
3878
3879
3880   for(UInt_t iPos=0; iPos < sizeP; iPos++){
3881     if(xPosBand[iPos]){
3882       new((*fCurrentEventPosElectronTClone)[fCurrentEventPosElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(positiveElectrons.At(iPos))));
3883       //      fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
3884     }
3885   }
3886   for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
3887     if(xNegBand[iNeg]){
3888       new((*fCurrentEventNegElectronTClone)[fCurrentEventNegElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(negativeElectrons.At(iNeg))));
3889       //      fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
3890     }
3891   }
3892   for(UInt_t iGam=0; iGam < sizeG; iGam++){
3893     if(gammaBand[iGam]){
3894       new((*fKFReconstructedGammasCutTClone)[fKFReconstructedGammasCutTClone->GetEntriesFast()]) AliKFParticle((*(AliKFParticle*)(gammas.At(iGam))));
3895       //fKFReconstructedGammasCut.push_back(*(AliKFParticle*)gammas->At(iGam));
3896     }
3897   }
3898 }
3899
3900
3901 void  AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)
3902 {
3903   // see header file for documentation
3904   pid = -1;
3905   weight = -1;
3906
3907   double wpart[5];
3908   double wpartbayes[5];
3909
3910   //get probability of the diffenrent particle types
3911   track->GetESDpid(wpart);
3912
3913   // Tentative particle type "concentrations"
3914   double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
3915
3916   //Bayes' formula
3917   double rcc = 0.;
3918   for (int i = 0; i < 5; i++)
3919     {
3920       rcc+=(c[i] * wpart[i]);
3921     }
3922
3923
3924
3925   for (int i=0; i<5; i++) {
3926     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)
3927       wpartbayes[i] = c[i] * wpart[i] / rcc;
3928     }
3929   }
3930
3931
3932
3933   Float_t max=0.;
3934   int ipid=-1;
3935   //find most probable particle in ESD pid
3936   //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
3937   for (int i = 0; i < 5; i++)
3938     {
3939       if (wpartbayes[i] > max)
3940         {
3941           ipid = i;
3942           max = wpartbayes[i];
3943         }
3944     }
3945
3946   pid = ipid;
3947   weight = max;
3948 }
3949 double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)
3950 {
3951   // Calculates the number of sigma to the vertex.
3952
3953   Float_t b[2];
3954   Float_t bRes[2];
3955   Float_t bCov[3];
3956   t->GetImpactParameters(b,bCov);
3957   if (bCov[0]<=0 || bCov[2]<=0) {
3958     AliDebug(1, "Estimated b resolution lower or equal zero!");
3959     bCov[0]=0; bCov[2]=0;
3960   }
3961   bRes[0] = TMath::Sqrt(bCov[0]);
3962   bRes[1] = TMath::Sqrt(bCov[2]);
3963
3964   // -----------------------------------
3965   // How to get to a n-sigma cut?
3966   //
3967   // The accumulated statistics from 0 to d is
3968   //
3969   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
3970   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
3971   //
3972   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
3973   // Can this be expressed in a different way?
3974
3975   if (bRes[0] == 0 || bRes[1] ==0)
3976     return -1;
3977
3978   double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
3979
3980   // stupid rounding problem screws up everything:
3981   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
3982   if (TMath::Exp(-d * d / 2) < 1e-10)
3983     return 1000;
3984
3985
3986   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
3987   return d;
3988 }
3989
3990 //vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
3991 TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *const esdTrack){
3992   //Return TLoresntz vector of track?
3993   //  vector <TLorentzVector> tlVtrack(0);
3994   TClonesArray array("TLorentzVector",0);
3995
3996   for(Int_t itrack=0; itrack < esdTrack->GetEntriesFast(); itrack++){
3997     double p[3];
3998     //esdTrack[itrack]->GetConstrainedPxPyPz(p);
3999     ((AliESDtrack*)(esdTrack->At(itrack)))->GetConstrainedPxPyPz(p);
4000     TLorentzVector currentTrack;
4001     currentTrack.SetXYZM(p[0],p[1],p[2],fElectronMass);
4002     new((array)[array.GetEntriesFast()])  TLorentzVector(currentTrack);
4003     //    tlVtrack.push_back(currentTrack);
4004   }
4005
4006   return array;
4007
4008   //  return tlVtrack;
4009 }
4010
4011 Int_t AliAnalysisTaskGammaConversion::GetProcessType(AliMCEvent * mcEvt) {
4012
4013   // Determine if the event was generated with pythia or phojet and return the process type
4014
4015   // Check if mcEvt is fine
4016   if (!mcEvt) {
4017     AliFatal("NULL mc event");
4018   }
4019
4020   // Determine if it was a pythia or phojet header, and return the correct process type
4021   AliGenPythiaEventHeader * headPy  = 0;
4022   AliGenDPMjetEventHeader * headPho = 0;
4023   AliGenEventHeader * htmp = mcEvt->GenEventHeader();
4024   if(!htmp) {
4025     AliFatal("Cannot Get MC Header!!");
4026   }
4027   if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
4028     headPy =  (AliGenPythiaEventHeader*) htmp;
4029   } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
4030     headPho = (AliGenDPMjetEventHeader*) htmp;
4031   } else {
4032     AliError("Unknown header");
4033   }
4034
4035   // Determine process type
4036   if(headPy)   {
4037     if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
4038       // single difractive
4039       return kProcSD;
4040     } else if (headPy->ProcessType() == 94) {
4041       // double diffractive
4042       return kProcDD;
4043     }
4044     else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {
4045       // non difractive
4046       return kProcND;
4047     }
4048   } else if (headPho) {
4049     if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
4050       // single difractive
4051       return kProcSD;
4052     } else if (headPho->ProcessType() == 7) {
4053       // double diffractive
4054       return kProcDD;
4055     } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6  && headPho->ProcessType() != 7 ) {
4056       // non difractive
4057       return kProcND;
4058     }
4059   }
4060
4061
4062   // no process type found?
4063   AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
4064   return kProcUnknown;
4065 }