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