Corrected coding violations
authorgustavo <gustavo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 29 Oct 2007 13:49:42 +0000 (13:49 +0000)
committergustavo <gustavo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 29 Oct 2007 13:49:42 +0000 (13:49 +0000)
25 files changed:
PWG4/AliAnaCaloTriggerMC.cxx
PWG4/AliAnaGamma.cxx
PWG4/AliAnaGamma.h
PWG4/AliAnaGammaCorrelation.h
PWG4/AliAnaGammaDirect.cxx
PWG4/AliAnaGammaDirect.h
PWG4/AliAnaGammaHadron.h
PWG4/AliAnaGammaJetFinder.cxx
PWG4/AliAnaGammaJetFinder.h
PWG4/AliAnaGammaJetLeadCone.cxx
PWG4/AliAnaGammaJetLeadCone.h
PWG4/AliAnaGammaParton.h
PWG4/AliAnaGammaPhos.cxx
PWG4/AliAnalysisTaskGamma.cxx
PWG4/AliAnalysisTaskGamma.h
PWG4/AliGammaDataReader.cxx
PWG4/AliGammaDataReader.h
PWG4/AliGammaMCDataReader.cxx
PWG4/AliGammaMCDataReader.h
PWG4/AliGammaMCReader.cxx
PWG4/AliGammaMCReader.h
PWG4/AliGammaReader.cxx
PWG4/AliGammaReader.h
PWG4/AliNeutralMesonSelection.cxx
PWG4/AliNeutralMesonSelection.h

index d69b214..67796f6 100644 (file)
@@ -67,7 +67,7 @@ AliAnaCaloTriggerMC::AliAnaCaloTriggerMC(const char *name) :
 }
 //____________________________________________________________________________
 AliAnaCaloTriggerMC::AliAnaCaloTriggerMC(const AliAnaCaloTriggerMC & ct) : 
-  fChain(ct.fChain), fESD(ct.fESD),
+  AliAnalysisTask(ct),fChain(ct.fChain), fESD(ct.fESD),
   fOutputContainer(ct.fOutputContainer), fCalorimeter(ct. fCalorimeter),
   fNtTrigger22(ct.fNtTrigger22), fNtTriggerNN(ct.fNtTriggerNN)
 {
index e34ac63..ba18027 100644 (file)
  */
 
 //_________________________________________________________________________
-// Base class for prompt gamma and correlation analysis
+// Base class for gamma and correlation analysis
+// It is called by the task class AliAnalysisGammaTask and it connects the input (ESD/AOD/MonteCarlo)
+// got with AliGammaReader (produces TClonesArrays of TParticles), with the analysis classes 
+// AliAnaGammaDirect, AliAnaGammaCorrelation ....
+//
 //*-- Author: Gustavo Conesa (INFN-LNF)
 
 // --- ROOT system ---
@@ -40,6 +44,7 @@
 #include "AliGammaReader.h" 
 #include "AliAnaGammaDirect.h" 
 #include "AliAnaGammaCorrelation.h" 
+#include "AliAnaGammaSelection.h" 
 #include "AliNeutralMesonSelection.h"
 #include "AliAODCaloCluster.h"
 #include "AliAODEvent.h"
@@ -54,7 +59,7 @@ ClassImp(AliAnaGamma)
     TObject(),
     fOutputContainer(0x0), 
     fAnaType(0),  fCalorimeter(0), fData(0x0), fKine(0x0), 
-    fReader(0x0), fGammaDirect(0x0), fGammaCorrelation(0x0),
+    fReader(0x0), fGammaDirect(0x0), fGammaCorrelation(0x0), fGammaSelection(0x0),
     fNeutralMesonSelection(0x0), fAODclusters(0x0), fNAODclusters(0)
 {
   //Default Ctor
@@ -66,6 +71,8 @@ ClassImp(AliAnaGamma)
     fGammaDirect = new AliAnaGammaDirect();
   if(!fGammaCorrelation)
     fGammaCorrelation = new AliAnaGammaCorrelation();
+  if(!fGammaSelection)
+    fGammaSelection = new AliAnaGammaSelection();
   if(!fNeutralMesonSelection)
     fNeutralMesonSelection = new AliNeutralMesonSelection();
 
@@ -82,6 +89,7 @@ AliAnaGamma::AliAnaGamma(const AliAnaGamma & g) :
   fAnaType(g.fAnaType),  fCalorimeter(g.fCalorimeter), 
   fData(g.fData), fKine(g.fKine),fReader(g.fReader),
   fGammaDirect(g.fGammaDirect), fGammaCorrelation(g.fGammaCorrelation),
+  fGammaSelection(g.fGammaSelection),
   fNeutralMesonSelection(g.fNeutralMesonSelection),  
   fAODclusters(g. fAODclusters), fNAODclusters(g.fNAODclusters)
 {
@@ -105,6 +113,7 @@ AliAnaGamma & AliAnaGamma::operator = (const AliAnaGamma & source)
   fReader = source.fReader ;
   fGammaDirect = source.fGammaDirect ;
   fGammaCorrelation = source.fGammaCorrelation ;
+  fGammaSelection = source.fGammaSelection ;
   fNeutralMesonSelection = source.fNeutralMesonSelection ;
 
   return *this;
@@ -124,6 +133,7 @@ AliAnaGamma::~AliAnaGamma()
   delete fReader ;
   delete fGammaDirect ;
   delete fGammaCorrelation ;
+  delete fGammaSelection ;
   delete fNeutralMesonSelection ;
 
 }
@@ -139,14 +149,28 @@ void AliAnaGamma::Init()
   
   //Fill container with appropriate histograms
   
+  //Set selection  analysis histograms
+  TList * selectcontainer =  fGammaSelection->GetCreateOutputObjects(); 
+  for(Int_t i = 0; i < selectcontainer->GetEntries(); i++){
+    Bool_t  add = kTRUE ;
+    TString name = (selectcontainer->At(i))->GetName();   
+    if(!fReader->IsEMCALOn() && name.Contains("EMCAL")) add = kFALSE;
+    if(!fReader->IsPHOSOn() && name.Contains("PHOS"))   add = kFALSE;
+    if(!fReader->IsCTSOn() &&  !fGammaSelection->FillCTS() && name.Contains("CTS"))   add = kFALSE;
+    if(add) fOutputContainer->Add(selectcontainer->At(i)) ;
+  }  //Set selection  analysis histograms
+  
   //Set prompt photon analysis histograms
   TList * promptcontainer =  fGammaDirect->GetCreateOutputObjects(); 
   for(Int_t i = 0; i < promptcontainer->GetEntries(); i++)
     fOutputContainer->Add(promptcontainer->At(i)) ;
   
   //Check if selected options are correct or set them when necessary
-  if(fReader->GetDataType() == AliGammaReader::kMCData)
+  if(fReader->GetDataType() == AliGammaReader::kMCData){
     fGammaDirect->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
+    fGammaSelection->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
+  }
   
   if(fAnaType == kCorrelation){
     
@@ -238,7 +262,7 @@ void AliAnaGamma::Print(const Option_t * opt) const
 Bool_t AliAnaGamma::ProcessEvent(Long64_t entry){
 
   AliDebug(1,Form("Entry %d",entry));
-  cout<<"Event >>>>>>>>>>>>> "<<entry<<endl;
+  //cout<<"Event >>>>>>>>>>>>> "<<entry<<endl;
   if(!fOutputContainer)
     AliFatal("Histograms not initialized");
 
@@ -274,6 +298,12 @@ Bool_t AliAnaGamma::ProcessEvent(Long64_t entry){
   //Temporal solution, just for testing
   //FillAODs(plPHOS,plEMCAL);  
 
+  //Select particles to do the final analysis.
+  if(fReader->IsEMCALOn()) fGammaSelection->Selection("EMCAL",plEMCAL,plPrimEMCAL);
+  if(fReader->IsPHOSOn())  fGammaSelection->Selection("PHOS",plPHOS,plPrimPHOS);
+  if(fReader->IsCTSOn() &&  fGammaSelection->FillCTS()) fGammaSelection->Selection("CTS",plCTS,plPrimCTS);
+
+
   //Search highest energy prompt gamma in calorimeter
   if(fCalorimeter == "PHOS")
     MakeAnalysis(plPHOS, plEMCAL, plCTS, plParton, plPrimPHOS) ; 
index 235cdee..cd6e303 100644 (file)
  */
 
 //_________________________________________________________________________
-// Base class for prompt gamma and correlation analysis
+// Base class for gamma and correlation analysis
+// It is called by the task class AliAnalysisGammaTask and it connects the input (ESD/AOD/MonteCarlo)
+// got with AliGammaReader (produces TClonesArrays of TParticles), with the analysis classes 
+// AliAnaGammaDirect, AliAnaGammaCorrelation ....
+//
 //*-- Author: Gustavo Conesa (INFN-LNF)
 
 // --- ROOT system ---
@@ -35,7 +39,7 @@
 class AliGammaReader ;
 class AliAnaGammaDirect ;
 class AliAnaGammaCorrelation ;
-class AliAnaGammaJetLeadCone ;
+class AliAnaGammaSelection ;
 class AliNeutralMesonSelection ;
 class AliAODEvent;
 
@@ -49,28 +53,34 @@ public:
   AliAnaGamma & operator = (const AliAnaGamma & g) ;//cpy assignment
   virtual ~AliAnaGamma() ; //virtual dtor
 
-  enum anatype_t {kPrompt, kCorrelation};
+  enum Anatype {kPrompt, kCorrelation};
 
   //Setter and getters
   TList * GetOutputContainer()      const {return fOutputContainer ; }
 
-  Int_t GetAnalysisType(){  return fAnaType ; }
+  Int_t GetAnalysisType() const {return fAnaType ; }
   void SetAnalysisType(Int_t ana ){  fAnaType = ana ; }
 
-  TString GetCalorimeter() {return fCalorimeter ; }
+  TString GetCalorimeter() const {return fCalorimeter ; }
   void SetCalorimeter(TString calo) {if (calo == "PHOS" || calo == "EMCAL") fCalorimeter = calo ;
     else AliFatal("Wrong calorimeter name") ; }
 
-  TObject * GetData() {return fData ; }
-  TObject * GetKine() {return fKine ;}
+  TObject * GetData() const {return fData ; }
+  TObject * GetKine() const {return fKine ;}
   void SetData(TObject * data) {fData = data ; }
   void SetKine(TObject * kine) {fKine = kine ; }
 
-  AliGammaReader * GetReader() {return fReader ; }
+  AliGammaReader * GetReader() const {return fReader ; }
   void SetReader(AliGammaReader * reader) { fReader = reader ; }
 
+  AliAnaGammaDirect * GetGammaDirect() const { return fGammaDirect ; }
+  AliAnaGammaCorrelation * GetGammaCorrelation() const { return fGammaCorrelation  ;}
+  AliAnaGammaSelection * GetGammaSelection() const { return fGammaSelection ;}
+  AliNeutralMesonSelection * GetNeutralMesonSelection() const { return fNeutralMesonSelection  ; }
+
   void SetGammaDirect(AliAnaGammaDirect * dg) { fGammaDirect = dg ; }
   void SetGammaCorrelation(AliAnaGammaCorrelation * gc) { fGammaCorrelation = gc ;}
+  void SetGammaSelection(AliAnaGammaSelection * gs) { fGammaSelection = gs ;}
   void SetNeutralMesonSelection(AliNeutralMesonSelection * nms) { fNeutralMesonSelection = nms ; }
 
   //AOD stuff  
@@ -99,6 +109,7 @@ public:
   AliGammaReader *      fReader ; //! Pointer to reader 
   AliAnaGammaDirect *   fGammaDirect ; //! Pointer to prompt gamma algorithm 
   AliAnaGammaCorrelation *   fGammaCorrelation ; //! Pointer to gamma correlation algorithm
+  AliAnaGammaSelection *   fGammaSelection ; //! Pointer to gamma selection algorithm
   AliNeutralMesonSelection *  fNeutralMesonSelection ; //! Pointer to pair selection for pi0 identification.
   TClonesArray* fAODclusters;        //! reconstructed jets
   Int_t         fNAODclusters;       //! number of reconstructed jets
index dd9f6fb..3cae78a 100644 (file)
@@ -38,7 +38,7 @@ public:
   AliAnaGammaCorrelation & operator = (const AliAnaGammaCorrelation & g) ;//cpy assignment
   virtual ~AliAnaGammaCorrelation() ; //virtual dtor
 
-  enum corrtype_t {kParton, kHadron, kJetLeadCone, kJetFinder};
+  enum Corrtype {kParton, kHadron, kJetLeadCone, kJetFinder};
 
   //General methods
 
@@ -47,12 +47,12 @@ public:
   void SetNeutralMesonSelection(AliNeutralMesonSelection * nms) 
   { fNeutralMesonSelection = nms ; } 
 
-  TList * GetOutputContainer() {return fOutputContainer ;} 
+  TList * GetOutputContainer() const {return fOutputContainer ;} 
   void SetOutputContainer(TList * oc) {fOutputContainer = oc ;}  
 
   void InitParameters();
 
-  Int_t GetCorrelationType(){  return fCorrelationType ; }
+  Int_t GetCorrelationType() const {  return fCorrelationType ; }
   void SetCorrelationType(Int_t ana ){  fCorrelationType = ana ; }
 
   void Print(const Option_t * opt) const;
@@ -61,7 +61,7 @@ public:
   void SetJetsOnlyInCTS(Bool_t opt){fJetsOnlyInCTS = opt; }
 
   virtual TList * GetCreateOutputObjects() {return fOutputContainer ;}
-  virtual void MakeGammaCorrelation(TParticle * p,  TClonesArray * ob1, TClonesArray * ob2)  {;}
+  virtual void MakeGammaCorrelation(TParticle * ,  TClonesArray *, TClonesArray *)  {;}
 
   //Gamma hadron correlations methods: kHadron
   Float_t    GetMinPtHadron() const {return fMinPtHadron ; }
index 9024ef9..7eae44b 100644 (file)
@@ -54,7 +54,7 @@ ClassImp(AliAnaGammaDirect)
     TObject(),
     fMinGammaPt(0.),
     fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
-    fICMethod(0), fAnaMC(0),
+    fICMethod(0), fAnaMC(0), fIsolatePi0(0),
     fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0), 
     fntuplePrompt(0),
     //kSeveralIC
@@ -76,7 +76,8 @@ AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
   fPtThreshold(g.fPtThreshold),
   fPtSumThreshold(g.fPtSumThreshold), 
   fICMethod(g.fICMethod),
-  fAnaMC( g.fAnaMC),
+  fAnaMC( g.fAnaMC), 
+  fIsolatePi0(g.fIsolatePi0),
   fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),
   fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),    
   fntuplePrompt(g.fntuplePrompt),
@@ -114,6 +115,7 @@ AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & sou
   fPtSumThreshold = source.fPtSumThreshold ; 
   fICMethod = source.fICMethod ;
   fAnaMC = source.fAnaMC ;
+  fIsolatePi0 =  source.fIsolatePi0 ;
 
   fhNGamma = source.fhNGamma ; 
   fhPhiGamma = source.fhPhiGamma ;
@@ -193,7 +195,7 @@ TList *  AliAnaGammaDirect::GetCreateOutputObjects()
   outputContainer->Add(fhConeSumPt) ;
   
   //NTUPLE
-  fntuplePrompt = new TNtuple("ntuplePromptGamma", "Tree of prompt #gamma", "ptcluster:phicluster:etacluster:ptprimary:phiprimary:etaprimary:pdgprimary:statusprimary");
+  fntuplePrompt = new TNtuple("ntuplePromptGamma", "Tree of prompt #gamma", "ptcluster:phicluster:etacluster:pdg:status:ptprimary:phiprimary:etaprimary:pdgprimary:statusprimary");
   outputContainer->Add(fntuplePrompt) ;
 
   if(fICMethod == kSeveralIC){
@@ -241,14 +243,17 @@ void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plC
   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)){
+    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)
     {
@@ -269,6 +274,8 @@ void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plC
     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);
@@ -294,7 +301,8 @@ void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plC
     }
 
     //Fill ntuple with cluster / MC data
-    fntuplePrompt->Fill(ptcluster,phicluster,etacluster,ptprimary,phiprimary, etaprimary,pdgprimary,statusprimary);
+    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)) ;
@@ -314,11 +322,12 @@ void AliAnaGammaDirect::InitParameters()
 
   fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
   fAnaMC = kFALSE ;
+  fIsolatePi0 = kFALSE ;
  //-----------kSeveralIC-----------------
   fNCones           = 5 ; 
-  fNPtThres         = 5 ; 
+  fNPtThres         = 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[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
 
 }
 
@@ -402,7 +411,8 @@ void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArr
   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)){
+    if((particle->Pt() > fMinGammaPt) && (particle->Pt() > ptC) && 
+       (particle->GetPdgCode() == 22 ||  (fIsolatePi0 && particle->GetPdgCode() == 111))){
       index = ipr ;
       ptC = particle->Pt();
       pCandidate = particle ;
@@ -410,13 +420,14 @@ void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArr
     }
   }
   
+  //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[fNCones][fNPtThres];
+    Int_t ncone[10][10];//[fNCones][fNPtThres];
     Bool_t  icPtThres   = kFALSE;
     Bool_t  icPtSum     = kFALSE;
     
@@ -427,22 +438,22 @@ void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArr
        ncone[icone][ipt]=0;
        fPtThreshold=fPtThresholds[ipt] ;
        MakeIsolationCut(plCTS,plCalo, pCandidate, index,  
-                        ncone[icone][ipt]=0, icPtThres, icPtSum,coneptsum);
+                        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;
-       }
+//     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(), coneptsum,type,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
-
 }
 
 //__________________________________________________________________
index 4985f5d..0515353 100644 (file)
@@ -43,7 +43,7 @@ public:
   AliAnaGammaDirect & operator = (const AliAnaGammaDirect & g) ;//cpy assignment
   virtual ~AliAnaGammaDirect() ; //virtual dtor
   
-  enum anatype_t {kNoIC, kPtIC, kSumPtIC, kSeveralIC};
+  enum Anatype {kNoIC, kPtIC, kSumPtIC, kSeveralIC};
   
   Double_t  GetMinGammaPt()    const {return fMinGammaPt ; }
   Float_t     GetConeSize()          const {return fConeSize ; }
@@ -68,6 +68,7 @@ public:
   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 ; }
 
   Int_t    GetNCones()                  const {return fNCones ; }
   Int_t    GetNPtThresholds()                const {return fNPtThres ; }
@@ -94,7 +95,7 @@ public:
                                            // 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.
   //Histograms  
   TH1F * fhNGamma    ;  //Number of (isolated) gamma identified
   TH2F * fhPhiGamma  ; // Phi of identified gamma
index 5f8d567..b81d6f3 100644 (file)
@@ -44,15 +44,15 @@ public:
   private:
   
   //Histograms
-  TH2F * fhPhiCharged  ; 
-  TH2F * fhPhiNeutral   ; 
-  TH2F * fhEtaCharged  ; 
-  TH2F * fhEtaNeutral   ; 
-  TH2F * fhDeltaPhiGammaCharged  ;  
-  TH2F * fhDeltaPhiGammaNeutral   ; 
-  TH2F * fhDeltaEtaGammaCharged  ; 
-  TH2F * fhDeltaEtaGammaNeutral  ; 
-  TH2F * fhDeltaPhiChargedPt  ;
+  TH2F * fhPhiCharged  ; //Phi distribution of charged particles
+  TH2F * fhPhiNeutral   ;  //Phi distribution of neutral particles
+  TH2F * fhEtaCharged  ; //Eta distribution of charged particles
+  TH2F * fhEtaNeutral   ; //Eta distribution of neutral particles
+  TH2F * fhDeltaPhiGammaCharged  ;  //Difference of charged particle phi and prompt gamma phi as function of gamma pT
+  TH2F * fhDeltaPhiGammaNeutral   ;  //Difference of neutral particle phi and prompt gamma phi as function of gamma pT
+  TH2F * fhDeltaEtaGammaCharged  ;  //Difference of charged particle eta and prompt gamma eta as function of gamma pT
+  TH2F * fhDeltaEtaGammaNeutral  ;  //Difference of neutral particle eta and prompt gamma eta as function of gamma pT
+  TH2F * fhDeltaPhiChargedPt  ;  //Difference of charged particle eta and prompt gamma eta as function of charged pT
 
   TH2F * fhCorrelationGammaNeutral  ; //Neutral hadron correlation histogram 
   TH2F * fhCorrelationGammaCharged  ; //Charged hadron correlation histogram
index ae895c5..5da1b9c 100644 (file)
@@ -17,6 +17,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.2  2007/08/17 12:40:04  schutz
+ * New analysis classes by Gustavo Conesa
+ *
  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  * new analysis classes in the the new analysis framework
  *
@@ -160,5 +163,6 @@ void  AliAnaGammaJetFinder::MakeGammaCorrelation(TParticle * pGamma, TClonesArra
 {
   //Gamma -Jet  Correlation Analysis
   AliDebug(2, "Begin jet analysis");
+  cout<<pGamma<<" "<<pl<<endl;
   AliInfo("Not implemented");  
 } 
index 435eafb..def8b49 100644 (file)
@@ -41,10 +41,10 @@ class AliAnaGammaJetFinder : public AliAnaGammaCorrelation {
        
   private:
        
-       TH2F * fhDeltaEtaJet;
-       TH2F * fhDeltaPhiJet;
-       TH2F * fhDeltaPtJet;
-       TH2F * fhPtRatJet;
+       TH2F * fhDeltaEtaJet; // Difference of jet eta and prompt gamma eta as function of gamma pT
+       TH2F * fhDeltaPhiJet;  // Difference of jet phi and prompt gamma phi as function of gamma pT
+       TH2F * fhDeltaPtJet; // Difference of jet pT and prompt gamma pT as function of gamma pT
+       TH2F * fhPtRatJet; // Ratio of jet pT and prompt gamma pT as function of gamma pT
        
        ClassDef(AliAnaGammaJetFinder,1)
  } ;
index 1aa58b8..4867065 100644 (file)
@@ -17,6 +17,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.2  2007/08/17 12:40:04  schutz
+ * New analysis classes by Gustavo Conesa
+ *
  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  * new analysis classes in the the new analysis framework
  *
  */
 
 //_________________________________________________________________________
-// Class for the analysis of gamma-jet correlations, jet reconstructed in cone around leading
+// Class for the analysis of gamma-jet correlations:
+// 1)Take the prompt photon found with AliAnaGammaDirect,
+// 2) Search for the highest pt leading particle opposite to the photon within a phi, pt window
+// 3) Take all particles around leading in a cone R with pt larger than threshold and construct the jet
 //
 //  Class created from old AliPHOSGammaJet 
 //  (see AliRoot versions previous Release 4-09)
 
 
 // --- ROOT system ---
-#include "Riostream.h"
 
 //---- AliRoot system ----
-#include "AliLog.h"
 #include "AliNeutralMesonSelection.h"
 #include "AliAnaGammaJetLeadCone.h"  
+#include "AliLog.h"
 
 ClassImp(AliAnaGammaJetLeadCone)
 
index 4f0f268..9ae910a 100644 (file)
 
 //_________________________________________________________________________
 // Class that contains the algorithm for the reconstruction of jet, cone around leading particle
+// 1)Take the prompt photon found with AliAnaGammaDirect,
+// 2) Search for the highest pt leading particle opposite to the photon within a phi, pt window
+// 3) Take all particles around leading in a cone R with pt larger than threshold and construct the jet
+//
 //  Class created from old AliPHOSGammaJet
 //  (see AliRoot versions previous Release 4-09)
 //-- Author: Gustavo Conesa (INFN-LNF)
@@ -144,45 +148,50 @@ public:
   UInt_t       fSelect  ;   //kTRUE: Selects all jets, no limits.
  
   //Histograms
-
-  TH2F * fhPhiCharged  ; 
-  TH2F * fhPhiNeutral   ; 
-  TH2F * fhEtaCharged  ; 
-  TH2F * fhEtaNeutral   ; 
-  TH2F * fhDeltaPhiGammaCharged  ;  
-  TH2F * fhDeltaPhiGammaNeutral   ; 
-  TH2F * fhDeltaEtaGammaCharged  ; 
-  TH2F * fhDeltaEtaGammaNeutral  ; 
-
-  TH2F * fhAnglePairLeading  ; 
-  TH2F * fhInvMassPairLeading  ; 
-  TH2F * fhChargedRatio  ; 
-  TH2F * fhNeutralRatio   ; 
-  TH1F * fhNBkg   ; 
-  TH2F * fhNLeading  ; 
-  TH1F * fhNJet  ; 
-  TH2F * fhJetRatio  ; 
-  TH2F * fhJetPt   ; 
-  TH2F * fhBkgRatio   ; 
-  TH2F * fhBkgPt  ; 
-  TH2F * fhJetFragment  ; 
-  TH2F * fhBkgFragment  ; 
-  TH2F * fhJetPtDist  ; 
-  TH2F * fhBkgPtDist  ; 
-
-  TH2F * fhJetRatios[5][5];
-  TH2F * fhJetPts[5][5];
-  TH2F * fhBkgRatios[5][5];  
-  TH2F * fhBkgPts[5][5];
+  //Particle distributions
+  TH2F * fhPhiCharged  ; //Phi distribution of charged particles
+  TH2F * fhPhiNeutral   ;  //Phi distribution of neutral particles
+  TH2F * fhEtaCharged  ;  //Eta distribution of charged particles
+  TH2F * fhEtaNeutral   ;  //Eta distribution of neutral particles
+  //Leading particle distributions
+  TH2F * fhDeltaPhiGammaCharged  ;   //Difference of charged particle phi and prompt gamma phi as function of gamma pT
+  TH2F * fhDeltaPhiGammaNeutral   ;  //Difference of neutral particle phi and prompt gamma phi as function of gamma pT
+  TH2F * fhDeltaEtaGammaCharged  ;  //Difference of charged particle eta and prompt gamma eta as function of gamma pT
+  TH2F * fhDeltaEtaGammaNeutral  ;   //Difference of charged particle eta and prompt gamma eta as function of charged pT
+
+  TH2F * fhAnglePairLeading  ; //Aperture angle of decay photons of leading pi0
+  TH2F * fhInvMassPairLeading  ; //Invariant mass of decay photons of leading pi0
+  TH2F * fhChargedRatio  ; //Ratio of leading charge and prompt gamma
+  TH2F * fhNeutralRatio   ;  //Ratio of leading neutral and prompt gamma
+  TH1F * fhNBkg   ; //Bakground multiplicity
+  TH2F * fhNLeading  ; //Accepted leading particle pt distribution
+
+  //Jet distributions
+  //Fixed cone and pt threshold
+  TH1F * fhNJet  ; //Accepted reconstructed Jet pt distribution 
+  TH2F * fhJetRatio  ; //Ratio of pt jet and pt gamma
+  TH2F * fhJetPt   ; //reconstructed pt jet vs prompt pt gamma
+  TH2F * fhBkgRatio   ; //leading pt bakground / pt gamma 
+  TH2F * fhBkgPt  ; //leading pt bakground vs pt gamma 
+  TH2F * fhJetFragment  ; //Accepted reconstructed jet fragmentation function
+  TH2F * fhBkgFragment  ;  //Background "fragmentation function"
+  TH2F * fhJetPtDist  ; //Jet particle pt distribution
+  TH2F * fhBkgPtDist  ; //Background jet particle pt distribution
+
+  //Variable cone and pt threshold
+  TH2F * fhJetRatios[5][5]; //Ratio of pt jet and pt gamma
+  TH2F * fhJetPts[5][5]; //reconstructed pt jet vs prompt pt gamma
+  TH2F * fhBkgRatios[5][5];  //leading pt bakground / pt gamma 
+  TH2F * fhBkgPts[5][5]; //leading pt bakground vs pt gamma 
   
-  TH2F * fhNLeadings[5][5];
-  TH1F * fhNJets[5][5];
-  TH1F * fhNBkgs[5][5];
+  TH2F * fhNLeadings[5][5];  //Accepted leading particle pt distribution
+  TH1F * fhNJets[5][5]; //Accepted reconstructed Jet pt distribution 
+  TH1F * fhNBkgs[5][5]; //Bakground multiplicity
   
-  TH2F * fhJetFragments[5][5];
-  TH2F * fhBkgFragments[5][5];
-  TH2F * fhJetPtDists[5][5];
-  TH2F * fhBkgPtDists[5][5];
+  TH2F * fhJetFragments[5][5];//Accepted reconstructed jet fragmentation function
+  TH2F * fhBkgFragments[5][5];  //Background "fragmentation function"
+  TH2F * fhJetPtDists[5][5]; //Jet particle pt distribution
+  TH2F * fhBkgPtDists[5][5]; //Background jet particle pt distribution
   
 
   ClassDef(AliAnaGammaJetLeadCone,1)
index 8642ea9..e2342ff 100644 (file)
@@ -41,10 +41,10 @@ class AliAnaGammaParton : public AliAnaGammaCorrelation {
        
   private:
        
-       TH2F * fhDeltaEtaParton;
-       TH2F * fhDeltaPhiParton;
-       TH2F * fhDeltaPtParton;
-       TH2F * fhPtRatParton;
+       TH2F * fhDeltaEtaParton; //Difference of parton eta and prompt gamma eta
+       TH2F * fhDeltaPhiParton; //Difference of parton phi and prompt gamma phi
+       TH2F * fhDeltaPtParton; //Difference of parton pT and prompt gamma pT
+       TH2F * fhPtRatParton; //Ratio of parton pT and prompt gamma pT
        
        ClassDef(AliAnaGammaParton,1)
  } ;
index 9a09834..4fc33eb 100644 (file)
@@ -52,6 +52,7 @@ AliAnaGammaPhos::AliAnaGammaPhos() :
   fTreeA(0x0),
   fPhotonId(1.0),
   fOutputList(0x0), 
+  fhPHOSPos(0),
   fhPHOS(0),
   fhPHOSEnergy(0),
   fhPHOSDigits(0),
@@ -74,6 +75,7 @@ AliAnaGammaPhos::AliAnaGammaPhos(const char *name) :
   fTreeA(0x0),
   fPhotonId(1.0),
   fOutputList(0x0), 
+  fhPHOSPos(0),
   fhPHOS(0),
   fhPHOSEnergy(0),
   fhPHOSDigits(0),
index d680967..94337ce 100644 (file)
@@ -44,7 +44,8 @@ AliAnalysisTaskGamma::AliAnalysisTaskGamma():
     fESD(0x0),
     fAOD(0x0),
     fTreeG(0x0),
-    fOutputContainer(0x0)
+    fOutputContainer(0x0),
+    fConfigName(0)
 {
   // Default constructor
 }
@@ -57,7 +58,8 @@ AliAnalysisTaskGamma::AliAnalysisTaskGamma(const char* name):
     fESD(0x0),
     fAOD(0x0),
     fTreeG(0x0),
-    fOutputContainer(0x0)
+    fOutputContainer(0x0),
+    fConfigName("ConfigGammaAnalysis")
 {
   // Default constructor
  
@@ -106,7 +108,13 @@ void AliAnalysisTaskGamma::Init()
   AliDebug(1,"Begin");
   
   // Call configuration file
-  gROOT->LoadMacro("ConfigGammaAnalysis.C");
+
+  if(fConfigName == ""){
+    fConfigName="ConfigGammaAnalysis";
+  }
+  AliInfo(Form("### Configuration file is %s.C ###", fConfigName.Data()));
+  gROOT->LoadMacro(fConfigName+".C");
   fAna = (AliAnaGamma*) gInterpreter->ProcessLine("ConfigGammaAnalysis()");
   
   if(!fAna)
index 7fe978b..fb2bde4 100644 (file)
@@ -25,7 +25,10 @@ class AliAnalysisTaskGamma : public AliAnalysisTask
     virtual void LocalInit() {Init();}
     virtual void Exec(Option_t *option);
     virtual void Terminate(Option_t *option);
-   
+
+    void SetConfigFileName(TString name ) {fConfigName = name ; }
+    TString GetConfigFileName() const {return fConfigName ; }
+
  private:
 
     AliAnaGamma* fAna; //  Pointer to the jet finder 
@@ -34,6 +37,7 @@ class AliAnalysisTaskGamma : public AliAnalysisTask
     AliAODEvent*       fAOD;       //! AOD
     TTree*        fTreeG;     //  tree of prompt gamma, does nothing for the moment 
     TList * fOutputContainer ; // Histogram container
+    TString fConfigName ; //Configuration file name
 
     ClassDef(AliAnalysisTaskGamma, 1); // Analysis task for standard gamma correlation analysis
 };
index 3ba3e68..5dbdd83 100644 (file)
 
 
 // --- ROOT system ---
-#include "Riostream.h"
 #include <TParticle.h>
+#include <TFormula.h>
 
 //---- ANALYSIS system ----
 #include "AliGammaDataReader.h" 
-#include "AliLog.h"
 #include "AliESDEvent.h"
 #include "AliESDVertex.h"
 #include "AliESDCaloCluster.h"
+#include "AliLog.h"
 
 ClassImp(AliGammaDataReader)
 
@@ -88,8 +88,7 @@ void AliGammaDataReader::CreateParticleList(TObject * data, TObject *,
                                            TClonesArray * ){
   
   //Create a list of particles from the ESD. These particles have been measured 
-  //by the Central Tracking system (TPC+ITS), PHOS and EMCAL 
-  //Also create particle list with mothers.
+  //by the Central Tracking system (TPC+ITS+...), PHOS and EMCAL 
 
   AliESDEvent* esd = (AliESDEvent*) data;
 
@@ -101,189 +100,181 @@ void AliGammaDataReader::CreateParticleList(TObject * data, TObject *,
   Double_t v[3] ; //vertex ;
   esd->GetVertex()->GetXYZ(v) ; 
   
-  //########### PHOS ##############
-  
-  Int_t begphos = 0;  
-  Int_t endphos = esd->GetNumberOfCaloClusters() ;  
+  //########### CALORIMETERS ##############
+    
+  Int_t nCaloCluster = esd->GetNumberOfCaloClusters() ;  
   Int_t indexPH = plPHOS->GetEntries() ;
-  Int_t begem = 0;
-  Int_t endem = endphos;
-
-  AliDebug(3,Form("First Calorimeter particle %d, last particle %d", begphos,endphos));
+  Int_t indexEM = plEMCAL->GetEntries() ;
   
-  for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
-    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
 
-    //Check if it is PHOS cluster
-    if(!clus->IsEMCAL())
-      begem=npar+1;
-    else
-      continue;
-
-    AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
-    if(clus->GetTrackMatched()==-1){
-      //Create a TParticle to fill the particle list
-      TLorentzVector momentum ;
-      clus->GetMomentum(momentum, v);      
-      Double_t phi = momentum.Phi();
-      if(phi<0) phi+=TMath::TwoPi() ;
-      if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
-        phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
-       
-       pid=clus->GetPid();     
-       Int_t pdg = 22;
-
-       if(IsPHOSPIDOn()){
-         AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
-                         momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
-                         pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
+  for (npar =  0; npar <  nCaloCluster; npar++) {//////////////CaloCluster loop
+    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve cluster from esd
+    Int_t type = clus->GetClusterType();
+    
+    //########### PHOS ##############
+    if(fSwitchOnPHOS && type ==  AliESDCaloCluster::kPHOSCluster){
+      AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
+      
+      if(clus->GetTrackMatched()==-1){
+       TLorentzVector momentum ;
+       clus->GetMomentum(momentum, v);      
+       Double_t phi = momentum.Phi();
+       if(phi<0) phi+=TMath::TwoPi() ;
+       if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
+          phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
          
-         Float_t wPhoton =  fPHOSPhotonWeight;
-         Float_t wPi0 =  fPHOSPi0Weight;
+         pid=clus->GetPid();   
+         Int_t pdg = 22;
+           
+         if(IsPHOSPIDOn()){
+           AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
+                           momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
+                           pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
+           
+           Float_t wPhoton =  fPHOSPhotonWeight;
+           Float_t wPi0 =  fPHOSPi0Weight;
+           
+           if(fPHOSWeightFormula){
+             wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ;
+             wPi0 =    fPHOSPi0WeightFormula->Eval(momentum.E());
+           }
+           
+           if(pid[AliPID::kPhoton] > wPhoton) 
+             pdg = kPhoton ;
+           else if(pid[AliPID::kPi0] > wPi0) 
+             pdg = kPi0 ; 
+           else if(pid[AliPID::kElectron] > fPHOSElectronWeight)  
+             pdg = kElectron ;
+           else if(pid[AliPID::kEleCon] > fPHOSElectronWeight) 
+             pdg = kEleCon ;
+           else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight) 
+             pdg = kChargedHadron ;  
+           else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight) 
+             pdg = kNeutralHadron ; 
+           
+           else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
+                   pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
+             pdg = kChargedUnknown  ; 
+           else 
+             pdg = kNeutralUnknown ; 
+           //neutral cluster, unidentifed.
+         }
          
-         if(fPHOSWeightFormula){
-           wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ;
-           wPi0 =    fPHOSPi0WeightFormula->Eval(momentum.E());
+         if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array
+           TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
+                                                momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
+           
+           AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
+           
+           new((*plPHOS)[indexPH++])   TParticle(*particle) ;
          }
+         else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n", 
+                              pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));      
          
-         if(pid[AliPID::kPhoton] > wPhoton) 
-           pdg = kPhoton ;
-         else if(pid[AliPID::kPi0] > wPi0) 
-           pdg = kPi0 ; 
-         else if(pid[AliPID::kElectron] > fPHOSElectronWeight)  
-           pdg = kElectron ;
-         else if(pid[AliPID::kEleCon] > fPHOSElectronWeight) 
-           pdg = kEleCon ;
-         else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight) 
-           pdg = kChargedHadron ;  
-         else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight) 
-           pdg = kNeutralHadron ; 
+       }//pt, eta, phi cut
+       else    AliDebug(4,"Particle not added");
+      }//track-match?
+    }//PHOS cluster
+
+    //################ EMCAL ##############
+    else if(fSwitchOnEMCAL &&  type ==  AliESDCaloCluster::kEMCALClusterv1){
+      AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
+      
+      if(clus->GetTrackMatched()==-1 ){
+       TLorentzVector momentum ;
+       clus->GetMomentum(momentum, v); 
+       Double_t phi = momentum.Phi();
+       if(phi<0) phi+=TMath::TwoPi() ;
+       if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
+          phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
          
-         else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
-                 pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
-           pdg = kChargedUnknown  ; 
-         else 
-           pdg = kNeutralUnknown ; 
-         //neutral cluster, unidentifed.
-       }
-       
-       if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array
-         TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
-                                              momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
+         pid=clus->GetPid();   
+         Int_t pdg = 22;
          
-         AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
+         if(IsEMCALPIDOn()){
+           AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
+                           momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
+                           pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
+           
+           if(pid[AliPID::kPhoton] > fEMCALPhotonWeight) 
+             pdg = kPhoton ;
+           else if(pid[AliPID::kPi0] > fEMCALPi0Weight) 
+             pdg = kPi0 ; 
+           else if(pid[AliPID::kElectron] > fEMCALElectronWeight)  
+             pdg = kElectron ;
+           else if(pid[AliPID::kEleCon] > fEMCALElectronWeight) 
+             pdg = kEleCon ;
+           else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight) 
+             pdg = kChargedHadron ;  
+           else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight) 
+             pdg = kNeutralHadron ; 
+           else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
+                   pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
+             pdg = kChargedUnknown ; 
+           else 
+             pdg = kNeutralUnknown ;
+         }
+         
+         if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array
+           
+           TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
+                                                momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
+           AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
+           
+           new((*plEMCAL)[indexEM++])   TParticle(*particle) ;
+         }
+         else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f", 
+                              pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
          
-         new((*plPHOS)[indexPH++])   TParticle(*particle) ;
-       }
-       else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n", 
-                            pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));        
+       }//pt, phi, eta cut
+       else    AliDebug(4,"Particle not added");
+      }//track-matched
+    }//EMCAL cluster
 
-      }//pt, eta, phi cut
-      else     AliDebug(4,"Particle not added");
-    }//not charged
   }//cluster loop
-  
+
+
   //########### CTS (TPC+ITS) #####################
-  Int_t begtpc   = 0 ;  
-  Int_t endtpc   = esd->GetNumberOfTracks() ;
+  Int_t nTracks   = esd->GetNumberOfTracks() ;
   Int_t indexCh  = plCTS->GetEntries() ;
-  AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
   
-  for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
-    AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
-    
-    //We want tracks fitted in the detectors:
-    ULong_t status=AliESDtrack::kTPCrefit;
-    status|=AliESDtrack::kITSrefit;
+  if(fSwitchOnCTS){
+    AliDebug(3,Form("Number of tracks %d",nTracks));
     
-    //We want tracks whose PID bit is set:
-    //     ULong_t status =AliESDtrack::kITSpid;
-    //     status|=AliESDtrack::kTPCpid;
-    
-    if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
-      // Do something with the tracks which were successfully
-      // re-fitted 
-      Double_t en = 0; //track ->GetTPCsignal() ;
-      Double_t mom[3];
-      track->GetPxPyPz(mom) ;
-      Double_t px = mom[0];
-      Double_t py = mom[1];
-      Double_t pz = mom[2]; //Check with TPC people if this is correct.
-      Int_t pdg = 11; //Give any charged PDG code, in this case electron.
-      //I just want to tag the particle as charged
-       TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
-                                           px, py, pz, en, v[0], v[1], v[2], 0);
-  
-      //TParticle * particle = new TParticle() ;
-      //particle->SetMomentum(px,py,pz,en) ;
-       if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut)
-        new((*plCTS)[indexCh++])       TParticle(*particle) ;    
-    }
-  }
-  
-  //################ EMCAL ##############
-  
-  Int_t indexEM  = plEMCAL->GetEntries() ; 
-  
-  for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
-    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
-    Int_t clustertype= clus->GetClusterType();
-    AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
-
-    if(clustertype == AliESDCaloCluster::kEMCALClusterv1 && clus->GetTrackMatched()==-1 ){
+    for (npar =  0; npar <  nTracks; npar++) {////////////// track loop
+      AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
       
-      TLorentzVector momentum ;
-      clus->GetMomentum(momentum, v); 
-      Double_t phi = momentum.Phi();
-      if(phi<0) phi+=TMath::TwoPi() ;
-      if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
-       phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
+      //We want tracks fitted in the detectors:
+      ULong_t status=AliESDtrack::kTPCrefit;
+      status|=AliESDtrack::kITSrefit;
       
-       pid=clus->GetPid();     
-       Int_t pdg = 22;
-
-       if(IsEMCALPIDOn()){
-         AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
-                         momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
-                         pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
-         
-         if(pid[AliPID::kPhoton] > fEMCALPhotonWeight) 
-           pdg = kPhoton ;
-         else if(pid[AliPID::kPi0] > fEMCALPi0Weight) 
-           pdg = kPi0 ; 
-         else if(pid[AliPID::kElectron] > fEMCALElectronWeight)  
-           pdg = kElectron ;
-         else if(pid[AliPID::kEleCon] > fEMCALElectronWeight) 
-           pdg = kEleCon ;
-         else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight) 
-           pdg = kChargedHadron ;  
-         else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight) 
-           pdg = kNeutralHadron ; 
-         else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
-                 pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
-           pdg = kChargedUnknown ; 
-         else 
-           pdg = kNeutralUnknown ;
-       }
-       
-       if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array
-         
-         TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
-                                              momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
-         AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
-         
-         new((*plEMCAL)[indexEM++])   TParticle(*particle) ;
-       }
-       else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f", 
-                            pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
+      //We want tracks whose PID bit is set:
+      //     ULong_t status =AliESDtrack::kITSpid;
+      //     status|=AliESDtrack::kTPCpid;
+      
+      if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
+       // Do something with the tracks which were successfully
+       // re-fitted 
+       Double_t en = 0; //track ->GetTPCsignal() ;
+       Double_t mom[3];
+       track->GetPxPyPz(mom) ;
+       Double_t px = mom[0];
+       Double_t py = mom[1];
+       Double_t pz = mom[2]; //Check with TPC people if this is correct.
+       Int_t pdg = 11; //Give any charged PDG code, in this case electron.
+       //I just want to tag the particle as charged
+       TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
+                                            px, py, pz, en, v[0], v[1], v[2], 0);
        
-      }//pt, phi, eta cut
-      else     AliDebug(4,"Particle not added");
-    }//not charged, not pseudocluster
-  }//Cluster loop
+       //TParticle * particle = new TParticle() ;
+       //particle->SetMomentum(px,py,pz,en) ;
+       if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut)
+         new((*plCTS)[indexCh++])       TParticle(*particle) ;    
+      }// select track from refit
+    }//track loop
+  }//CTS
+  
+  AliDebug(3,Form("Particle lists filled, tracks  %d , clusters: EMCAL %d, PHOS %d", indexCh,indexEM,indexPH));
   
-  AliDebug(3,Form("Particle lists filled, PHOS clusters %d, EMCAL clusters %d, CTS tracks %d", 
-                 indexPH,indexEM,indexCh));
-
 }
 
index f4a62f5..19ee021 100644 (file)
 
 //*-- Author: Gustavo Conesa (INFN-LNF)
 
-// --- ROOT system ---
-#include <TParticle.h> 
-#include <TClonesArray.h> 
-#include "AliGammaReader.h" 
+// --- ROOT system --- 
 
-class AliESDEvent ;
+// --- AliRoot system ---
+#include "AliGammaReader.h" 
 
 class AliGammaDataReader : public AliGammaReader {
 
index 474ebd0..bb44f47 100644 (file)
 
 
 // --- ROOT system ---
-#include "Riostream.h"
+#include <TFormula.h>
 #include <TParticle.h>
-
 //---- ANALYSIS system ----
 #include "AliGammaMCDataReader.h" 
-#include "AliLog.h"
 #include "AliESDEvent.h"
 #include "AliESDVertex.h"
 #include "AliESDCaloCluster.h"
 #include "AliStack.h"
+#include "AliLog.h"
 
 ClassImp(AliGammaMCDataReader)
 
@@ -90,7 +90,8 @@ void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
                                              TClonesArray * plPrimPHOS){
   
   //Create a list of particles from the ESD. These particles have been measured 
-  //by the Central Tracking system (TPC+ITS), PHOS and EMCAL 
+  //by the Central Tracking system (TPC+ITS+...), PHOS and EMCAL 
+  //Also create particle list with mothers.
 
   AliESDEvent* esd = (AliESDEvent*) data;
   AliStack* stack = (AliStack*) kine;
@@ -103,221 +104,207 @@ void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
   Double_t v[3] ; //vertex ;
   esd->GetVertex()->GetXYZ(v) ; 
   
-  //########### PHOS ##############
-  
-  Int_t begphos = 0;  
-  Int_t endphos = esd->GetNumberOfCaloClusters() ;  
+  //########### CALORIMETERS ##############  
+  Int_t nCaloCluster = esd->GetNumberOfCaloClusters() ;  
   Int_t indexPH = plPHOS->GetEntries() ;
-  Int_t begem = 0;
-  Int_t endem = endphos;
-  
-  AliDebug(3,Form("First Calorimeter particle %d, last particle %d", begphos,endphos));
+  Int_t indexEM = plEMCAL->GetEntries() ;
   
-  for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
-    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
-
-    //Check if it is PHOS cluster
-    if(!clus->IsEMCAL())
-      begem=npar+1;
-    else
-      continue;
-
-    AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
-    if(clus->GetTrackMatched()==-1 ){
-      //Create a TParticle to fill the particle list
-      TLorentzVector momentum ;
-      clus->GetMomentum(momentum, v);
-      Double_t phi = momentum.Phi();
-      if(phi<0) phi+=TMath::TwoPi() ;
-      if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
-        phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
-       
-       pid=clus->GetPid();     
-       Int_t pdg = 22;
-
-       if(IsPHOSPIDOn()){
-         AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
-                         momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
-                         pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
+  for (npar =  0; npar <  nCaloCluster; npar++) {//////////////CaloCluster loop
+    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve cluster from esd
+    Int_t type = clus->GetClusterType();
+    
+    //########### PHOS ##############
+    if(fSwitchOnPHOS && type ==  AliESDCaloCluster::kPHOSCluster){
+      AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
+      
+      if(clus->GetTrackMatched()==-1){
+       TLorentzVector momentum ;
+       clus->GetMomentum(momentum, v);      
+       Double_t phi = momentum.Phi();
+       if(phi<0) phi+=TMath::TwoPi() ;
+       if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
+          phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
          
-         Float_t wPhoton =  fPHOSPhotonWeight;
-         Float_t wPi0 =  fPHOSPi0Weight;
+         pid=clus->GetPid();   
+         Int_t pdg = 22;
          
-         if(fPHOSWeightFormula){
-           wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ;
-           wPi0 =    fPHOSPi0WeightFormula->Eval(momentum.E());
+         if(IsPHOSPIDOn()){
+           AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
+                           momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
+                           pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
+           
+           Float_t wPhoton =  fPHOSPhotonWeight;
+           Float_t wPi0 =  fPHOSPi0Weight;
+           
+           if(fPHOSWeightFormula){
+             wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ;
+             wPi0 =    fPHOSPi0WeightFormula->Eval(momentum.E());
+           }
+           
+           if(pid[AliPID::kPhoton] > wPhoton) 
+             pdg = kPhoton ;
+           else if(pid[AliPID::kPi0] > wPi0) 
+             pdg = kPi0 ; 
+           else if(pid[AliPID::kElectron] > fPHOSElectronWeight)  
+             pdg = kElectron ;
+           else if(pid[AliPID::kEleCon] > fPHOSElectronWeight) 
+             pdg = kEleCon ;
+           else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight) 
+             pdg = kChargedHadron ;  
+           else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight) 
+             pdg = kNeutralHadron ; 
+           
+           else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
+                   pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
+             pdg = kChargedUnknown  ; 
+           else 
+             pdg = kNeutralUnknown ; 
+           //neutral cluster, unidentifed.
          }
          
-         if(pid[AliPID::kPhoton] > wPhoton) 
-           pdg = kPhoton ;
-         else if(pid[AliPID::kPi0] > wPi0) 
-           pdg = kPi0 ; 
-         else if(pid[AliPID::kElectron] > fPHOSElectronWeight)  
-           pdg = kElectron ;
-         else if(pid[AliPID::kEleCon] > fPHOSElectronWeight) 
-           pdg = kEleCon ;
-         else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight) 
-           pdg = kChargedHadron ;  
-         else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight) 
-           pdg = kNeutralHadron ; 
+         if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array
+           
+           //###############
+           //Check kinematics
+           //###############
+           Int_t label = clus->GetLabel();
+           TParticle * pmother = GetMotherParticle(label,stack, "PHOS",momentum);
+           
+           TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
+                                                momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
+           
+           AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
+           new((*plPHOS)[indexPH])   TParticle(*particle) ;
+           new((*plPrimPHOS)[indexPH])   TParticle(*pmother) ;
+           indexPH++;
+           
+         }
+         else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n", 
+                              pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));      
          
-         else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
-                 pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
-           pdg = kChargedUnknown  ; 
-         else 
-           pdg = kNeutralUnknown ; 
-         //neutral cluster, unidentifed.
-       }
-       
-       if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array
-
-         //###############
-         //Check kinematics
-         //###############
-         Int_t label = clus->GetLabel();
-         TParticle * pmother = GetMotherParticle(label,stack, "PHOS",momentum);
-
-         TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
-                                              momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
+       }//pt, eta, phi cut
+       else    AliDebug(4,"Particle not added");
+      }//track-match?
+    }//PHOS cluster
+
+    //################ EMCAL ##############
+    else if(fSwitchOnEMCAL &&  type ==  AliESDCaloCluster::kEMCALClusterv1){
+      AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
+      
+      if(clus->GetTrackMatched()==-1 ){
+       TLorentzVector momentum ;
+       clus->GetMomentum(momentum, v); 
+       Double_t phi = momentum.Phi();
+       if(phi<0) phi+=TMath::TwoPi() ;
+       if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
+          phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
          
-         AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
-         //printf("PHOS added: pdg %d, pt %f, phi %f, eta %f\n", pdg, particle->Pt(),particle->Phi(),particle->Eta());
-         new((*plPHOS)[indexPH])   TParticle(*particle) ;
-         new((*plPrimPHOS)[indexPH])   TParticle(*pmother) ;
-         indexPH++;
+         pid=clus->GetPid();   
+         Int_t pdg = 22;
          
-       }
-       else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n", pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));     
+         if(IsEMCALPIDOn()){
+           AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
+                           momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
+                           pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
+           
+           if(pid[AliPID::kPhoton] > fEMCALPhotonWeight) 
+             pdg = kPhoton ;
+           else if(pid[AliPID::kPi0] > fEMCALPi0Weight) 
+             pdg = kPi0 ; 
+           else if(pid[AliPID::kElectron] > fEMCALElectronWeight)  
+             pdg = kElectron ;
+           else if(pid[AliPID::kEleCon] > fEMCALElectronWeight) 
+             pdg = kEleCon ;
+           else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight) 
+             pdg = kChargedHadron ;  
+           else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight) 
+             pdg = kNeutralHadron ; 
+           else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
+                   pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
+             pdg = kChargedUnknown ; 
+           else 
+             pdg = kNeutralUnknown ;
+         }
+         
+         if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array
+           
+           //###############
+           //Check kinematics
+           //###############
+           Int_t label = clus->GetLabel();
+           TParticle * pmother = GetMotherParticle(label,stack, "EMCAL",momentum);
+           
+           TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
+                                                momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
+           AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
+           
+           new((*plEMCAL)[indexEM])   TParticle(*particle) ;
+           new((*plPrimEMCAL)[indexEM])   TParticle(*pmother) ;
+           indexEM++;
+         }
+         else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f", 
+                              pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
+         
+       }//pt, phi, eta cut
+       else    AliDebug(4,"Particle not added");
+      }//track-matched
+    }//EMCAL cluster
 
-      }//pt, eta, phi cut
-    }//not charged
   }//cluster loop
 
+
   //########### CTS (TPC+ITS) #####################
-  Int_t begtpc   = 0 ;  
-  Int_t endtpc   = esd->GetNumberOfTracks() ;
+  Int_t nTracks   = esd->GetNumberOfTracks() ;
   Int_t indexCh  = plCTS->GetEntries() ;
-  AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
-  
-  for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
-    AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
-    
-    //We want tracks fitted in the detectors:
-    ULong_t status=AliESDtrack::kTPCrefit;
-    status|=AliESDtrack::kITSrefit;
-    
-    //We want tracks whose PID bit is set:
-    //     ULong_t status =AliESDtrack::kITSpid;
-    //     status|=AliESDtrack::kTPCpid;
-    
-    if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
-      // Do something with the tracks which were successfully
-      // re-fitted 
-      Double_t en = 0; //track ->GetTPCsignal() ;
-      Double_t mom[3];
-      track->GetPxPyPz(mom) ;
-      Double_t px = mom[0];
-      Double_t py = mom[1];
-      Double_t pz = mom[2]; //Check with TPC people if this is correct.
-      Int_t pdg = 11; //Give any charged PDG code, in this case electron.
-      //I just want to tag the particle as charged
-
-      //Check kinematics
-      Int_t label = track->GetLabel();
-      TParticle * pmother = new TParticle();
-      if(label>=0){
-       pmother = stack->Particle(label);
-       //          AliDebug(4,("Primary Track: Pt %2.2f, mother: label %d, pdg  %d, status %d,  mother 2 %d, grandmother %d",
-       //                      TMath::Sqrt(px*px+py*py), label, pmother->GetPdgCode(),pmother->GetStatusCode(), pmother->GetMother(0), stack->GetPrimary(label)));
-      }
-      else AliDebug(3, Form("Negative Kinematic label of track:  %d",label));
-      
-      TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
-                                          px, py, pz, en, v[0], v[1], v[2], 0);
-      
-      if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut){
-       new((*plCTS)[indexCh])       TParticle(*particle) ;   
-       new((*plPrimCTS)[indexCh])       TParticle(*pmother) ;
-       indexCh++;    
-       
-      }
-    }
-  }
-  
-  //################ EMCAL ##############
-  
-  Int_t indexEM  = plEMCAL->GetEntries() ; 
-  
+
+  if(fSwitchOnCTS){
+    AliDebug(3,Form("Number of tracks %d",nTracks));
   
-  for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
-    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
-    Int_t clustertype= clus->GetClusterType();
-    if(clustertype == AliESDCaloCluster::kEMCALClusterv1 && clus->GetTrackMatched()==-1 ){
+    for (npar =  0; npar <  nTracks; npar++) {////////////// track loop
+      AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
       
-      TLorentzVector momentum ;
-      clus->GetMomentum(momentum, v);          
-      Double_t phi = momentum.Phi();
-      if(phi<0) phi+=TMath::TwoPi() ;
-      if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
-        phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
-       
-      pid=clus->GetPid();      
-      Int_t pdg = 22;
-
-      if(IsEMCALPIDOn()){
-       AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
-                       momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
-                       pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
-       
-       if(pid[AliPID::kPhoton] > fEMCALPhotonWeight) 
-           pdg = kPhoton ;
-       else if(pid[AliPID::kPi0] > fEMCALPi0Weight) 
-         pdg = kPi0 ; 
-       else if(pid[AliPID::kElectron] > fEMCALElectronWeight)  
-         pdg = kElectron ;
-       else if(pid[AliPID::kEleCon] > fEMCALElectronWeight) 
-         pdg = kEleCon ;
-       else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight) 
-         pdg = kChargedHadron ;  
-       else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight) 
-         pdg = kNeutralHadron ; 
-       else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
-               pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
-         pdg = kChargedUnknown ; 
-       else 
-         pdg = kNeutralUnknown ;
-      }
+      //We want tracks fitted in the detectors:
+      ULong_t status=AliESDtrack::kTPCrefit;
+      status|=AliESDtrack::kITSrefit;
+    
+      //We want tracks whose PID bit is set:
+      //     ULong_t status =AliESDtrack::kITSpid;
+      //     status|=AliESDtrack::kTPCpid;
       
-      if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array
+      if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
+       // Do something with the tracks which were successfully
+       // re-fitted 
+       Double_t en = 0; //track ->GetTPCsignal() ;
+       Double_t mom[3];
+       track->GetPxPyPz(mom) ;
+       Double_t px = mom[0];
+       Double_t py = mom[1];
+       Double_t pz = mom[2]; //Check with TPC people if this is correct.
+       Int_t pdg = 11; //Give any charged PDG code, in this case electron.
+       //I just want to tag the particle as charged
        
-       //###############
        //Check kinematics
-       //###############
-       Int_t label = clus->GetLabel();
-       TParticle * pmother = GetMotherParticle(label,stack, "EMCAL",momentum);
+       Int_t label = TMath::Abs(track->GetLabel());
+       TParticle * pmother = new TParticle();
+       pmother = stack->Particle(label);
        
        TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
-                                            momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
-       AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
+                                            px, py, pz, en, v[0], v[1], v[2], 0);
        
-       new((*plEMCAL)[indexEM])   TParticle(*particle) ;
-       new((*plPrimEMCAL)[indexEM])   TParticle(*pmother) ;
-       indexEM++;
-      }
-      else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f", 
-                          pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
-      
-      }//pt, phi, eta cut
-      else     AliDebug(4,"Particle not added");
-      
-    }//not charged, not pseudocluster
-  }//Cluster loop
-  
-  AliDebug(3,Form("Particle lists filled, PHOS clusters %d, EMCAL clusters %d, CTS tracks %d", 
-                 indexPH,indexEM,indexCh));
-  
-}
+       if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut){
+         new((*plCTS)[indexCh])       TParticle(*particle) ;   
+         new((*plPrimCTS)[indexCh])       TParticle(*pmother) ;
+         indexCh++;    
+         
+       }//kinematic selection
+      }//select track from refit
+    }//track loop    
+  }//CTS
 
+  AliDebug(3,Form("Particle lists filled, tracks  %d , clusters: EMCAL %d, PHOS %d", indexCh,indexEM,indexPH));
+
+}
 
 TParticle * AliGammaMCDataReader:: GetMotherParticle(Int_t label, AliStack *stack, TString calo,  TLorentzVector momentum)
 {
index f6c1c1d..88af3d2 100644 (file)
 //*-- Author: Gustavo Conesa (INFN-LNF)
 
 // --- ROOT system ---
-#include <TParticle.h> 
-#include <TClonesArray.h> 
-#include "AliGammaReader.h" 
 
-class AliESDEvent ;
+// --- AliRoot system ---
+#include "AliGammaReader.h" 
 
 class AliGammaMCDataReader : public AliGammaReader {
 
index 8195314..14633e4 100644 (file)
 #include <TH2.h>
 #include <TChain.h>
 #include <TRandom.h>
+#include <TArrayI.h>
+#include <TClonesArray.h>
 
 //---- ANALYSIS system ----
 #include "AliGammaMCReader.h" 
 #include "Riostream.h"
 #include "AliLog.h"
+#include "AliStack.h"
 
 ClassImp(AliGammaMCReader)
 
 //____________________________________________________________________________
 AliGammaMCReader::AliGammaMCReader() : 
-  AliGammaReader(), fDecayPi0(0), fCheckOverlapping(kFALSE)
+  AliGammaReader(), fDecayPi0(0), fCheckOverlapping(kFALSE),  fNeutralParticlesArray(0x0)
 {
   //Ctor
   
@@ -58,7 +61,8 @@ AliGammaMCReader::AliGammaMCReader() :
 
 //____________________________________________________________________________
 AliGammaMCReader::AliGammaMCReader(const AliGammaMCReader & g) :   
-  AliGammaReader(g), fDecayPi0(g.fDecayPi0), fCheckOverlapping(g.fCheckOverlapping)
+  AliGammaReader(g), fDecayPi0(g.fDecayPi0), fCheckOverlapping(g.fCheckOverlapping),
+  fNeutralParticlesArray(g.fNeutralParticlesArray?new TArrayI(*g.fNeutralParticlesArray):0x0)
 {
   // cpy ctor
 }
@@ -72,37 +76,54 @@ AliGammaMCReader & AliGammaMCReader::operator = (const AliGammaMCReader & source
 
   fDecayPi0 = source.fDecayPi0; 
   fCheckOverlapping = source.fCheckOverlapping;
-
+  delete fNeutralParticlesArray;
+  fNeutralParticlesArray = source.fNeutralParticlesArray?new TArrayI(*source.fNeutralParticlesArray):0x0;
   return *this;
 
 }
 
+//_________________________________
+AliGammaMCReader::~AliGammaMCReader() {
+  //Dtor
+
+  delete fNeutralParticlesArray ;
+
+}
+
 //___________________________________________________________________________
 void AliGammaMCReader::CaseDecayGamma(Int_t index, TParticle * particle, AliStack * sta,
                                  TClonesArray * plEMCAL, Int_t &indexEMCAL,
                                  TClonesArray * plPHOS, Int_t &indexPHOS){
   //In case pi0 are decayed by pythia, check if mother is pi0 and in such case look if 
   //there is overlapping. Send particle=pi0 if there is overlapping. 
-  
+
   TParticle * pmother =sta->Particle(particle->GetFirstMother());
   if(pmother->GetPdgCode() == 111 && pmother->GetNDaughters() == 2) {//Do not consider 3 particle decay case
     Int_t idaug0 = pmother->GetDaughter(0);
     Int_t idaug1 = pmother->GetDaughter(1);
     TParticle * pdaug0 = sta -> Particle(idaug0);
     TParticle * pdaug1 = sta -> Particle(idaug1);
-    
+
     if((index ==  idaug0 &&  pdaug0->Pt() > fNeutralPtCut) ||
-       (index ==  idaug1 &&  pdaug0->Pt() <= fNeutralPtCut))//Check decay when first daughter arrives, do nothing with second.
+       (index ==  idaug1 &&  pdaug0->Pt() <= fNeutralPtCut)){//Check decay when first daughter arrives, do nothing with second.
+
       FillListWithDecayGammaOrPi0(pmother, pdaug0, pdaug1, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
     
+    }
   }//mother is a pi0 with 2 daughters
   else{
-    if(IsInPHOS(particle->Phi(),particle->Eta()))
+    
+    if(IsInPHOS(particle->Phi(),particle->Eta())) {
+      TParticle * pmother =sta->Particle(particle->GetFirstMother());     
+      SetPhotonStatus(particle,pmother);     
       new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
-    else if(IsInEMCAL(particle->Phi(),particle->Eta()))
+    }
+    else if(IsInEMCAL(particle->Phi(),particle->Eta())){
+      TParticle * pmother =sta->Particle(particle->GetFirstMother());     
+      SetPhotonStatus(particle,pmother);
       new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
-  }
-
+    }
+  } 
 }
 
 //___________________________________________________________________________
@@ -158,6 +179,10 @@ void AliGammaMCReader::InitParameters()
 
   fDecayPi0 = kGeantDecay;
   fCheckOverlapping = kTRUE ;
+
+  Int_t pdgarray[]={12,14,16};// skip neutrinos
+  fNeutralParticlesArray = new TArrayI(3, pdgarray);
+
 }
 
 //____________________________________________________________________________
@@ -191,16 +216,17 @@ void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
       charge = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode())->Charge();
       
       //---------- Charged particles ----------------------
-      if((charge != 0) && (particle->Pt() > fChargedPtCut)){
+      if( fSwitchOnCTS && (charge != 0) && (particle->Pt() > fChargedPtCut)){
        //Particles in CTS acceptance
        if(TMath::Abs(particle->Eta())<fCTSEtaCut){  
          //Fill lists
          new((*plCh)[indexCh++])       TParticle(*particle) ;
        }
       }
+
       //-------------Neutral particles ----------------------
-      else if((charge == 0) && particle->Pt() > fNeutralPtCut &&  
-             TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos
+      else if((charge == 0) && particle->Pt() > fNeutralPtCut ){ //&&  
+       //TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos
        
        if(particle->GetPdgCode()!=111){
          //In case that we force PYTHIA to decay pi0, and we want to check the overlapping of 
@@ -209,11 +235,19 @@ void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
            CaseDecayGamma(iParticle,particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS); //If pythia decays pi0
          }
          else{
-           if(IsInPHOS(particle->Phi(),particle->Eta()))
-             new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
-           else if(IsInEMCAL(particle->Phi(),particle->Eta()))
-             new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
-         }
+           //Take out some particles like neutrinos
+           if(!SkipNeutralParticles(particle->GetPdgCode())){
+             TParticle * pmother =stack->Particle(particle->GetFirstMother());                
+             if(IsInPHOS(particle->Phi(),particle->Eta())){
+               if(particle->GetPdgCode()==22)SetPhotonStatus(particle,pmother);     
+               new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
+             }
+             else if(IsInEMCAL(particle->Phi(),particle->Eta())){
+               if(particle->GetPdgCode()==22) SetPhotonStatus(particle,pmother); 
+               new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
+             }
+           }//skip neutral particles
+         }//Case kDecayGamma
        }//no pi0
        else{
          if(fDecayPi0 == kNoDecay){//keep the pi0 do not decay
@@ -221,7 +255,7 @@ void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
              new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
            else if(IsInEMCAL(particle->Phi(),particle->Eta()))
              new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
-         }
+         }//nodecay
          else if(fDecayPi0 == kDecay)
            CasePi0Decay(particle,plEMCAL,indexEMCAL,plPHOS, indexPHOS);
          else if(fDecayPi0 == kGeantDecay)
@@ -239,31 +273,31 @@ void AliGammaMCReader::FillListWithDecayGammaOrPi0(TParticle * pPi0, TParticle *
                                   TClonesArray * plPHOS, Int_t &indexPHOS){
 
   //Check if decay gamma overlapp in calorimeter, in such case keep the pi0, if not keep both photons.
-  
   Bool_t  overlap = kFALSE ;
   TLorentzVector lv1 , lv2 ;
   pdaug0->Momentum(lv1);
   pdaug1->Momentum(lv2);
   Double_t angle = lv1.Angle(lv2.Vect());
-  
-  if(fCheckOverlapping){//Check if decay products overlapp
-    if(IsInEMCAL(pPi0->Phi(), pPi0->Eta())){
+  if(IsInEMCAL(pPi0->Phi(), pPi0->Eta()))
+    {
       if (angle < fEMCALMinAngle){
        new((*plEMCAL)[indexEMCAL++])       TParticle(*pPi0) ;
        overlap = kTRUE;
       }
     }
-    else if(IsInPHOS(pPi0->Phi(), pPi0->Eta())){
-      if (angle < fPHOSMinAngle){
-       new((*plPHOS)[indexPHOS++])       TParticle(*pPi0) ;
-       overlap = kTRUE;
-      }
+
+  else if(IsInPHOS(pPi0->Phi(), pPi0->Eta())){
+    if (angle < fPHOSMinAngle){
+      new((*plPHOS)[indexPHOS++])       TParticle(*pPi0) ;
+      overlap = kTRUE;
     }
   }//fCheckOverlapping
-  
+
   //Fill with gammas if not overlapp
   if(!overlap){
     if(pdaug0->GetPdgCode() == 22 || TMath::Abs(pdaug0->GetPdgCode() ) == 11 ){
+      pdaug0->SetStatusCode(kPi0DecayPhoton);
       if(IsInEMCAL(pdaug0->Phi(),pdaug0->Eta()) &&  pdaug0->Pt() > fNeutralPtCut)
        new((*plEMCAL)[indexEMCAL++])       TParticle(*pdaug0) ;
       else if(IsInPHOS(pdaug0->Phi(),pdaug0->Eta()) &&  pdaug0->Pt() > fNeutralPtCut)
@@ -271,36 +305,43 @@ void AliGammaMCReader::FillListWithDecayGammaOrPi0(TParticle * pPi0, TParticle *
     }
     
     if(pdaug1->GetPdgCode() == 22 || TMath::Abs(pdaug1->GetPdgCode() ) == 11 ){
+      pdaug1->SetStatusCode(kPi0DecayPhoton);
       if(IsInEMCAL(pdaug1->Phi(),pdaug1->Eta()) &&  pdaug1->Pt() > fNeutralPtCut)
        new((*plEMCAL)[indexEMCAL++])       TParticle(*pdaug1) ;
       else if(IsInPHOS(pdaug1->Phi(),pdaug1->Eta()) &&  pdaug1->Pt() > fNeutralPtCut)
        new((*plPHOS)[indexPHOS++])       TParticle(*pdaug1) ;
     }
   }// overlap?
+
  }
 
 
  
 //___________________________________________________________________________
-Bool_t  AliGammaMCReader::IsInEMCAL(Double_t phi, Double_t eta){
+Bool_t  AliGammaMCReader::IsInEMCAL(Double_t phi, Double_t eta) const {
   //Check if particle is in EMCAL acceptance
-  if(phi<0)
-     phi+=TMath::TwoPi();
-     if( phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] && 
-       TMath::Abs(eta)<fEMCALEtaCut) return kTRUE ;
-  else  return kFALSE;     
-  
+  if(fSwitchOnEMCAL){
+    if(phi<0)
+      phi+=TMath::TwoPi();
+    if( phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] && 
+       TMath::Abs(eta)<fEMCALEtaCut) return kTRUE ;
+    else  return kFALSE;     
+  }
   return kFALSE ;
 }
 
 //___________________________________________________________________________
-Bool_t  AliGammaMCReader::IsInPHOS(Double_t phi, Double_t eta){
-  //Check if particle is in EMCAL acceptance
-  if(phi<0)
-    phi+=TMath::TwoPi();
-  if( phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] && 
-      TMath::Abs(eta)<fPHOSEtaCut) return kTRUE ;
-  else  return kFALSE;
+Bool_t  AliGammaMCReader::IsInPHOS(Double_t phi, Double_t eta) const {
+  //Check if particle is in PHOS acceptance
+  
+  if(fSwitchOnPHOS){
+    if(phi<0)
+      phi+=TMath::TwoPi();
+    if( phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] && 
+       TMath::Abs(eta)<fPHOSEtaCut) return kTRUE ;
+    else  return kFALSE;
+  }
   
   return kFALSE ;
 }
@@ -363,4 +404,33 @@ void AliGammaMCReader::Print(const Option_t * opt) const
 
 
 
+//________________________________________________________________
+void AliGammaMCReader::SetPhotonStatus(TParticle* pphoton, TParticle* pmother)
+{
 
+  //Check the origin of the photon and tag it  as decay from pi0, from eta, from other mesons, or prompt or fragmentation.
+  
+  if(pmother->GetStatusCode() != 21 ){
+    if(pmother->GetStatusCode() ==11){
+      if(pmother->GetPdgCode()==111) pphoton->SetStatusCode(kPi0DecayPhoton);//decay gamma from pi0
+      if(pmother->GetPdgCode()==221) pphoton->SetStatusCode(kEtaDecayPhoton);//decay gamma from eta
+      else pphoton->SetStatusCode(kOtherDecayPhoton);// decay gamma from other pphotons        
+    }
+    else        pphoton->SetStatusCode(kUnknown);
+  }
+  else if(pmother->GetPdgCode()==22)   pphoton->SetStatusCode(kPromptPhoton);//Prompt photon
+  else pphoton->SetStatusCode(kFragmentPhoton); //Fragmentation photon
+  
+}
+
+//________________________________________________________________
+Bool_t AliGammaMCReader::SkipNeutralParticles(Int_t pdg) const {
+  //Check if pdg is equal to one of the neutral particles list
+  //These particles will be skipped from analysis.
+
+  for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
+    if(TMath::Abs(pdg) ==  fNeutralParticlesArray->At(i)) return kTRUE ;
+  
+  return kFALSE; 
+  
+}
index cfbade3..a0c2306 100644 (file)
 //*-- Author: Gustavo Conesa (INFN-LNF)
 
 // --- ROOT system ---
-#include <TParticle.h> 
-#include <TClonesArray.h> 
-#include "AliStack.h"
+class TArrayI ;
+
+// --- AliRoot system ---
 #include "AliGammaReader.h" 
  
-class TH2F ; 
-
 class AliGammaMCReader : public AliGammaReader {
 
 public: 
@@ -32,23 +30,29 @@ public:
   AliGammaMCReader() ; // ctor
   AliGammaMCReader(const AliGammaMCReader & g) ; // cpy ctor
   AliGammaMCReader & operator = (const AliGammaMCReader & g) ;//cpy assignment
-  virtual ~AliGammaMCReader() {;} //virtual dtor
+  virtual ~AliGammaMCReader() ;//virtual dtor
 
-  enum decay_t {kNoDecay, kGeantDecay, kDecay, kDecayGamma};
+  enum DecayType {kNoDecay, kGeantDecay, kDecay, kDecayGamma};
  
   void InitParameters();
 
-  Bool_t  IsInEMCAL(Double_t phi, Double_t eta) ;
-  Bool_t  IsInPHOS(Double_t phi, Double_t eta) ;
+  Bool_t  IsInEMCAL(Double_t phi, Double_t eta) const ;
+  Bool_t  IsInPHOS(Double_t phi, Double_t eta) const ;
 
   Int_t    GetDecayPi0Flag() const {return fDecayPi0 ; }
 
-  void Print(const Option_t * opt)const;
+  void Print(const Option_t * opt) const;
   
   void SetDecayPi0Flag(Int_t d){ fDecayPi0 = d ; }
 
   void SetCheckOverlapping(Bool_t check){fCheckOverlapping = check ;}
-  Bool_t IsCheckOverlappingOn() {return fCheckOverlapping ;}
+  Bool_t IsCheckOverlappingOn() const {return fCheckOverlapping ;}
+
+  void SetPhotonStatus(TParticle* pphoton, TParticle* mother);
+
+  void AddNeutralParticlesArray(TArrayI & array)  
+  { fNeutralParticlesArray   = new TArrayI(array) ; }
+  TArrayI * GetNeutralParticlesArray() const   {return  fNeutralParticlesArray;}
 
   private:
   
@@ -72,12 +76,14 @@ public:
                                   TClonesArray * plPHOS, Int_t &indexPHOS);  
   void MakePi0Decay(TLorentzVector &p0, TLorentzVector &p1, TLorentzVector &p2);//, Double_t &angle);
   
-  
+  Bool_t SkipNeutralParticles(Int_t pdg) const ;
   private:
   
-  Int_t      fDecayPi0; //Decay Pi0.
+  Int_t      fDecayPi0; //Options to study pi0 decay
   Bool_t   fCheckOverlapping; // if True, check if gammas from decay overlapp in calorimeters.
-  
+  TArrayI * fNeutralParticlesArray ; //Do not keep neutral particles of the list.
   ClassDef(AliGammaMCReader,1)
 } ;
 
index 8a42521..7f635ff 100644 (file)
  */
 
 //_________________________________________________________________________
-// Base class for reading data in order to do prompt gamma correlations
+// Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and 
+// Central Barrel Tracking detectors.
+// Not all MC particles/tracks/clusters are kept, some kinematical restrictions are done.
+// Mother class of : AliGammaDataReader: Fills ESD data in 3 TClonesArrays (PHOS, EMCAL, CTS)
+//                 : AliGammaMCReader: Fills Kinematics data in 3 TClonesArrays (PHOS, EMCAL, CTS)
+//                 : AliGammaMCDataReader: Fills ESD data in 3 TClonesArrays (PHOS, EMCAL, CTS) 
+//                             and MC data in other 3 TClonesArray
 //*-- Author: Gustavo Conesa (LNF-INFN) 
 //////////////////////////////////////////////////////////////////////////////
 
 
 // --- ROOT system ---
+#include <TFormula.h>
+#include <TMath.h>
 
 //---- ANALYSIS system ----
-#include "Riostream.h"
-#include "AliLog.h"
 #include "AliGammaReader.h"
 
 ClassImp(AliGammaReader)
@@ -45,6 +51,7 @@ ClassImp(AliGammaReader)
 //____________________________________________________________________________
 AliGammaReader::AliGammaReader() : 
   TObject(), fDataType(0),
+  fSwitchOnEMCAL(0),  fSwitchOnPHOS(0),  fSwitchOnCTS(0),
   fCTSEtaCut(0.), fEMCALEtaCut(0.), fPHOSEtaCut(0.),
   fNeutralPtCut(0.), fChargedPtCut(0.),
   fEMCALIPDistance(0.),   fPHOSIPDistance(0.), 
@@ -70,6 +77,7 @@ AliGammaReader::AliGammaReader() :
 //____________________________________________________________________________
 AliGammaReader::AliGammaReader(const AliGammaReader & g) :   
   TObject(g), fDataType(g.fDataType),
+  fSwitchOnEMCAL(g.fSwitchOnEMCAL),  fSwitchOnPHOS(g.fSwitchOnPHOS),  fSwitchOnCTS(g.fSwitchOnCTS),
   fCTSEtaCut(g.fCTSEtaCut),  fEMCALEtaCut(g.fEMCALEtaCut),  fPHOSEtaCut(g.fPHOSEtaCut),
   fNeutralPtCut(g.fNeutralPtCut), fChargedPtCut(g.fChargedPtCut),
   fEMCALIPDistance(g.fEMCALIPDistance),  fPHOSIPDistance(g.fPHOSIPDistance),  
@@ -106,6 +114,11 @@ AliGammaReader & AliGammaReader::operator = (const AliGammaReader & source)
   if(&source == this) return *this;
 
   fDataType = source.fDataType ;
+
+  fSwitchOnEMCAL = source.fSwitchOnEMCAL;  
+  fSwitchOnPHOS = source.fSwitchOnPHOS;
+  fSwitchOnCTS = source.fSwitchOnCTS;
+
   fCTSEtaCut = source.fCTSEtaCut;  
   fEMCALEtaCut = source.fEMCALEtaCut;  
   fPHOSEtaCut = source.fPHOSEtaCut;
@@ -152,6 +165,11 @@ void AliGammaReader::InitParameters()
  
   //Initialize the parameters of the analysis.
   fDataType = kData ;
+
+  fSwitchOnEMCAL = kTRUE ;  
+  fSwitchOnPHOS = kTRUE ;
+  fSwitchOnCTS = kTRUE ;
+
   fCTSEtaCut         = 0.7 ;  
   fEMCALEtaCut         = 0.7 ;  
   fPHOSEtaCut         = 0.12 ;
@@ -203,6 +221,11 @@ void AliGammaReader::Print(const Option_t * opt) const
 
   Info("Print", "%s %s", GetName(), GetTitle() ) ;
   printf("Data type           : %d\n", fDataType) ;
+
+  printf(" EMCAL on?          : %d\n", fSwitchOnEMCAL) ;
+  printf(" PHOS on?          : %d\n", fSwitchOnPHOS) ;
+  printf(" CTS on?          : %d\n", fSwitchOnCTS) ;
+
   printf("CTS Eta cut           : %f\n", fCTSEtaCut) ;
   printf("EMCAL Eta cut           : %f\n", fEMCALEtaCut) ;
   printf("PHOS Eta cut           : %f\n", fPHOSEtaCut) ;
index 08fcb66..b8d6f15 100644 (file)
  */
 
 //_________________________________________________________________________
-// Base class for reading data in order to do prompt gamma correlations 
+// Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and 
+// Central Barrel Tracking detectors.
+// Not all MC particles/tracks/clusters are kept, some kinematical restrictions are done.
+// Mother class of : AliGammaDataReader: Fills ESD data in 3 TClonesArrays (PHOS, EMCAL, CTS)
+//                 : AliGammaMCReader: Fills Kinematics data in 3 TClonesArrays (PHOS, EMCAL, CTS)
+//                 : AliGammaMCDataReader: Fills ESD data in 3 TClonesArrays (PHOS, EMCAL, CTS) 
+//                             and MC data in other 3 TClonesArray
 //*-- Author: Gustavo Conesa (INFN-LNF)
 
 // --- ROOT system ---
-#include <TParticle.h> 
-#include <TFormula.h>
-#include <TClonesArray.h> 
-#include "AliStack.h"
 #include "TObject.h" 
-class AliESD ; 
 
-class TH2F ; 
+class TClonesArray ; 
+class TFormula ;
+class TParticle ; 
+class Riostream ;
+//--- AliRoot system ---
+
+class AliStack ;
+class AliESDEvent ; 
+class AliLog ;
 
 class AliGammaReader : public TObject {
 
@@ -52,42 +60,51 @@ public:
     kChargedUnknown=321
   };
 
-  enum datatype_t {kData, kMC, kMCData};
+  enum PhotonStatusType {
+    kPromptPhoton=2,
+    kFragmentPhoton=3,
+    kPi0DecayPhoton=4, 
+    kEtaDecayPhoton=5, 
+    kOtherDecayPhoton=6,
+    kUnknown=7
+  };
+
+  enum Datatype {kData, kMC, kMCData};
   
   void InitParameters();
 
-  Int_t GetDataType(){ return fDataType ; }
+  Int_t GetDataType() const { return fDataType ; }
   void SetDataType(Int_t data ){fDataType = data ; }
 
-  virtual Float_t   GetCTSEtaCut() const {return fCTSEtaCut ; }
-  virtual Float_t   GetEMCALEtaCut() const {return fEMCALEtaCut ; }
-  virtual Float_t   GetPHOSEtaCut() const {return fPHOSEtaCut ; }
+  virtual Float_t  GetCTSEtaCut() const {return fCTSEtaCut ; }
+  virtual Float_t  GetEMCALEtaCut() const {return fEMCALEtaCut ; }
+  virtual Float_t  GetPHOSEtaCut() const {return fPHOSEtaCut ; }
   virtual Float_t  GetPhiEMCALCut(Int_t i) { return  fPhiEMCALCut[i]  ; }
   virtual Float_t  GetPhiPHOSCut(Int_t i) { return  fPhiPHOSCut[i]  ; }
   virtual Float_t  GetNeutralPtCut()    {  return fNeutralPtCut  ; }
   virtual Float_t  GetChargedPtCut()  {  return fChargedPtCut  ; }
 
-  virtual Float_t  GetEMCALIPDistance()  {  return fEMCALIPDistance ; }
-  virtual Float_t  GetPHOSIPDistance()  {  return fPHOSIPDistance ; }
-  virtual Float_t  GetEMCALMinAngle()  {  return fEMCALMinAngle ; }
-  virtual Float_t  GetPHOSMinAngle()  {  return fPHOSMinAngle ; }
+  virtual Float_t  GetEMCALIPDistance()  const {  return fEMCALIPDistance ; }
+  virtual Float_t  GetPHOSIPDistance()  const {  return fPHOSIPDistance ; }
+  virtual Float_t  GetEMCALMinAngle()  const {  return fEMCALMinAngle ; }
+  virtual Float_t  GetPHOSMinAngle()  const {  return fPHOSMinAngle ; }
 
   virtual Bool_t   IsEMCALPIDOn() const {return fEMCALPID ; }
   virtual Bool_t   IsPHOSPIDOn() const {return fPHOSPID ; }
-  virtual Float_t  GetEMCALPhotonWeight() { return  fEMCALPhotonWeight  ; }
-  virtual Float_t  GetEMCALPi0Weight()    {  return fEMCALPi0Weight  ; }
-  virtual Float_t  GetEMCALElectronWeight() { return  fEMCALElectronWeight  ; }
-  virtual Float_t  GetEMCALChargeWeight()    {  return fEMCALChargeWeight  ; }
-  virtual Float_t  GetEMCALNeutralWeight()    {  return fEMCALNeutralWeight  ; }
-  virtual Float_t  GetPHOSPhotonWeight()  {  return fPHOSPhotonWeight  ; }
-  virtual Float_t  GetPHOSPi0Weight()  {  return fPHOSPi0Weight  ; }
-  virtual Float_t  GetPHOSElectronWeight()  {  return fPHOSElectronWeight  ; }
-  virtual Float_t  GetPHOSChargeWeight()  {  return fPHOSChargeWeight  ; }
-  virtual Float_t  GetPHOSNeutralWeight()  {  return fPHOSNeutralWeight  ; }
-
-  virtual Bool_t  IsPHOSPIDWeightFormulaOn()  {  return fPHOSWeightFormula  ; } 
-  virtual TFormula * GetPHOSPhotonWeightFormula()    {  return fPHOSPhotonWeightFormula  ; } 
-  virtual TFormula * GetPHOSPi0WeightFormula()   {  return fPHOSPi0WeightFormula  ; }
+  virtual Float_t  GetEMCALPhotonWeight() const  { return  fEMCALPhotonWeight  ; }
+  virtual Float_t  GetEMCALPi0Weight() const     {  return fEMCALPi0Weight  ; }
+  virtual Float_t  GetEMCALElectronWeight() const  { return  fEMCALElectronWeight  ; }
+  virtual Float_t  GetEMCALChargeWeight() const     {  return fEMCALChargeWeight  ; }
+  virtual Float_t  GetEMCALNeutralWeight() const     {  return fEMCALNeutralWeight  ; }
+  virtual Float_t  GetPHOSPhotonWeight() const   {  return fPHOSPhotonWeight  ; }
+  virtual Float_t  GetPHOSPi0Weight() const   {  return fPHOSPi0Weight  ; }
+  virtual Float_t  GetPHOSElectronWeight() const   {  return fPHOSElectronWeight  ; }
+  virtual Float_t  GetPHOSChargeWeight() const   {  return fPHOSChargeWeight  ; }
+  virtual Float_t  GetPHOSNeutralWeight() const   {  return fPHOSNeutralWeight  ; }
+
+  virtual Bool_t  IsPHOSPIDWeightFormulaOn() const   {  return fPHOSWeightFormula  ; } 
+  virtual TFormula * GetPHOSPhotonWeightFormula() const     {  return fPHOSPhotonWeightFormula  ; } 
+  virtual TFormula * GetPHOSPi0WeightFormula() const    {  return fPHOSPi0WeightFormula  ; }
 
   virtual void Print(const Option_t * opt)const;
   
@@ -123,24 +140,42 @@ public:
   virtual void SetPHOSPhotonWeightFormula(TFormula * photon)    {  fPHOSPhotonWeightFormula  = photon; } 
   virtual void SetPHOSPi0WeightFormula(TFormula * pi0)   {  fPHOSPi0WeightFormula  = pi0; }
 
-  virtual void CreateParticleList(TObject* data, TObject * data2, 
-                                 TClonesArray * plCh, TClonesArray * plEMCAL, TClonesArray * plPHOS,     
-                                 TClonesArray * parton, TClonesArray * plPrimEMCAL, TClonesArray * plPrimPHOS) {;}
+  virtual Bool_t IsEMCALOn() const  { return fSwitchOnEMCAL ; }
+  virtual Bool_t IsPHOSOn() const  { return fSwitchOnPHOS ; }
+  virtual Bool_t IsCTSOn() const  { return fSwitchOnCTS ; }
+  
+  virtual void SwitchOnEMCAL(Bool_t sw) {fSwitchOnEMCAL = sw ; }
+  virtual void SwitchOnPHOS(Bool_t sw) {fSwitchOnPHOS = sw ; }
+  virtual void SwitchOnCTS(Bool_t sw) {fSwitchOnCTS = sw ; }
+
+  virtual void CreateParticleList(TObject* , TObject * , 
+                                 TClonesArray * , TClonesArray * , TClonesArray *,  
+                                 TClonesArray * , TClonesArray * , TClonesArray * ) {;}
+  
  protected:
-  Int_t        fDataType ;
+  
+  Int_t        fDataType ; //Select MC:Kinematics, Data:ESD/AOD, MCData:Both
+  
+  Bool_t     fSwitchOnEMCAL ; //Do not consider EMCAL neutral particles/clusters
+  Bool_t     fSwitchOnPHOS ;//Do not consider PHOS neutral particles/clusters
+  Bool_t     fSwitchOnCTS ;//Do not consider tracks/ charged particles
+
+  //Kinematical cuts
   Float_t      fCTSEtaCut ;//CTS  pseudorapidity acceptance
   Float_t      fEMCALEtaCut ;//EMCAL pseudorapidity acceptance
   Float_t      fPHOSEtaCut ;//PHOS pseudorapidity acceptance
   Float_t      fPhiEMCALCut[2]; //EMCAL phi acceptance 
   Float_t      fPhiPHOSCut[2];  //PHOS phi acceptance
-  Float_t      fNeutralPtCut; //
-  Float_t      fChargedPtCut;  // 
+  Float_t      fNeutralPtCut; //pT Threshold on neutral particles
+  Float_t      fChargedPtCut;  // pT  Threshold on charged particles
 
+  //Overlapping 
   Float_t      fEMCALIPDistance; //Calorimeter IP distance.
   Float_t      fPHOSIPDistance; //Calorimeter IP distance
   Float_t      fEMCALMinAngle; //Gamma decay minimum aperture angle for overlapping.
   Float_t      fPHOSMinAngle; //Gamma decay minimum aperture angle for overlapping.
 
+  //PID
   Bool_t       fEMCALPID ;//Fill EMCAL particle lists with particles with corresponding pid
   Bool_t       fPHOSPID;  //Fill PHOS particle lists with particles with corresponding pid
   Float_t      fEMCALPhotonWeight; //Bayesian PID weight for photons in EMCAL 
index dfdcfa9..e773e3a 100644 (file)
@@ -17,6 +17,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.2  2007/08/17 12:40:04  schutz
+ * New analysis classes by Gustavo Conesa
+ *
  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  * new analysis classes in the the new analysis framework
  *
@@ -25,6 +28,9 @@
 
 //_________________________________________________________________________
 // 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 <TLorentzVector.h>
 #include <TH2.h>
 #include <TList.h>
-
+#include <TArrayD.h>
 //---- AliRoot system ----
 #include "AliNeutralMesonSelection.h" 
-#include "Riostream.h"
 #include "AliLog.h"
 
 ClassImp(AliNeutralMesonSelection)
@@ -240,7 +246,7 @@ void AliNeutralMesonSelection::InitParameters()
 }
 
 //__________________________________________________________________________-
-Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) {
+Bool_t AliNeutralMesonSelection::IsAngleInWindow(const Float_t angle,const Float_t e) const {
   //Check if the opening angle of the candidate pairs is inside 
   //our selection windowd
 
@@ -259,7 +265,7 @@ 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)  
+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 ;
index 73421fe..4767e8a 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<TObject.h>
-#include <TArrayD.h>
+#include<TArrayD.h>
 
 class TLorentzVector ;
 class TParticle ;
 class TList ;
 class TH2F ;
+class Riostream ;
+
+//--- AliRoot system ---
+class AliLog ;
 
 class AliNeutralMesonSelection : public TObject {
 
@@ -37,7 +44,7 @@ class AliNeutralMesonSelection : public TObject {
   AliNeutralMesonSelection & operator = (const AliNeutralMesonSelection & g) ;//cpy assignment
   virtual ~AliNeutralMesonSelection() ; //virtual dtor
 
-  enum type_t {kSelectPhiMinPt, kSelectPhiPtRatio, kNoSelectPhiPt};
+  enum Type {kSelectPhiMinPt, kSelectPhiPtRatio, kNoSelectPhiPt};
 
   TList * GetCreateOutputObjects();
   
@@ -65,17 +72,17 @@ class AliNeutralMesonSelection : public TObject {
   Double_t GetMass() const {return fM ; }
   void SetMass(Double_t m) { fM =m ; }
   
-  Int_t GetPhiPtSelection(){  return fSelect ; }
+  Int_t GetPhiPtSelection() const { return fSelect ; }
   void SetPhiPtSelection(Int_t ana ){  fSelect = ana ; }
 
-  Bool_t AreNeutralMesonSelectionHistosKept() { return fKeepNeutralMesonHistos ; }
+  Bool_t AreNeutralMesonSelectionHistosKept() const { return fKeepNeutralMesonHistos ; }
   void KeepNeutralMesonSelectionHistos(Bool_t keep) { fKeepNeutralMesonHistos = keep ; }
 
   void InitParameters();       
-  Bool_t IsAngleInWindow(const Float_t angle, const Float_t e);
+  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) ;
+  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)  ;
   
   private:
@@ -85,21 +92,21 @@ class AliNeutralMesonSelection : public TObject {
   Double_t   fInvMassMinCut ;  // Invariant Masscut minimun
   TArrayD    fAngleMaxParam ; //Max opening angle selection parameters
   Double_t   fMinPt;       // Minimum pt 
-  Double_t   fDeltaPhiMaxCut ;      // 
-  Double_t   fDeltaPhiMinCut ;      // 
+  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  ; 
-  TH2F * fhAnglePairCorrelationCut  ; 
-  TH2F * fhAnglePairOpeningAngleCut   ; 
-  TH2F * fhAnglePairAllCut   ; 
-  TH2F * fhInvMassPairNoCut    ; 
-  TH2F * fhInvMassPairCorrelationCut    ;
-  TH2F * fhInvMassPairOpeningAngleCut  ; 
-  TH2F * fhInvMassPairAllCut   ; 
+  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
   
   ClassDef(AliNeutralMesonSelection,1)