New more general analysis implemention for particle identification and correlation...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 19 May 2008 09:36:55 +0000 (09:36 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 19 May 2008 09:36:55 +0000 (09:36 +0000)
PWG4/AliAnaGammaDirect.cxx
PWG4/AliAnaGammaDirect.h
PWG4/AliAnaScale.cxx
PWG4/AliNeutralMesonSelection.cxx
PWG4/AliNeutralMesonSelection.h

index 2f0531e..cf53dab 100644 (file)
@@ -38,7 +38,7 @@
 //  Class created from old AliPHOSGammaJet 
 //  (see AliRoot versions previous Release 4-09)
 //
-//*-- Author: Gustavo Conesa (LNF-INFN) 
+// -- Author: Gustavo Conesa (LNF-INFN) 
 //////////////////////////////////////////////////////////////////////////////
   
   
 #include <TH2.h>
 #include <TList.h>
 #include "Riostream.h"
-#include "TROOT.h"
 
-// --- AliRoot system --- 
+// --- Analysis system --- 
 #include "AliAnaGammaDirect.h" 
 #include "AliLog.h"
+#include "AliCaloTrackReader.h"
+#include "AliIsolationCut.h"
+#include "AliNeutralMesonSelection.h"
 
 ClassImp(AliAnaGammaDirect)
   
 //____________________________________________________________________________
   AliAnaGammaDirect::AliAnaGammaDirect() : 
-    TObject(),
-    fMinGammaPt(0.),
-    fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
-    fICMethod(0), fAnaMC(0), fIsolatePi0(0),
-    fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0), 
-    fntuplePrompt(0),
-    //kSeveralIC
-    fNCones(0),fNPtThres(0), fConeSizes(),  fPtThresholds(), 
-    fhPtThresIsolated(), fhPtSumIsolated(), fntSeveralIC()
+    AliAnaBaseClass(), fDetector(""), fMakeIC(0),  fReMakeIC(0), 
+    fMakeSeveralIC(0), fMakeInvMass(0),
+    fhPtGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0),
+    //Several IC
+    fNCones(0),fNPtThresFrac(0), fConeSizes(),  fPtThresholds(),  fPtFractions(), 
+    fhPtThresIsolated(), fhPtFracIsolated(),  fhPtSumIsolated(),
+    //MC
+    fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0), 
+    fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(),  fhPtSumIsolatedPrompt(),
+    fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0), 
+    fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(),  fhPtSumIsolatedFragmentation(),
+    fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0), 
+    fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(),  fhPtSumIsolatedPi0Decay(),
+    fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0), 
+    fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(),  fhPtSumIsolatedOtherDecay(),
+    fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0), 
+    fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(),  fhPtSumIsolatedConversion(),
+    fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0), 
+    fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(),  fhPtSumIsolatedUnknown()
 {
   //default ctor
   
@@ -76,74 +88,147 @@ ClassImp(AliAnaGammaDirect)
 
 //____________________________________________________________________________
 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) : 
-  TObject(g),
-  fMinGammaPt(g.fMinGammaPt), 
-  fConeSize(g.fConeSize),
-  fPtThreshold(g.fPtThreshold),
-  fPtSumThreshold(g.fPtSumThreshold), 
-  fICMethod(g.fICMethod),
-  fAnaMC( g.fAnaMC), 
-  fIsolatePi0(g.fIsolatePi0),
-  fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),
-  fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),    
-  fntuplePrompt(g.fntuplePrompt),
-  //kSeveralIC
-  fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(), 
-  fhPtThresIsolated(), fhPtSumIsolated()
+  AliAnaBaseClass(g), fDetector(g.fDetector),
+  fMakeIC(g.fMakeIC),   fReMakeIC(g.fReMakeIC), 
+  fMakeSeveralIC(g.fMakeSeveralIC),  fMakeInvMass(g.fMakeInvMass),
+  fhPtGamma(g.fhPtGamma),fhPhiGamma(g.fhPhiGamma),
+  fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),  
+  //Several IC
+  fNCones(g.fNCones),fNPtThresFrac(g.fNPtThresFrac), fConeSizes(), fPtThresholds(),  fPtFractions(), 
+  fhPtThresIsolated(),  fhPtFracIsolated(), fhPtSumIsolated(),
+  //MC
+  fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt), 
+  fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(),  fhPtSumIsolatedPrompt(),
+  fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation), 
+  fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(),  fhPtSumIsolatedFragmentation(),
+  fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay), 
+  fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(),  fhPtSumIsolatedPi0Decay(),
+  fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay), 
+  fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(),  fhPtSumIsolatedOtherDecay(),
+  fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion), 
+  fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(),  fhPtSumIsolatedConversion(),
+  fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown), 
+  fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(),  fhPtSumIsolatedUnknown()
 {
   // cpy ctor
   
-  //kSeveralIC
+  //Several IC
   for(Int_t i = 0; i < fNCones ; i++){
     fConeSizes[i] =  g.fConeSizes[i];
     fhPtSumIsolated[i] = g.fhPtSumIsolated[i]; 
-    fntSeveralIC[i]= g.fntSeveralIC[i];  
-    for(Int_t j = 0; j < fNPtThres ; j++)
+
+   fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i]; 
+   fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i]; 
+   fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i]; 
+   fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i]; 
+   fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i]; 
+   fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i]; 
+
+    for(Int_t j = 0; j < fNPtThresFrac ; j++){
       fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j]; 
+      fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j];
+
+      fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j]; 
+      fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j]; 
+      fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j]; 
+      fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j]; 
+      fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j]; 
+      fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j]; 
+      fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j]; 
+      fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j]; 
+      fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j]; 
+      fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j]; 
+      fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j]; 
+      fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j]; 
+
+    } 
   }
   
-  for(Int_t i = 0; i < fNPtThres ; i++)
+  for(Int_t i = 0; i < fNPtThresFrac ; i++){
+    fPtFractions[i]=   g.fPtFractions[i];
     fPtThresholds[i]=   g.fPtThresholds[i];
-
+  }
 
 }
 
 //_________________________________________________________________________
-AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & source)
+AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & g)
 {
   // assignment operator
   
-  if(&source == this) return *this;
-
-  fMinGammaPt = source.fMinGammaPt ;   
-  fConeSize = source.fConeSize ;
-  fPtThreshold = source.fPtThreshold ;
-  fPtSumThreshold = source.fPtSumThreshold ; 
-  fICMethod = source.fICMethod ;
-  fAnaMC = source.fAnaMC ;
-  fIsolatePi0 =  source.fIsolatePi0 ;
-
-  fhNGamma = source.fhNGamma ; 
-  fhPhiGamma = source.fhPhiGamma ;
-  fhEtaGamma = source.fhEtaGamma ;
-  fhConeSumPt = source.fhConeSumPt ;
-
-  fntuplePrompt = source.fntuplePrompt ;
-
-  //kSeveralIC
-  fNCones = source.fNCones ;
-  fNPtThres = source.fNPtThres ; 
+  if(&g == this) return *this;
+
+  fMakeIC = g.fMakeIC ;
+  fReMakeIC = g.fReMakeIC ;
+  fMakeSeveralIC = g.fMakeSeveralIC ;
+  fMakeInvMass = g.fMakeInvMass ;
+  fDetector = g.fDetector ;
+
+  fhPtGamma = g.fhPtGamma ; 
+  fhPhiGamma = g.fhPhiGamma ;
+  fhEtaGamma = g.fhEtaGamma ;
+  fhConeSumPt = g.fhConeSumPt ;
+
+  fhPtPrompt = g.fhPtPrompt;
+  fhPhiPrompt = g.fhPhiPrompt;
+  fhEtaPrompt = g.fhEtaPrompt; 
+  fhPtFragmentation = g.fhPtFragmentation;
+  fhPhiFragmentation = g.fhPhiFragmentation;
+  fhEtaFragmentation = g.fhEtaFragmentation; 
+  fhPtPi0Decay = g.fhPtPi0Decay;
+  fhPhiPi0Decay = g.fhPhiPi0Decay;
+  fhEtaPi0Decay = g.fhEtaPi0Decay; 
+  fhPtOtherDecay = g.fhPtOtherDecay;
+  fhPhiOtherDecay = g.fhPhiOtherDecay;
+  fhEtaOtherDecay = g.fhEtaOtherDecay; 
+  fhPtConversion = g. fhPtConversion;
+  fhPhiConversion = g.fhPhiConversion;
+  fhEtaConversion = g.fhEtaConversion; 
+  fhPtUnknown = g.fhPtUnknown;
+  fhPhiUnknown = g.fhPhiUnknown;
+  fhEtaUnknown = g.fhEtaUnknown; 
+
+  //Several IC
+  fNCones = g.fNCones ;
+  fNPtThresFrac = g.fNPtThresFrac ; 
    
   for(Int_t i = 0; i < fNCones ; i++){
-    fConeSizes[i] =  source.fConeSizes[i];
-    fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
-    fntSeveralIC[i]= source.fntSeveralIC[i];  
-    for(Int_t j = 0; j < fNPtThres ; j++)
-      fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
+    fConeSizes[i] =  g.fConeSizes[i];
+    fhPtSumIsolated[i] = g.fhPtSumIsolated[i] ;
+
+   fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i]; 
+   fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i]; 
+   fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i]; 
+   fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i]; 
+   fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i]; 
+   fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i]; 
+
+    for(Int_t j = 0; j < fNPtThresFrac ; j++){
+      fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j] ;
+      fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j] ;
+
+      fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j]; 
+      fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j]; 
+      fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j]; 
+      fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j]; 
+      fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j]; 
+      fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j]; 
+      fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j]; 
+      fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j]; 
+      fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j]; 
+      fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j]; 
+      fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j]; 
+      fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j]; 
+
+    }
   }
   
-  for(Int_t i = 0; i < fNPtThres ; i++)
-    fPtThresholds[i]=   source.fPtThresholds[i];
+  for(Int_t i = 0; i < fNPtThresFrac ; i++){
+    fPtThresholds[i]=   g.fPtThresholds[i];
+    fPtFractions[i]=   g.fPtFractions[i];
+  }
 
   return *this;
   
@@ -152,302 +237,705 @@ AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & sou
 //____________________________________________________________________________
 AliAnaGammaDirect::~AliAnaGammaDirect() 
 {
-  // Remove all pointers except analysis output pointers.
+  //dtor
+  //do not delete histograms
+
+ delete [] fConeSizes ; 
+ delete [] fPtThresholds ; 
+ delete [] fPtFractions ; 
+
+}
+//_________________________________________________________________________
+Bool_t AliAnaGammaDirect::CheckInvMass(const Int_t icalo,const TLorentzVector mom, Double_t *vertex, TClonesArray * pl){
+  //Search if there is a companion decay photon to the candidate 
+  // and discard it in such case
+  TLorentzVector mom2 ;
+  for(Int_t jcalo = 0; jcalo < pl->GetEntries(); jcalo++){
+    if(icalo==jcalo) continue ;
+    AliAODCaloCluster * calo =  dynamic_cast<AliAODCaloCluster*> (pl->At(jcalo));
+  
+    //Cluster selection, not charged, with photon id and in fidutial cut, fill TLorentz
+    if(!SelectCluster(calo, vertex, mom2)) continue ;
+
+    //Select good pair (good phit, pt cuts, aperture and invariant mass)
+    if(GetNeutralMesonSelection()->SelectPair(mom, mom2)){
+      if(GetDebug()>1)printf("Selected gamma pair: pt %f, phi %f, eta%f",(mom+mom2).Pt(), (mom+mom2).Phi(), (mom+mom2).Eta());
+      return kTRUE ;
+    }
+  }//loop
+
+  return kFALSE;
 
+}
+//_________________________________________________________________________
+Int_t AliAnaGammaDirect::CheckOrigin(const Int_t label){
+  //Play with the MC stack if available
+  //Check origin of the candidates, good for PYTHIA
+  
+  AliStack * stack =  GetMCStack() ;
+  if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
+  
+  if(label >= 0 && label <  stack->GetNtrack()){
+    //Mother
+    TParticle * mom = stack->Particle(label);
+    Int_t mPdg = TMath::Abs(mom->GetPdgCode());
+    Int_t mStatus =  mom->GetStatusCode() ;
+    Int_t iParent =  mom->GetFirstMother() ;
+    if(label < 8 ) printf("Mother is parton %d\n",iParent);
+    
+    //GrandParent
+    TParticle * parent = new TParticle ;
+    Int_t pPdg = -1;
+    Int_t pStatus =-1;
+    if(iParent > 0){
+      parent = stack->Particle(iParent);
+      pPdg = TMath::Abs(parent->GetPdgCode());
+      pStatus = parent->GetStatusCode();  
+    }
+    else
+      printf("Parent with label %d\n",iParent);
+    
+    //return tag
+    if(mPdg == 22){
+      if(mStatus == 1){
+       if(iParent < 8) {
+         if(pPdg == 22) return kPrompt;
+         else  return kFragmentation;
+       }
+       else if(pStatus == 11){
+         if(pPdg == 111) return kPi0Decay ;
+         else if (pPdg == 321)  return kEtaDecay ;
+         else  return kOtherDecay ;
+       }
+      }//Status 1 : Pythia generated
+      else if(mStatus == 0){
+       if(pPdg ==22 || pPdg ==11) return kConversion ;
+       if(pPdg == 111) return kPi0Decay ;
+       else if (pPdg == 221)  return kEtaDecay ;
+       else  return kOtherDecay ;
+      }//status 0 : geant generated
+    }//Mother gamma
+    else if(mPdg == 111)  return kPi0 ;
+    else if(mPdg == 221)  return kEta ;
+    else if(mPdg ==11){
+      if(mStatus == 0) return kConversion ;
+      else return kElectron ;
+    }
+    else return kUnknown;
+  }//Good label value
+  else{
+    if(label < 0 ) printf("*** bad label or no stack ***:  label %d \n", label);
+    if(label >=  stack->GetNtrack()) printf("*** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
+    return kUnknown;
+  }//Bad label
+  
+  return kUnknown;
+  
 }
 
 //________________________________________________________________________
 TList *  AliAnaGammaDirect::GetCreateOutputObjects()
 {  
-
   // Create histograms to be saved in output file and 
   // store them in outputContainer
   TList * outputContainer = new TList() ; 
   outputContainer->SetName("DirectGammaHistos") ; 
-  
+
   //Histograms of highest gamma identified in Event
-  fhNGamma  = new TH1F("NGamma","Number of #gamma over calorimeter",240,0,120); 
-  fhNGamma->SetYTitle("N");
-  fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
-  outputContainer->Add(fhNGamma) ; 
+  fhPtGamma  = new TH1F("hPtGamma","Number of #gamma over calorimeter",240,0,120); 
+  fhPtGamma->SetYTitle("N");
+  fhPtGamma->SetXTitle("p_{T #gamma}(GeV/c)");
+  outputContainer->Add(fhPtGamma) ; 
   
   fhPhiGamma  = new TH2F
-    ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7); 
+    ("hPhiGamma","#phi_{#gamma}",200,0,120,200,0,7); 
   fhPhiGamma->SetYTitle("#phi");
   fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
   outputContainer->Add(fhPhiGamma) ; 
   
   fhEtaGamma  = new TH2F
-    ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
+    ("hEtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
   fhEtaGamma->SetYTitle("#eta");
   fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
   outputContainer->Add(fhEtaGamma) ;
 
-  fhConeSumPt  = new TH2F
-    ("ConePtSum","#Sigma p_{T}  in cone ",200,0,120,100,0,100);
-  fhConeSumPt->SetYTitle("#Sigma p_{T}");
-  fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
-  outputContainer->Add(fhConeSumPt) ;
+  if(!fMakeSeveralIC){
+    fhConeSumPt  = new TH2F
+      ("hConePtSum","#Sigma p_{T}  in cone ",200,0,120,100,0,100);
+    fhConeSumPt->SetYTitle("#Sigma p_{T}");
+    fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhConeSumPt) ;
+  }
   
-  //NTUPLE
-  fntuplePrompt = new TNtuple("ntuplePromptGamma", "Tree of prompt #gamma", "ptcluster:phicluster:etacluster:pdg:status:ptprimary:phiprimary:etaprimary:pdgprimary:statusprimary");
-  outputContainer->Add(fntuplePrompt) ;
-
-  if(fICMethod == kSeveralIC){
+  if(IsDataMC()){
+    
+    fhPtPrompt  = new TH1F("hPtPrompt","Number of #gamma over calorimeter",240,0,120); 
+    fhPtPrompt->SetYTitle("N");
+    fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtPrompt) ; 
+    
+    fhPhiPrompt  = new TH2F
+      ("hPhiPrompt","#phi_{#gamma}",200,0,120,200,0,7); 
+    fhPhiPrompt->SetYTitle("#phi");
+    fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPhiPrompt) ; 
+    
+    fhEtaPrompt  = new TH2F
+      ("hEtaPrompt","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
+    fhEtaPrompt->SetYTitle("#eta");
+    fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhEtaPrompt) ;
+    
+    fhPtFragmentation  = new TH1F("hPtFragmentation","Number of #gamma over calorimeter",240,0,120); 
+    fhPtFragmentation->SetYTitle("N");
+    fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtFragmentation) ; 
+    
+    fhPhiFragmentation  = new TH2F
+      ("hPhiFragmentation","#phi_{#gamma}",200,0,120,200,0,7); 
+    fhPhiFragmentation->SetYTitle("#phi");
+    fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPhiFragmentation) ; 
+    
+    fhEtaFragmentation  = new TH2F
+      ("hEtaFragmentation","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
+    fhEtaFragmentation->SetYTitle("#eta");
+    fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhEtaFragmentation) ;
+    
+    fhPtPi0Decay  = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",240,0,120); 
+    fhPtPi0Decay->SetYTitle("N");
+    fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtPi0Decay) ; 
+    
+    fhPhiPi0Decay  = new TH2F
+      ("hPhiPi0Decay","#phi_{#gamma}",200,0,120,200,0,7); 
+    fhPhiPi0Decay->SetYTitle("#phi");
+    fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPhiPi0Decay) ; 
+    
+    fhEtaPi0Decay  = new TH2F
+      ("hEtaPi0Decay","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
+    fhEtaPi0Decay->SetYTitle("#eta");
+    fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhEtaPi0Decay) ;
+    
+    fhPtOtherDecay  = new TH1F("hPtOtherDecay","Number of #gamma over calorimeter",240,0,120); 
+    fhPtOtherDecay->SetYTitle("N");
+    fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtOtherDecay) ; 
+    
+    fhPhiOtherDecay  = new TH2F
+      ("hPhiOtherDecay","#phi_{#gamma}",200,0,120,200,0,7); 
+    fhPhiOtherDecay->SetYTitle("#phi");
+    fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPhiOtherDecay) ; 
+    
+    fhEtaOtherDecay  = new TH2F
+      ("hEtaOtherDecay","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
+    fhEtaOtherDecay->SetYTitle("#eta");
+    fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhEtaOtherDecay) ;
+    
+    fhPtConversion  = new TH1F("hPtConversion","Number of #gamma over calorimeter",240,0,120); 
+    fhPtConversion->SetYTitle("N");
+    fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtConversion) ; 
+    
+    fhPhiConversion  = new TH2F
+      ("hPhiConversion","#phi_{#gamma}",200,0,120,200,0,7); 
+    fhPhiConversion->SetYTitle("#phi");
+    fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPhiConversion) ; 
+    
+    fhEtaConversion  = new TH2F
+      ("hEtaConversion","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
+    fhEtaConversion->SetYTitle("#eta");
+    fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhEtaConversion) ;
+    
+    fhPtUnknown  = new TH1F("hPtUnknown","Number of #gamma over calorimeter",240,0,120); 
+    fhPtUnknown->SetYTitle("N");
+    fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtUnknown) ; 
+    
+    fhPhiUnknown  = new TH2F
+      ("hPhiUnknown","#phi_{#gamma}",200,0,120,200,0,7); 
+    fhPhiUnknown->SetYTitle("#phi");
+    fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPhiUnknown) ; 
+    
+    fhEtaUnknown  = new TH2F
+      ("hEtaUnknown","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
+    fhEtaUnknown->SetYTitle("#eta");
+    fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhEtaUnknown) ;
+  }//Histos with MC
+  
+  if(fMakeSeveralIC){
     char name[128];
     char title[128];
     for(Int_t icone = 0; icone<fNCones; icone++){
-      sprintf(name,"PtSumIsolated_Cone_%d",icone);
+      sprintf(name,"hPtSumIsolated_Cone_%d",icone);
       sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
       fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
       fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
       fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
       outputContainer->Add(fhPtSumIsolated[icone]) ; 
-      
-      sprintf(name,"nt_Cone_%d",icone);
-      sprintf(title,"ntuple for cone size %d",icone);
-      fntSeveralIC[icone] = new TNtuple(name, title, "ptcand:phicand:etacand:pdg:status:ptsum:ncone0:ncone1:ncone2:ncone3:ncone4:ncone5");
-      outputContainer->Add(fntSeveralIC[icone]) ; 
-
-      for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
-       sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
+    
+      if(IsDataMC()){
+       sprintf(name,"hPtSumIsolatedPrompt_Cone_%d",icone);
+       sprintf(title,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+       fhPtSumIsolatedPrompt[icone]  = new TH2F(name, title,240,0,120,120,0,10);
+       fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+       fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
+       outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ; 
+
+       sprintf(name,"hPtSumIsolatedFragmentation_Cone_%d",icone);
+       sprintf(title,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+       fhPtSumIsolatedFragmentation[icone]  = new TH2F(name, title,240,0,120,120,0,10);
+       fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+       fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
+       outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ; 
+
+       sprintf(name,"hPtSumIsolatedPi0Decay_Cone_%d",icone);
+       sprintf(title,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+       fhPtSumIsolatedPi0Decay[icone]  = new TH2F(name, title,240,0,120,120,0,10);
+       fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+       fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
+       outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ; 
+
+       sprintf(name,"hPtSumIsolatedOtherDecay_Cone_%d",icone);
+       sprintf(title,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+       fhPtSumIsolatedOtherDecay[icone]  = new TH2F(name, title,240,0,120,120,0,10);
+       fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+       fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
+       outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ; 
+
+       sprintf(name,"hPtSumIsolatedConversion_Cone_%d",icone);
+       sprintf(title,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+       fhPtSumIsolatedConversion[icone]  = new TH2F(name, title,240,0,120,120,0,10);
+       fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+       fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
+       outputContainer->Add(fhPtSumIsolatedConversion[icone]) ; 
+
+       sprintf(name,"hPtSumIsolatedUnknown_Cone_%d",icone);
+       sprintf(title,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+       fhPtSumIsolatedUnknown[icone]  = new TH2F(name, title,240,0,120,120,0,10);
+       fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+       fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");
+       outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ; 
+
+      }//Histos with MC
+
+      for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){ 
+       sprintf(name,"hPtThresIsol_Cone_%d_Pt%d",icone,ipt);
        sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
        fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
        fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
        outputContainer->Add(fhPtThresIsolated[icone][ipt]) ; 
+
+       sprintf(name,"hPtFracIsol_Cone_%d_Pt%d",icone,ipt);
+       sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+       fhPtFracIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
+       fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+       outputContainer->Add(fhPtFracIsolated[icone][ipt]) ; 
+
+       if(IsDataMC()){
+         sprintf(name,"hPtThresIsolPrompt_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtThresIsolatedPrompt[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ; 
+         
+         sprintf(name,"hPtFracIsolPrompt_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtFracIsolatedPrompt[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ; 
+
+         sprintf(name,"hPtThresIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtThresIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ; 
+         
+         sprintf(name,"hPtFracIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtFracIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ; 
+
+         sprintf(name,"hPtThresIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtThresIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ; 
+         
+         sprintf(name,"hPtFracIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtFracIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ; 
+
+         sprintf(name,"hPtThresIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtThresIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ; 
+         
+         sprintf(name,"hPtFracIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtFracIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
+
+         sprintf(name,"hPtThresIsolConversion_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtThresIsolatedConversion[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ; 
+         
+         sprintf(name,"hPtFracIsolConversion_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtFracIsolatedConversion[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
+
+         sprintf(name,"hPtThresIsolUnknown_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtThresIsolatedUnknown[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ; 
+         
+         sprintf(name,"hPtFracIsolUnknown_Cone_%d_Pt%d",icone,ipt);
+         sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+         fhPtFracIsolatedUnknown[icone][ipt]  = new TH1F(name, title,240,0,120);
+         fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+         outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;  
+       }//Histos with MC
+
       }//icone loop
     }//ipt loop
   }
 
-//  gROOT->cd();
+  //Keep neutral meson selection histograms if requiered
+  //Setting done in AliNeutralMesonSelection
+  if(fMakeInvMass && GetNeutralMesonSelection()){
+    TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
+    cout<<"NMSHistos "<< nmsHistos<<endl;
+    if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
+      for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
+  }
 
   return outputContainer ;
 
 }
 
-//____________________________________________________________________________
-void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plCTS, TClonesArray * plPrimCalo, TParticle *pGamma, Bool_t &found) const 
+//__________________________________________________________________
+void  AliAnaGammaDirect::MakeAnalysisFillAOD() 
 {
-  //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
+  //Do analysis and fill aods
+  //Search for the isolated photon in fDetector with pt > GetMinPt()
 
-  Double_t pt = 0;
-  Int_t n = 0;
-  Int_t index = -1; 
-  Float_t coneptsum = 0 ;
+  //Fill CaloClusters if working with ESDs
+  //if(GetReader()->GetDataType() == AliCaloTrackReader::kESD) ConnectAODCaloClusters(); 
 
-  for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
-    TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
-
-    if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt) && 
-       (particle->GetPdgCode() == 22 || (fIsolatePi0 && particle->GetPdgCode() == 111))){
-      index = ipr ;
-      pt = particle->Pt();
-      pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
-      pGamma->SetStatusCode(particle->GetStatusCode());
-      pGamma->SetPdgCode(particle->GetPdgCode());
-      found  = kTRUE;
-    }
-  }
-  //Do Isolation?
-  if( ( fICMethod == kPtIC  ||  fICMethod == kSumPtIC )  && found)
-    {
-      Bool_t  icPtThres   = kFALSE;
-      Bool_t  icPtSum     = kFALSE;
-      MakeIsolationCut(plCTS,plCalo, pGamma, index,n,
-                      icPtThres, icPtSum,coneptsum);
-      if(fICMethod == kPtIC) //Pt thres method
-       found = icPtThres ;
-      if(fICMethod == kSumPtIC) //Pt cone sum method
-       found = icPtSum ;
-    }
+  Int_t n = 0, nfrac = 0;
+  Bool_t isolated = kFALSE ; 
+  Float_t coneptsum = 0 ;
+  TClonesArray * pl = new TClonesArray; 
+
+  //Get vertex for photon momentum calculation
+  Double_t vertex[]={0,0,0} ; //vertex ;
+  if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
+
+  //Select the detector of the photon
+  if(fDetector == "PHOS")
+    pl = GetAODPHOS();
+  else if (fDetector == "EMCAL")
+    pl = GetAODEMCAL();
+  //cout<<"Number of entries "<<pl->GetEntries()<<endl;
   
-  if(found){
-    AliDebug(1,Form("Cluster with p_{T} larger than the pt treshold %f found in calorimeter ", fMinGammaPt)) ;
-    AliDebug(1,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
-    //Fill prompt gamma histograms 
-    Float_t ptcluster = pGamma->Pt();
-    Float_t phicluster = pGamma->Phi();
-    Float_t etacluster = pGamma->Eta();
-    Int_t statuscluster = pGamma->GetStatusCode();
-    Int_t pdgcluster = pGamma->GetPdgCode();
-
-    fhNGamma   ->Fill(ptcluster);
-    fhPhiGamma ->Fill(ptcluster,phicluster);
-    fhEtaGamma ->Fill(ptcluster,etacluster);
-    fhConeSumPt->Fill(ptcluster,coneptsum);
-
-    Float_t ptprimary = 0 ;
-    Float_t phiprimary = 0 ;
-    Float_t etaprimary = 0 ;
-    Int_t pdgprimary = 0 ;
-    Int_t statusprimary = 0 ;
-
-    if(fAnaMC){
-      TParticle * primary = dynamic_cast<TParticle *>(plPrimCalo->At(index)) ;
-      ptprimary = primary->Pt();
-      phiprimary = primary->Phi();
-      etaprimary = primary->Eta();
-      pdgprimary =  TMath::Abs(primary->GetPdgCode()) ;
-      statusprimary = primary->GetStatusCode(); // = 2 means decay gamma!!!
-      AliDebug(1, Form("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f",ptcluster,pdgprimary,ptprimary));
-      //printf("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f \n",ptcluster,pdgprimary,ptprimary);
+  //Fill AODCaloClusters and AODParticleCorrelation with PHOS aods
+  TLorentzVector mom ;
+  for(Int_t icalo = 0; icalo < pl->GetEntries(); icalo++){
+    AliAODCaloCluster * calo =  dynamic_cast<AliAODCaloCluster*> (pl->At(icalo));
+  
+    //Cluster selection, not charged, with photon id and in fidutial cut
+    if(!SelectCluster(calo,vertex,mom)) continue ;
+    
+    //If too small pt, skip
+    if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; 
+
+    //Play with the MC stack if available
+    Int_t tag =-1;
+    if(IsDataMC()){
+      //Check origin of the candidates
+      tag = CheckOrigin(calo->GetLabel(0));
+      if(GetDebug() > 0) printf("Origin of candidate %d\n",tag);
+    }//Work with stack also   
+
+    //Check invariant mass
+    Bool_t decay = kFALSE ;
+    if(fMakeInvMass) decay = CheckInvMass(icalo,mom,vertex,pl);
+    if(decay) continue ;
+
+    //Create AOD for analysis
+    AliAODParticleCorrelation ph = AliAODParticleCorrelation(mom);
+    ph.SetLabel(calo->GetLabel(0));
+    ph.SetPdg(AliCaloPID::kPhoton);
+    ph.SetDetector(fDetector);
+    ph.SetTag(tag);  
+    if(fMakeIC){
+      n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
+      GetIsolationCut()->MakeIsolationCut((TSeqCollection*)GetAODCTS(), (TSeqCollection*)pl, 
+                                         vertex, kTRUE, &ph,icalo,-1,
+                                         n,nfrac,coneptsum, isolated);
+      if(isolated){
+       //cout<<"Isolated : E "<<mom.E()<<" ; Phi"<<mom.Phi()<< " ; Eta "<<mom.Eta()<<endl;
+       AddAODParticleCorrelation(ph);
+      }
     }
+    else //Fill all if not isolation requested
+      AddAODParticleCorrelation(ph);
 
-    //Fill ntuple with cluster / MC data
-//    gROOT->cd();
-    fntuplePrompt->Fill(ptcluster,phicluster,etacluster,pdgcluster,statuscluster,ptprimary,phiprimary, etaprimary,pdgprimary,statusprimary);
-  }
-  else
-    AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
+  }//loop
+  
+  if(GetDebug() > 1) printf("End fill AODs ");  
+  
 }
 
-  //____________________________________________________________________________
+//__________________________________________________________________
+void  AliAnaGammaDirect::MakeAnalysisFillHistograms() 
+{
+  //Do analysis and fill histograms
+  Int_t n = 0, nfrac = 0;
+  Bool_t isolated = kFALSE ; 
+  Float_t coneptsum = 0 ;
+  //cout<<"MakeAnalysisFillHistograms"<<endl;
+
+  //Get vertex for photon momentum calculation
+  Double_t v[]={0,0,0} ; //vertex ;
+  if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(v);
+
+  //Loop on stored AOD photons
+  Int_t naod = GetAODBranch()->GetEntriesFast();
+  if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
+  for(Int_t iaod = 0; iaod < naod ; iaod++){
+   AliAODParticleCorrelation* ph =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
+   Int_t pdg = ph->GetPdg();
+
+   //Only isolate photons in detector fDetector
+   if(pdg != AliCaloPID::kPhoton || ph->GetDetector() != fDetector) continue;
+   
+   if(fMakeSeveralIC) {
+     //Analysis of multiple IC at same time
+     MakeSeveralICAnalysis(ph,v);
+     continue;
+   }
+   else if(fReMakeIC){
+     //In case a more strict IC is needed in the produced AOD
+     n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
+     GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(), 
+                                        (TSeqCollection*)ph->GetRefIsolationConeClusters(), 
+                                        v, kFALSE, ph,-1,-1,
+                                        n,nfrac,coneptsum, isolated);
+   }
+
+   //Fill histograms (normal case)
+   if(!fReMakeIC || (fReMakeIC && isolated)){
+     
+     //printf("Isolated Gamma: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta()) ;
+     
+     //Fill prompt gamma histograms 
+     Float_t ptcluster = ph->Pt();
+     Float_t phicluster = ph->Phi();
+     Float_t etacluster = ph->Eta();
+     
+     fhPtGamma   ->Fill(ptcluster);
+     fhPhiGamma ->Fill(ptcluster,phicluster);
+     fhEtaGamma ->Fill(ptcluster,etacluster);
+     fhConeSumPt->Fill(ptcluster,coneptsum);
+
+     if(IsDataMC()){
+       Int_t tag =ph->GetTag();
+           
+       if(tag == kPrompt){
+        fhPtPrompt   ->Fill(ptcluster);
+        fhPhiPrompt ->Fill(ptcluster,phicluster);
+        fhEtaPrompt ->Fill(ptcluster,etacluster);
+       }
+       else if(tag==kFragmentation)
+        {
+          fhPtFragmentation   ->Fill(ptcluster);
+          fhPhiFragmentation ->Fill(ptcluster,phicluster);
+          fhEtaFragmentation ->Fill(ptcluster,etacluster);
+        }
+       else if(tag==kPi0Decay)
+        {
+          fhPtPi0Decay   ->Fill(ptcluster);
+          fhPhiPi0Decay ->Fill(ptcluster,phicluster);
+          fhEtaPi0Decay ->Fill(ptcluster,etacluster);
+        }
+       else if(tag==kEtaDecay || tag==kOtherDecay)
+        {
+          fhPtOtherDecay   ->Fill(ptcluster);
+          fhPhiOtherDecay ->Fill(ptcluster,phicluster);
+          fhEtaOtherDecay ->Fill(ptcluster,etacluster);
+        }
+       else if(tag==kConversion)
+        {
+          fhPtConversion   ->Fill(ptcluster);
+          fhPhiConversion ->Fill(ptcluster,phicluster);
+          fhEtaConversion ->Fill(ptcluster,etacluster);
+        }
+       else{
+        
+        fhPtUnknown   ->Fill(ptcluster);
+        fhPhiUnknown ->Fill(ptcluster,phicluster);
+        fhEtaUnknown ->Fill(ptcluster,etacluster);
+       }
+     }//Histograms with MC
+     
+   }
+  }// aod loop
+  
+}
+
+//____________________________________________________________________________
 void AliAnaGammaDirect::InitParameters()
 {
+  
   //Initialize the parameters of the analysis.
-  fMinGammaPt  = 5. ;
 
-  //Fill particle lists when PID is ok
-  fConeSize             = 0.4 ; 
-  fPtThreshold         = 1.; 
-  fPtSumThreshold  = 1.; 
+  fDetector = "PHOS" ;
+  fMakeIC = kTRUE;
+  fReMakeIC = kFALSE ;
+  fMakeSeveralIC = kFALSE ;
+  fMakeInvMass = kFALSE ;
 
-  fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
-  fAnaMC = kFALSE ;
-  fIsolatePi0 = kFALSE ;
- //-----------kSeveralIC-----------------
+ //----------- Several IC-----------------
   fNCones           = 5 ; 
-  fNPtThres         = 6 ; 
+  fNPtThresFrac         = 6 ; 
   fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
   fPtThresholds[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
+  fPtFractions[0]=0.05; fPtFractions[1]=0.075; fPtFractions[2]=0.1; fPtFractions[3]=1.25; fPtFractions[4]=1.5;fPtFractions[5]=2.;
 
 }
 
 //__________________________________________________________________
-void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS, 
-                                         TClonesArray * plNe, 
-                                         TParticle * pCandidate, 
-                                         Int_t index, Int_t & n,
-                                         Bool_t  &icmpt,  Bool_t  &icms, 
-                                         Float_t &coneptsum) const 
-{  
-  //Search in cone around a candidate particle if it is isolated 
-  Float_t phiC  = pCandidate->Phi() ;
-  Float_t etaC = pCandidate->Eta() ;
-  Float_t pt     = -100. ;
-  Float_t eta   = -100.  ;
-  Float_t phi    = -100.  ;
-  Float_t rad   = -100 ;
-  TParticle * particle  = new TParticle;
-
-  n = 0 ;
-  coneptsum = 0.; 
-  icmpt = kFALSE;
-  icms  = kFALSE;
-
-  //Check charged particles in cone.
-  for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
-    particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
-    pt    = particle->Pt();
-    eta  = particle->Eta();
-    phi  = particle->Phi() ;
-    
-    //Check if there is any particle inside cone with pt larger than  fPtThreshold
-    rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
-    if(rad < fConeSize){
-      AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
-      coneptsum+=pt;
-      if(pt > fPtThreshold ) n++;
+void  AliAnaGammaDirect::MakeSeveralICAnalysis(AliAODParticleCorrelation* ph, Double_t v[3]) 
+{
+  //Isolation Cut Analysis for both methods and different pt cuts and cones
+  Float_t ptC = ph->Pt();
+  Float_t phiC = ph->Phi();
+  Float_t etaC = ph->Eta();
+
+  fhPtGamma->Fill(ptC);
+  fhPhiGamma->Fill(ptC,ph->Phi());
+  fhEtaGamma->Fill(ptC,ph->Eta());
+  Int_t tag =ph->GetTag();
+
+  if(IsDataMC()){
+    if(tag == kPrompt){
+      fhPtPrompt   ->Fill(ptC);
+      fhPhiPrompt ->Fill(ptC,phiC);
+      fhEtaPrompt ->Fill(ptC,etaC);
     }
-  }// charged particle loop
-  
-  //Check neutral particles in cone.
-  for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
-    if(ipr != index){//Do not count the candidate
-      particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
-      pt    = particle->Pt();
-      eta  = particle->Eta();
-      phi  = particle->Phi() ;
+    else if(tag==kFragmentation)
+      {
+       fhPtFragmentation   ->Fill(ptC);
+       fhPhiFragmentation ->Fill(ptC,phiC);
+       fhEtaFragmentation ->Fill(ptC,etaC);
+      }
+    else if(tag==kPi0Decay)
+      {
+       fhPtPi0Decay   ->Fill(ptC);
+       fhPhiPi0Decay ->Fill(ptC,phiC);
+       fhEtaPi0Decay ->Fill(ptC,etaC);
+      }
+    else if(tag==kEtaDecay || tag==kOtherDecay)
+      {
+       fhPtOtherDecay   ->Fill(ptC);
+       fhPhiOtherDecay ->Fill(ptC,phiC);
+       fhEtaOtherDecay ->Fill(ptC,etaC);
+      }
+    else if(tag==kConversion)
+      {
+       fhPtConversion   ->Fill(ptC);
+       fhPhiConversion ->Fill(ptC,phiC);
+       fhEtaConversion ->Fill(ptC,etaC);
+      }
+    else{
       
-      //Check if there is any particle inside cone with pt larger than  fPtThreshold
-      rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
-      if(rad < fConeSize){
-       AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
-       coneptsum+=pt;
-       if(pt > fPtThreshold ) n++;
+      fhPtUnknown   ->Fill(ptC);
+      fhPhiUnknown ->Fill(ptC,phiC);
+      fhEtaUnknown ->Fill(ptC,etaC);
+    }
+  }//Histograms with MC
+  //Keep original setting used when filling AODs, reset at end of analysis  
+  Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
+  Float_t ptfracorg = GetIsolationCut()->GetPtFraction();
+  Float_t rorg = GetIsolationCut()->GetConeSize();
+
+  Float_t coneptsum = 0 ;  
+  Int_t n[10][10];//[fNCones][fNPtThresFrac];
+  Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];
+  Bool_t  isolated   = kFALSE;
+
+  for(Int_t icone = 0; icone<fNCones; icone++){
+    GetIsolationCut()->SetConeSize(fConeSizes[icone]);
+    coneptsum = 0 ;
+    for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){
+      n[icone][ipt]=0;
+      nfrac[icone][ipt]=0;
+      GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
+      GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(), 
+                                         (TSeqCollection*)ph->GetRefIsolationConeClusters(), 
+                                         v, kFALSE, ph,-1,-1,
+                                         n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
+      if(n[icone][ipt] == 0) {
+       fhPtThresIsolated[icone][ipt]->Fill(ptC);
+       if(IsDataMC()){
+         if(tag==kPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
+         else if(tag==kConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
+         else if(tag==kFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
+         else if(tag==kPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
+         else if(tag==kOtherDecay || tag==kEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
+         else  fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
+       }
+      }
+      if(nfrac[icone][ipt] == 0) {
+       fhPtFracIsolated[icone][ipt]->Fill(ptC);
+       if(IsDataMC()){
+         if(tag==kPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
+         else if(tag==kConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
+         else if(tag==kFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
+         else if(tag==kPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
+         else if(tag==kOtherDecay || tag==kEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
+         else  fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
+       }
       }
+    }//pt thresh loop
+    fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
+    if(IsDataMC()){
+      if(tag==kPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;
+      else if(tag==kConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;
+      else if(tag==kFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
+      else if(tag==kPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;
+      else if(tag==kOtherDecay || tag==kEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;
+      else  fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;
     }
-  }// neutral particle loop
-  
-  if(n == 0) 
-    icmpt =  kTRUE ;
-  if(coneptsum <= fPtSumThreshold)
-    icms  =  kTRUE ;
 
-}
+  }//cone size loop
+
+  //Reset original parameters for AOD analysis
+  GetIsolationCut()->SetPtThreshold(ptthresorg);
+  GetIsolationCut()->SetPtFraction(ptfracorg);
+  GetIsolationCut()->SetConeSize(rorg);
 
-//__________________________________________________________________
-void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS) 
-{
-  //Isolation Cut Analysis for both methods and different pt cuts and cones
-  if (fICMethod != kSeveralIC)
-    AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
-
-  //Search maximum energy photon in the event
-  Double_t ptC = 0;
-  Int_t index = -1; 
-  TParticle * pCandidate = new TParticle();
-  Bool_t found = kFALSE;
-  
-  for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
-    TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
-    
-    if((particle->Pt() > fMinGammaPt) && (particle->Pt() > ptC) && 
-       (particle->GetPdgCode() == 22 ||  (fIsolatePi0 && particle->GetPdgCode() == 111))){
-      index = ipr ;
-      ptC = particle->Pt();
-      pCandidate = particle ;
-      found  = kTRUE;
-    }
-  }
-  
-  //If there is a large cluster, larger than threshold, study isolation cut
-  if(found){
-    
-    fhNGamma->Fill(ptC);
-    fhPhiGamma->Fill(ptC,pCandidate->Phi());
-    fhEtaGamma->Fill(ptC,pCandidate->Eta());
-    
-    Int_t ncone[10][10];//[fNCones][fNPtThres];
-    Bool_t  icPtThres   = kFALSE;
-    Bool_t  icPtSum     = kFALSE;
-    
-    for(Int_t icone = 0; icone<fNCones; icone++){
-      fConeSize=fConeSizes[icone] ;
-      Float_t coneptsum = 0 ;
-      for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
-       ncone[icone][ipt]=0;
-       fPtThreshold=fPtThresholds[ipt] ;
-       MakeIsolationCut(plCTS,plCalo, pCandidate, index,  
-                        ncone[icone][ipt], icPtThres, icPtSum,coneptsum);
-       AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
-                       pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
-//     if(ptC >15 && ptC < 25 && (icPtThres || icPtSum) && ipt ==0){
-//       printf("R %0.1f, ptthres %1.1f, ptsum %1.1f, Candidate pt %2.2f,  eta %2.2f, phi %2.2f, pt in cone %2.2f, Isolated? ICPt %d, ICSum %d\n",
-//              fConeSize,  fPtThreshold, fPtSumThreshold, ptC, pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), coneptsum, icPtThres, icPtSum);
-//       //cout<<"mother label "<<pCandidate->GetMother(0)<<endl;
-//     }
-       
-       fhPtThresIsolated[icone][ipt]->Fill(ptC); 
-      }//pt thresh loop
-      fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
-//      gROOT->cd();
-      fntSeveralIC[icone]->Fill(ptC,pCandidate->Phi(),pCandidate->Eta(),  pCandidate->GetPdgCode(), pCandidate->GetStatusCode(), coneptsum,ncone[icone][0],ncone[icone][1],ncone[icone][2],ncone[icone][3],ncone[icone][4],ncone[icone][5]);
-    }//cone size loop
-  }//found high energy gamma in the event
 }
 
 //__________________________________________________________________
@@ -458,25 +946,74 @@ void AliAnaGammaDirect::Print(const Option_t * opt) const
   if(! opt)
     return;
   
-  Info("Print", "%s %s", GetName(), GetTitle() ) ;
+  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
   
-  printf("Min Gamma pT      =     %f\n",  fMinGammaPt) ;
-  printf("IC method               =     %d\n", fICMethod) ; 
-  printf("Cone Size               =     %f\n", fConeSize) ; 
-   if(fICMethod == kPtIC) printf("pT threshold           =     %f\n", fPtThreshold) ;
-   if(fICMethod == kSumPtIC) printf("pT sum threshold   =     %f\n", fPtSumThreshold) ;
-   
-  if(fICMethod == kSeveralIC){
-    printf("N Cone Sizes               =     %d\n", fNCones) ; 
-    printf("N pT thresholds           =     %d\n", fNPtThres) ;
-    printf("Cone Sizes                  =    \n") ;
+  printf("Min Gamma pT       =     %2.2f\n",  GetMinPt()) ;
+  printf("Max Gamma pT       =     %3.2f\n",  GetMaxPt()) ;
+
+//   if(IsCaloPIDOn())printf("Check PID \n") ;
+//   if(IsCaloPIDRecalculationOn())printf("Recalculate PID \n") ;
+//   if(IsFidutialCutOn())printf("Check Fidutial cut \n") ;
+  printf("Make Isolation     =     %d \n",  fMakeIC) ;
+  printf("ReMake Isolation   =     %d \n",  fReMakeIC) ;
+  printf("Make Several Isolation = %d \n",  fMakeSeveralIC) ;
+  
+  if(fMakeSeveralIC){
+    printf("N Cone Sizes       =     %d\n", fNCones) ; 
+    printf("Cone Sizes          =    \n") ;
     for(Int_t i = 0; i < fNCones; i++)
-      printf("   %f;",  fConeSizes[i]) ;
+      printf("  %1.2f;",  fConeSizes[i]) ;
     printf("    \n") ;
-    for(Int_t i = 0; i < fNPtThres; i++)
-      printf("   %f;",  fPtThresholds[i]) ;
-  }
-
+    
+    printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
+    printf(" pT thresholds         =    \n") ;
+    for(Int_t i = 0; i < fNPtThresFrac; i++)
+      printf("   %2.2f;",  fPtThresholds[i]) ;
+    
+    printf("    \n") ;
+    
+    printf(" pT fractions         =    \n") ;
+    for(Int_t i = 0; i < fNPtThresFrac; i++)
+      printf("   %2.2f;",  fPtFractions[i]) ;
+    
+  }  
+  
   printf("    \n") ;
   
 } 
+
+//____________________________________________________________________________
+Bool_t  AliAnaGammaDirect::SelectCluster(AliAODCaloCluster * calo, Double_t vertex[3], TLorentzVector & mom){
+   //Select cluster depending on its pid and acceptance selections
+   
+   //Skip matched clusters with tracks
+   if(calo->GetNTracksMatched() > 0) return kFALSE ;
+   
+   //Check PID
+   calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
+   Int_t pdg = AliCaloPID::kPhoton;   
+   if(IsCaloPIDOn()){
+     //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
+     //or redo PID, recommended option for EMCal.
+     if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
+       pdg = GetCaloPID()->GetPdg(fDetector,calo->PID(),mom.E());//PID with weights
+     else
+       pdg = GetCaloPID()->GetPdg(fDetector,mom,calo->GetM02(),0,0,0,0);//PID recalculated
+     
+     if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
+     
+     //If it does not pass pid, skip
+     if(pdg != AliCaloPID::kPhoton) return kFALSE ;
+   }
+   
+   //Check acceptance selection
+   if(IsFidutialCutOn()){
+     Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
+     if(! in ) return kFALSE ;
+   }
+   
+   if(GetDebug() > 1) printf("Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
+   
+   return kTRUE;
+   
+ }
index 0515353..b0d03dc 100644 (file)
@@ -1,5 +1,5 @@
-#ifndef ALIANAGAMMADIRECT_H
-#define ALIANAGAMMADIRECT_H
+#ifndef AliAnaGammaDirect_H
+#define AliAnaGammaDirect_H
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice     */
 /* $Id$ */
 //  Class created from old AliPHOSGammaJet
 //  (see AliRoot versions previous Release 4-09)
 
-//*-- Author: Gustavo Conesa (INFN-LNF)
+//-- Author: Gustavo Conesa (INFN-LNF)
 
 // --- ROOT system ---
 #include <TParticle.h> 
-#include <TClonesArray.h> 
-#include "TObject.h" 
+#include <TClonesArray.h>  
 #include <TH2F.h>
-#include <TNtuple.h>
+#include <TString.h>
+
+// --- ANALYSIS system ---
+#include "AliAnaBaseClass.h"
+class AliAODParticleCorrelations ;
 
 class TList ;
 
-class AliAnaGammaDirect : public TObject {
+class AliAnaGammaDirect : public AliAnaBaseClass {
 
 public: 
 
@@ -42,83 +45,125 @@ public:
   AliAnaGammaDirect(const AliAnaGammaDirect & g) ; // cpy ctor
   AliAnaGammaDirect & operator = (const AliAnaGammaDirect & g) ;//cpy assignment
   virtual ~AliAnaGammaDirect() ; //virtual dtor
-  
-  enum Anatype {kNoIC, kPtIC, kSumPtIC, kSeveralIC};
-  
-  Double_t  GetMinGammaPt()    const {return fMinGammaPt ; }
-  Float_t     GetConeSize()          const {return fConeSize ; }
-  Float_t     GetPtThreshold()      const {return fPtThreshold ; }
-  Float_t     GetPtSumThres()     const {return fPtSumThreshold ; }
-  Int_t        GetICMethod()          const {return fICMethod ; }
-  Bool_t     IsMC() const {return fAnaMC ; };
 
+  enum mcTypes {kPrompt, kFragmentation, kPi0Decay, kEtaDecay, kOtherDecay, kPi0, kEta, kElectron, kConversion, kUnknown};
+
+  Bool_t CheckInvMass(const Int_t icalo,const TLorentzVector mom, Double_t *v, TClonesArray * pl);
+  Int_t CheckOrigin(const Int_t label);
+  
   TList *  GetCreateOutputObjects();
-  void GetPromptGamma(TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plPrimNe,  TParticle * pGamma, Bool_t &Is)  const;
+
+  void MakeAnalysisFillAOD()  ;
   
-  void MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS); 
-  void MakeIsolationCut(TClonesArray * plCTS, TClonesArray * plNe, 
-                       TParticle *pCandidate, Int_t index, Int_t &n,
-                       Bool_t &imcpt, Bool_t &icms, Float_t &ptsum) const ;  
+  void MakeAnalysisFillHistograms() ; 
+
+  void MakeSeveralICAnalysis(AliAODParticleCorrelation * ph, Double_t v[3]); 
   
   void Print(const Option_t * opt)const;
   
-  void SetMinGammaPt(Double_t ptcut){fMinGammaPt =ptcut;}
-  void SetConeSize(Float_t r)              {fConeSize = r ; }
-  void SetPtThreshold(Float_t pt)        {fPtThreshold = pt; };
-  void SetPtSumThreshold(Float_t pt) {fPtSumThreshold = pt; };
-  void SetICMethod(Int_t i )          {fICMethod = i ; }
-  void SetMC()    {fAnaMC = kTRUE ; }
-  void SetIsolatePi0(Bool_t iso)    {fIsolatePi0 = iso ; }
+  TString GetDetector()   const {return fDetector ; }
+  void SetDetector(TString det)    {fDetector = det ; }
 
   Int_t    GetNCones()                  const {return fNCones ; }
-  Int_t    GetNPtThresholds()                const {return fNPtThres ; }
+  Int_t    GetNPtThresFrac()                const {return fNPtThresFrac ; }
   Float_t GetConeSizes(Int_t i)      const {return fConeSizes[i] ; }
   Float_t GetPtThresholds(Int_t i)  const {return fPtThresholds[i] ; }
-  
+  Float_t GetPtFractions(Int_t i)  const {return fPtFractions[i] ; }
+
   void InitParameters();
  
   void SetNCones(Int_t ncs)              {fNCones = ncs ; }
-  void SetNPtThresholds(Int_t npt)        {fNPtThres = npt; }
+  void SetNPtThresFrac(Int_t npt)        {fNPtThresFrac = npt; }
   void SetConeSizes(Int_t i, Float_t r)         {fConeSizes[i] = r ; }
   void SetPtThresholds(Int_t i, Float_t pt)   {fPtThresholds[i] = pt; }
+  void SetPtFractions(Int_t i, Float_t pt)   {fPtFractions[i] = pt; } 
+
+  Bool_t IsIsolationOn() {return fMakeIC ; }
+  void SwitchOnIsolation() { fMakeIC = kTRUE;}
+  void SwitchOffIsolation() { fMakeIC = kFALSE;}
+
+  Bool_t IsReIsolationOn() {return fReMakeIC ; }
+  void SwitchOnReIsolation() { fReMakeIC = kTRUE;}
+  void SwitchOffReIsolation() { fReMakeIC = kFALSE;}
+
+  Bool_t IsSeveralIsolationOn() {return fMakeSeveralIC ; }
+  void SwitchOnSeveralIsolation() { fMakeSeveralIC = kTRUE;}
+  void SwitchOffSeveralIsolation() { fMakeSeveralIC = kFALSE;}
+
+  Bool_t IsInvariantMassOn() {return fMakeInvMass ; }
+  void SwitchOnInvariantMass() { fMakeInvMass = kTRUE;}
+  void SwitchOffInvariantMass() { fMakeInvMass = kFALSE;}
+
+  Bool_t SelectCluster(AliAODCaloCluster * calo, Double_t vertex[3], TLorentzVector & mom);
 
-  
   private:
-     
-  Double_t    fMinGammaPt ;  // Min pt in Calorimeter
-  Float_t      fConeSize ; //Size of the isolation cone 
-  Float_t      fPtThreshold ; //Mimium pt of the particles in the cone to set isolation
-  Float_t      fPtSumThreshold ; //Mimium pt sum of the particles in the cone to set isolation  
-  Int_t        fICMethod ; //Isolation cut method to be used
-                                           // kNoIC: No isolation
-                                           // kPtIC: Pt threshold method
-                                           // kSumPtIC: Cone pt sum method
-                                           // kSeveralIC: Analysis for several cuts
-  Bool_t fAnaMC ; //Set in case of using MCData reader 
-  Bool_t fIsolatePi0 ; //Consider identified pi0 in the isolation study.
+  TString fDetector ; // Detector where the gamma is searched;
+  Bool_t fMakeIC ; //Do isolation analysis
+  Bool_t fReMakeIC ; //Do isolation analysis
+  Bool_t fMakeSeveralIC ; //Do analysis for different IC
+  Bool_t fMakeInvMass; //Select candidate if no pair from decay
+
   //Histograms  
-  TH1F * fhNGamma    ;  //Number of (isolated) gamma identified
-  TH2F * fhPhiGamma  ; // Phi of identified gamma
-  TH2F * fhEtaGamma  ; // eta of identified gamma
+  TH1F * fhPtGamma    ;  //Number of identified (isolated) gamma 
+  TH2F * fhPhiGamma  ; // Phi of identified  (isolated) gamma
+  TH2F * fhEtaGamma  ; // eta of identified  (isolated) gamma
   TH2F * fhConeSumPt ; // Sum Pt in the cone
 
-  TNtuple *    fntuplePrompt ; //List of found prompt photons, pt, eta and phi. Also primary information.
-
-  //Prompt photon analysis data members for multiple cones and pt thresholds kIsolationCut
+  //Prompt photon analysis data members for multiple cones and pt thresholds 
   Int_t         fNCones   ; //Number of cone sizes to test
-  Int_t         fNPtThres ; //Number of ptThres to test
-  Float_t     fConeSizes[10] ; // Array with cones to test
-  Float_t     fPtThresholds[10] ; // Array with pt thresholds to test
-  
-  TH1F* fhPtThresIsolated[20][20]; // Isolated gamma with pt threshold 
-  TH2F* fhPtSumIsolated[20] ;  //  Isolated gamma with threshold on cone pt sume
-  TNtuple *    fntSeveralIC[20] ; //ntuple 
-
+  Int_t         fNPtThresFrac ; //Number of ptThres and ptFrac to test
+  Float_t     fConeSizes[5] ; // Array with cones to test
+  Float_t     fPtThresholds[5] ; // Array with pt thresholds to test
+  Float_t     fPtFractions[5] ; // Array with pt thresholds to test
+
+  TH1F* fhPtThresIsolated[5][5]; // Isolated gamma with pt threshold 
+  TH1F* fhPtFracIsolated[5][5]; // Isolated gamma with pt threshold 
+  TH2F* fhPtSumIsolated[5] ;  //  Isolated gamma with threshold on cone pt sume
+
+  //MC
+  TH1F * fhPtPrompt; //Number of identified (isolated) prompt gamma 
+  TH2F * fhPhiPrompt;  // Phi of identified  (isolated) prompt gamma
+  TH2F * fhEtaPrompt;  // eta of identified  (isolated) prompt gamma
+  TH1F * fhPtThresIsolatedPrompt[5][5];  // Isolated prompt gamma with pt threshold 
+  TH1F * fhPtFracIsolatedPrompt[5][5];    // Isolated prompt gamma with pt frac
+  TH2F * fhPtSumIsolatedPrompt[5]; //  Isolated prompt gamma with threshold on cone pt sume
+  TH1F * fhPtFragmentation; //Number of identified (isolated) fragmentation gamma 
+  TH2F * fhPhiFragmentation;  // Phi of identified  (isolated) fragmentation gamma
+  TH2F * fhEtaFragmentation;  // eta of identified  (isolated) fragmentation gamma
+  TH1F * fhPtThresIsolatedFragmentation[5][5];  // Isolated fragmentation gamma with pt threshold 
+  TH1F * fhPtFracIsolatedFragmentation[5][5];    // Isolated fragmentation gamma with pt frac
+  TH2F * fhPtSumIsolatedFragmentation[5]; //  Isolated fragmentation gamma with threshold on cone pt sume
+  TH1F * fhPtPi0Decay; //Number of identified (isolated) Pi0Decay gamma 
+  TH2F * fhPhiPi0Decay;  // Phi of identified  (isolated) Pi0Decay gamma
+  TH2F * fhEtaPi0Decay;  // eta of identified  (isolated) Pi0Decay gamma
+  TH1F * fhPtThresIsolatedPi0Decay[5][5];  // Isolated Pi0Decay gamma with pt threshold 
+  TH1F * fhPtFracIsolatedPi0Decay[5][5];    // Isolated Pi0Decay gamma with pt frac
+  TH2F * fhPtSumIsolatedPi0Decay[5]; //  Isolated Pi0Decay gamma with threshold on cone pt sume
+  TH1F * fhPtOtherDecay; //Number of identified (isolated) OtherDecay gamma 
+  TH2F * fhPhiOtherDecay;  // Phi of identified  (isolated) OtherDecay gamma
+  TH2F * fhEtaOtherDecay;  // eta of identified  (isolated) OtherDecay gamma
+  TH1F * fhPtThresIsolatedOtherDecay[5][5];  // Isolated OtherDecay gamma with pt threshold 
+  TH1F * fhPtFracIsolatedOtherDecay[5][5];    // Isolated OtherDecay gamma with pt frac
+  TH2F * fhPtSumIsolatedOtherDecay[5]; //  Isolated OtherDecay gamma with threshold on cone pt sume    
+  TH1F * fhPtConversion; //Number of identified (isolated) Conversion gamma 
+  TH2F * fhPhiConversion;  // Phi of identified  (isolated) Conversion gamma
+  TH2F * fhEtaConversion;  // eta of identified  (isolated) Conversion gamma
+  TH1F * fhPtThresIsolatedConversion[5][5];  // Isolated Conversion gamma with pt threshold 
+  TH1F * fhPtFracIsolatedConversion[5][5];    // Isolated Conversion gamma with pt frac
+  TH2F * fhPtSumIsolatedConversion[5]; //  Isolated Conversion gamma with threshold on cone pt sume
+  TH1F * fhPtUnknown; //Number of identified (isolated) Unknown gamma 
+  TH2F * fhPhiUnknown;  // Phi of identified  (isolated) Unknown gamma
+  TH2F * fhEtaUnknown;  // eta of identified  (isolated) Unknown gamma
+  TH1F * fhPtThresIsolatedUnknown[5][5];  // Isolated Unknown gamma with pt threshold 
+  TH1F * fhPtFracIsolatedUnknown[5][5];    // Isolated Unknown gamma with pt frac
+  TH2F * fhPtSumIsolatedUnknown[5]; //  Isolated Unknown gamma with threshold on cone pt sume
+                                               
   ClassDef(AliAnaGammaDirect,1)
 } ;
  
 
-#endif //ALIANAGAMMADIRECT_H
+#endif //AliAnaGammaDirect_H
 
 
 
index f4a2f60..d01f687 100644 (file)
@@ -94,7 +94,7 @@ void AliAnaScale::Exec(Option_t *)
   while ( (h = next()) ) { 
     if(h){
       if ( strcmp(h->ClassName(),"TNtuple") ) {
-      char name[20] ; 
+      char name[128] ; 
       sprintf(name, "%sScaled", h->GetName()) ; 
       TH1 * hout = dynamic_cast<TH1*> (h->Clone(name)) ; 
       hout->Scale(fScale) ;  
@@ -103,7 +103,7 @@ void AliAnaScale::Exec(Option_t *)
       else  fOutputList->Add(h) ; 
     }
   }
-  cout<<"end"<<endl;
+  
   PostData(0, fOutputList);
 }
 
index f0cc48a..20065cb 100644 (file)
 // Class that contains methods to select candidate pairs to neutral meson 
 // 2 main selections, invariant mass around pi0 (also any other mass),
 // apperture angle to distinguish from combinatorial.
-// There is a 3rd cut based on the gamma correlation on phi or pt.
 //-- Author: Gustavo Conesa (INFN-LNF)
 
 // --- ROOT system ---
-#include <TParticle.h>
 #include <TLorentzVector.h>
 #include <TH2.h>
 #include <TList.h>
@@ -52,15 +50,12 @@ ClassImp(AliNeutralMesonSelection)
   
 //____________________________________________________________________________
   AliNeutralMesonSelection::AliNeutralMesonSelection() : 
-    TObject(), fSelect(0), fM(0),
+    TObject(), fM(0),
     fInvMassMaxCut(0.), fInvMassMinCut(0.),
-    fAngleMaxParam(),  fMinPt(0),
-    fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
-    fRatioMaxCut(0), fRatioMinCut(0),  fKeepNeutralMesonHistos(0),
-    fhAnglePairNoCut(0),  fhAnglePairCorrelationCut(0),
+    fAngleMaxParam(), fhAnglePairNoCut(0),
     fhAnglePairOpeningAngleCut(0), fhAnglePairAllCut(0), 
-    fhInvMassPairNoCut(0),   fhInvMassPairCorrelationCut(0), 
-    fhInvMassPairOpeningAngleCut(0), fhInvMassPairAllCut(0) 
+    fhInvMassPairNoCut(0), fhInvMassPairOpeningAngleCut(0), 
+    fhInvMassPairAllCut(0) 
 {
   //Default Ctor
   
@@ -76,19 +71,14 @@ ClassImp(AliNeutralMesonSelection)
 
 //____________________________________________________________________________
 AliNeutralMesonSelection::AliNeutralMesonSelection(const AliNeutralMesonSelection & g) :   
-  TObject(),  
-  fSelect(g.fSelect), fM(g.fM),
+  TObject(), fM(g.fM),
   fInvMassMaxCut(g.fInvMassMaxCut), fInvMassMinCut(g.fInvMassMinCut),
-  fAngleMaxParam(g.fAngleMaxParam), fMinPt(g.fMinPt),
-  fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut), 
-  fRatioMaxCut(g.fRatioMaxCut), fRatioMinCut(g.fRatioMinCut), 
+  fAngleMaxParam(g.fAngleMaxParam),
   fKeepNeutralMesonHistos(g.fKeepNeutralMesonHistos),
   fhAnglePairNoCut(g. fhAnglePairNoCut), 
-  fhAnglePairCorrelationCut(g. fhAnglePairCorrelationCut), 
   fhAnglePairOpeningAngleCut(g. fhAnglePairOpeningAngleCut), 
   fhAnglePairAllCut(g. fhAnglePairAllCut), 
   fhInvMassPairNoCut(g.fhInvMassPairNoCut),  
-  fhInvMassPairCorrelationCut(g.fhInvMassPairCorrelationCut), 
   fhInvMassPairOpeningAngleCut(g.fhInvMassPairOpeningAngleCut), 
   fhInvMassPairAllCut(g.fhInvMassPairAllCut)
 {
@@ -103,24 +93,16 @@ AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutra
   if(this == &source)return *this;
   ((TObject *)this)->operator=(source);
 
-  fSelect = source.fSelect ;
   fM = source.fM ;
   fInvMassMaxCut = source.fInvMassMaxCut ; 
   fInvMassMinCut = source.fInvMassMinCut ;
   fAngleMaxParam = source.fAngleMaxParam ;
-  fMinPt = source.fMinPt ;
-  fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ; 
-  fDeltaPhiMinCut = source.fDeltaPhiMinCut ; 
-  fRatioMaxCut = source.fRatioMaxCut ; 
-  fRatioMinCut = source.fRatioMinCut ;
   fKeepNeutralMesonHistos = source.fKeepNeutralMesonHistos;
  
   fhAnglePairNoCut = source. fhAnglePairNoCut ; 
-  fhAnglePairCorrelationCut = source. fhAnglePairCorrelationCut ; 
   fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ; 
   fhAnglePairAllCut = source. fhAnglePairAllCut ; 
   fhInvMassPairNoCut = source.fhInvMassPairNoCut ; 
-  fhInvMassPairCorrelationCut = source.fhInvMassPairCorrelationCut ; 
   fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ; 
   fhInvMassPairAllCut = source.fhInvMassPairAllCut ; 
   
@@ -131,18 +113,29 @@ AliNeutralMesonSelection & AliNeutralMesonSelection::operator = (const AliNeutra
 //____________________________________________________________________________
 AliNeutralMesonSelection::~AliNeutralMesonSelection() 
 {
- // Remove all pointers except analysis output pointers.
-
-}
-
+  //dtor
 
+  if(!fKeepNeutralMesonHistos){
+    //Histograms initialized and filled but not passed to output container
+    //delete here, I am not sure this is correct
+    
+    if(fhAnglePairNoCut) delete fhAnglePairNoCut;
+    if(fhAnglePairOpeningAngleCut) delete fhAnglePairOpeningAngleCut; 
+    if(fhAnglePairAllCut) delete fhAnglePairAllCut;
+    if(fhInvMassPairNoCut) delete fhInvMassPairNoCut;
+    if(fhInvMassPairOpeningAngleCut) delete fhInvMassPairOpeningAngleCut;
+    if(fhInvMassPairAllCut) delete fhInvMassPairAllCut; 
 
+  }
+  
+}
 //________________________________________________________________________
 TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
 {  
 
   // Create histograms to be saved in output file and 
-  // store them in outputContainer
+  // store them in outputContainer of the analysis class that calls this class.
+
   TList * outputContainer = new TList() ; 
   outputContainer->SetName("MesonDecayHistos") ; 
   
@@ -187,18 +180,6 @@ TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
   fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
   fhInvMassPairAllCut->SetXTitle("E_{#pi^{0}}(GeV)");
 
-  fhAnglePairCorrelationCut  = new TH2F
-    ("AnglePairCorrelationCut",
-     "Angle between correlated #gamma pair vs E_{#pi^{0}}",200,0,50,200,0,0.2); 
-  fhAnglePairCorrelationCut->SetYTitle("Angle (rad)");
-  fhAnglePairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");
-
-  fhInvMassPairCorrelationCut  = new TH2F
-    ("InvMassPairCorrelationCut","Invariant Mass of correlated #gamma pair vs E_{#pi^{0}}",
-     120,0,120,360,0,0.5); 
-  fhInvMassPairCorrelationCut->SetYTitle("Invariant Mass (GeV/c^{2})");
-  fhInvMassPairCorrelationCut->SetXTitle("E_{ #pi^{0}} (GeV)");
-  
   outputContainer->Add(fhAnglePairNoCut) ; 
   outputContainer->Add(fhAnglePairOpeningAngleCut) ;
   outputContainer->Add(fhAnglePairAllCut) ; 
@@ -206,9 +187,6 @@ TList *  AliNeutralMesonSelection::GetCreateOutputObjects()
   outputContainer->Add(fhInvMassPairNoCut) ; 
   outputContainer->Add(fhInvMassPairOpeningAngleCut) ; 
   outputContainer->Add(fhInvMassPairAllCut) ; 
-
-  outputContainer->Add(fhAnglePairCorrelationCut) ; 
-  outputContainer->Add(fhInvMassPairCorrelationCut) ; 
   
   return outputContainer;
 }
@@ -218,7 +196,7 @@ void AliNeutralMesonSelection::InitParameters()
 {
  
   //Initialize the parameters of the analysis.
-  fKeepNeutralMesonHistos = kTRUE ;
+  fKeepNeutralMesonHistos = kFALSE ;
 
   //-------------kHadron, kJetLeadCone-----------------
   fAngleMaxParam.Set(4) ;
@@ -231,12 +209,6 @@ void AliNeutralMesonSelection::InitParameters()
   fInvMassMinCut  = 0.11 ;
 
   fM = 0.1349766;//neutralMeson mass
-
-  fMinPt = 0.   ;
-  fDeltaPhiMaxCut      = 4.5;
-  fDeltaPhiMinCut      = 1.5 ;
-  fRatioMaxCut    = 1.0 ;
-  fRatioMinCut    = 0.1 ; 
 }
 
 //__________________________________________________________________________-
@@ -259,31 +231,10 @@ Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float
 }
 
 //____________________________________________________________________________
-Bool_t  AliNeutralMesonSelection::CutPtPhi(Double_t ptg, Double_t phig, Double_t pt, Double_t phi)  const
-{ 
-  //Select pair if delta
-  Bool_t cut = kFALSE ;
-  if(fSelect == kNoSelectPhiPt) cut = kTRUE ;
-  else if((phig-phi) > fDeltaPhiMinCut && ((phig-phi) < fDeltaPhiMaxCut)){
-    //Cut on pt
-    if((fSelect == kSelectPhiPtRatio && ptg > 0. && pt/ptg  > fRatioMinCut &&  pt/ptg  < fRatioMaxCut) ||
-       (fSelect == kSelectPhiMinPt && pt > fMinPt)  )  cut = kTRUE ;
-  }
-  else cut = kFALSE ;
-  
-  return cut ;
-  
-}
-
-//____________________________________________________________________________
-Bool_t  AliNeutralMesonSelection::SelectPair(TParticle * pGamma, TLorentzVector gammai, TLorentzVector gammaj)  
+Bool_t  AliNeutralMesonSelection::SelectPair(TLorentzVector gammai, TLorentzVector gammaj)  
 {  
   
   //Search for the neutral pion within selection cuts
-  
-  Double_t ptg = pGamma->Pt();
-  Double_t phig = pGamma->Phi() ;
   Bool_t goodpair = kFALSE ;
   
   Double_t pt  = (gammai+gammaj).Pt();
@@ -297,29 +248,21 @@ Bool_t  AliNeutralMesonSelection::SelectPair(TParticle * pGamma, TLorentzVector
   //Fill histograms with no cuts applied.
   fhAnglePairNoCut->Fill(e,angle);
   fhInvMassPairNoCut->Fill(e,invmass);
-  
-  //Cut on phig-phi meson
-  if(CutPtPhi(ptg, phig, pt, phi)){
     
-    fhAnglePairCorrelationCut     ->Fill(e,angle);
-    fhInvMassPairCorrelationCut->Fill(e,invmass);
+  //Cut on the aperture of the pair
+  if(IsAngleInWindow(angle,e)){
+    fhAnglePairOpeningAngleCut     ->Fill(e,angle);
+    fhInvMassPairOpeningAngleCut->Fill(e,invmass);
+    AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
     
-    //Cut on the aperture of the pair
-    if(IsAngleInWindow(angle,e)){
-      fhAnglePairOpeningAngleCut     ->Fill(e,angle);
-      fhInvMassPairOpeningAngleCut->Fill(e,invmass);
-      AliDebug(2,Form("Angle cut: pt %f, phi %f",pt,phi));
-      
-      //Cut on the invariant mass of the pair
-      if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ 
-       fhInvMassPairAllCut  ->Fill(e,invmass);
-       fhAnglePairAllCut       ->Fill(e,angle);
-       goodpair = kTRUE;
-       AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
-      }//(invmass>0.125) && (invmass<0.145)
-    }//Opening angle cut
-  } // cut on pt and phi
-  
+    //Cut on the invariant mass of the pair
+    if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ 
+      fhInvMassPairAllCut  ->Fill(e,invmass);
+      fhAnglePairAllCut       ->Fill(e,angle);
+      goodpair = kTRUE;
+      AliDebug(2,Form("IM cut: pt %f, phi %f",pt,phi));
+    }//(invmass>0.125) && (invmass<0.145)
+  }//Opening angle cut
   
   return goodpair; 
   
@@ -343,11 +286,6 @@ void AliNeutralMesonSelection::Print(const Option_t * opt) const
   printf("p2 :     %f", fAngleMaxParam.At(2));
   printf("p3 :     %f", fAngleMaxParam.At(3));
 
-  printf("pT meson       >    %f\n", fMinPt) ; 
-  printf("Phi gamma-meson      <     %f\n", fDeltaPhiMaxCut) ; 
-  printf("Phi gamma-meson      >     %f\n", fDeltaPhiMinCut) ;
-  printf("pT meson / pT Gamma             <     %f\n", fRatioMaxCut) ; 
-  printf("pT meson / pT Gamma             >     %f\n", fRatioMinCut) ;
   printf("Keep Neutral Meson Histos = %d\n",fKeepNeutralMesonHistos);
 
 } 
index 4767e8a..dde3462 100644 (file)
 #include<TArrayD.h>
 
 class TLorentzVector ;
-class TParticle ;
 class TList ;
 class TH2F ;
 class Riostream ;
 
-//--- AliRoot system ---
+//--- ANALYSIS system ---
 class AliLog ;
 
 class AliNeutralMesonSelection : public TObject {
@@ -44,8 +43,6 @@ class AliNeutralMesonSelection : public TObject {
   AliNeutralMesonSelection & operator = (const AliNeutralMesonSelection & g) ;//cpy assignment
   virtual ~AliNeutralMesonSelection() ; //virtual dtor
 
-  enum Type {kSelectPhiMinPt, kSelectPhiPtRatio, kNoSelectPhiPt};
-
   TList * GetCreateOutputObjects();
   
   Double_t GetAngleMaxParam(Int_t i) const {return fAngleMaxParam.At(i) ; }
@@ -55,25 +52,9 @@ class AliNeutralMesonSelection : public TObject {
   Double_t GetInvMassMinCut() const {return fInvMassMinCut ; }
   void SetInvMassCutRange(Double_t invmassmin, Double_t invmassmax)
   {fInvMassMaxCut =invmassmax;  fInvMassMinCut =invmassmin;}   
-  
-  Double_t GetDeltaPhiMaxCut() const {return fDeltaPhiMaxCut ; }
-  Double_t GetDeltaPhiMinCut() const {return fDeltaPhiMinCut ; }
-  void SetDeltaPhiCutRange(Double_t phimin, Double_t phimax)
-  {fDeltaPhiMaxCut =phimax;  fDeltaPhiMinCut =phimin;}
-
-  Double_t GetRatioMaxCut() const {return fRatioMaxCut ; }
-  Double_t GetRatioMinCut() const {return fRatioMinCut ; }
-  void SetRatioCutRange(Double_t ratiomin, Double_t ratiomax)
-  {fRatioMaxCut = ratiomax;  fRatioMinCut = ratiomin;}
-
-  Float_t    GetMinPt() const {return fMinPt ; }
-  void SetMinPt(Float_t pt){fMinPt = pt; };
 
   Double_t GetMass() const {return fM ; }
   void SetMass(Double_t m) { fM =m ; }
-  
-  Int_t GetPhiPtSelection() const { return fSelect ; }
-  void SetPhiPtSelection(Int_t ana ){  fSelect = ana ; }
 
   Bool_t AreNeutralMesonSelectionHistosKept() const { return fKeepNeutralMesonHistos ; }
   void KeepNeutralMesonSelectionHistos(Bool_t keep) { fKeepNeutralMesonHistos = keep ; }
@@ -82,29 +63,21 @@ class AliNeutralMesonSelection : public TObject {
   Bool_t IsAngleInWindow(const Float_t angle, const Float_t e) const ;
   void Print(const Option_t * opt) const;
 
-  Bool_t  CutPtPhi(Double_t ptg, Double_t phig, Double_t pt, Double_t phi) const ;
-  Bool_t  SelectPair(TParticle * photon, TLorentzVector particlei,  TLorentzVector particlej)  ;
+  Bool_t  SelectPair(TLorentzVector particlei,  TLorentzVector particlej)  ;
   
   private:
-  Int_t fSelect; //Pair selection depends on analysis
   Double_t fM ; //mass of the neutral meson
   Double_t   fInvMassMaxCut ;  // Invariant Mass cut maximum
   Double_t   fInvMassMinCut ;  // Invariant Masscut minimun
   TArrayD    fAngleMaxParam ; //Max opening angle selection parameters
   Double_t   fMinPt;       // Minimum pt 
-  Double_t   fDeltaPhiMaxCut ;      // Leading particle - gamma Ratio cut maximum
-  Double_t   fDeltaPhiMinCut ;      // Leading particle - gamma Ratio cut maximum
-  Double_t   fRatioMaxCut ;    // Leading particle/gamma Ratio cut maximum
-  Double_t   fRatioMinCut ;    // Leading particle/gamma Ratio cut minimum
   Bool_t  fKeepNeutralMesonHistos ; // Keep neutral meson selection histograms
 
   //Histograms
   TH2F * fhAnglePairNoCut  ;  //Aperture angle of decay photons, no cuts
-  TH2F * fhAnglePairCorrelationCut  ;  //Aperture angle of decay photons, cut on phi/pT correlation with prompt gamma
   TH2F * fhAnglePairOpeningAngleCut   ;  //Aperture angle of decay photons, cut on opening angle
   TH2F * fhAnglePairAllCut   ;  //Aperture angle of decay photons, all cuts
   TH2F * fhInvMassPairNoCut    ;  //Invariant mass of decay photons, no cuts
-  TH2F * fhInvMassPairCorrelationCut    ;   //Invariant mass of decay photons, cut on phi/pT correlation with prompt gamma
   TH2F * fhInvMassPairOpeningAngleCut  ;  //Invariant mass of decay photons, cut on opening angle
   TH2F * fhInvMassPairAllCut   ;  //Invariant mass of decay photons, all cuts