Made ready for the analysis train
authorschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Oct 2007 17:25:12 +0000 (17:25 +0000)
committerschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Oct 2007 17:25:12 +0000 (17:25 +0000)
31 files changed:
PWG4/AliAnaCaloTrigger.cxx
PWG4/AliAnaCaloTrigger.h
PWG4/AliAnaCaloTriggerMC.cxx
PWG4/AliAnaCaloTriggerMC.h
PWG4/AliAnaGamma.cxx
PWG4/AliAnaGamma.h
PWG4/AliAnaGammaCorrelation.h
PWG4/AliAnaGammaDirect.cxx
PWG4/AliAnaGammaDirect.h
PWG4/AliAnaGammaHadron.cxx
PWG4/AliAnaGammaHadron.h
PWG4/AliAnaGammaJetFinder.h
PWG4/AliAnaGammaJetLeadCone.h
PWG4/AliAnaGammaParton.h
PWG4/AliAnaGammaPhos.cxx
PWG4/AliAnaGammaPhos.h
PWG4/AliAnaScale.cxx
PWG4/AliAnaScale.h
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.h
PWG4/ConfigGammaAnalysis.C
PWG4/ana.C

index e230d5f..b33f0d1 100644 (file)
@@ -226,7 +226,7 @@ void AliAnaCaloTrigger::Exec(Option_t *)
        phiMax = vpos.Phi() ; 
       }
 
-      Float_t * pid = cluster->GetPid() ;
+      Double_t * pid = cluster->GetPid() ;
       
       if(pid[AliPID::kPhoton] > 0.9) {
        if ( cluEnergy > phEnMax) { 
index fb6b112..163a412 100644 (file)
@@ -46,6 +46,6 @@ private:
   TNtuple * fNtTrigger22 ;
   TNtuple * fNtTriggerNN ;
 
-  ClassDef(AliAnaCaloTrigger, 0); // a PHOS photon analysis task 
+  ClassDef(AliAnaCaloTrigger, 1); // a photon analysis task 
 };
 #endif // ALIANACALOTRIGGER_H
index e29b9e4..d69b214 100644 (file)
@@ -27,6 +27,7 @@
 
 #include "AliAnalysisManager.h"
 #include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
 #include "AliAnaCaloTriggerMC.h" 
 #include "AliESDEvent.h" 
 #include "AliESDCaloCluster.h"
@@ -169,7 +170,7 @@ void AliAnaCaloTriggerMC::Exec(Option_t *)
     ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
   
   if(mctruth)
-    stack = mctruth->Stack();
+    stack = mctruth->MCEvent()->Stack();
   
   if (!stack) {
     AliError("Stack not found") ; 
@@ -245,7 +246,7 @@ void AliAnaCaloTriggerMC::Exec(Option_t *)
          labelmax = cluster->GetLabel();
        }
        
-       Float_t * pid = cluster->GetPid() ;
+       Double_t * pid = cluster->GetPid() ;
        
        if(pid[AliPID::kPhoton] > 0.9) {
          if ( cluEnergy > phEnMax) { 
index cd1bdd5..7912467 100644 (file)
@@ -45,6 +45,6 @@ private:
   TNtuple * fNtTrigger22 ;
   TNtuple * fNtTriggerNN ;
 
-  ClassDef(AliAnaCaloTriggerMC, 0); // a PHOS photon analysis task 
+  ClassDef(AliAnaCaloTriggerMC, 1); // a photon analysis task 
 };
 #endif // ALIANACALOTRIGGERMC_H
index 68d5179..e34ac63 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
  *
@@ -38,6 +41,8 @@
 #include "AliAnaGammaDirect.h" 
 #include "AliAnaGammaCorrelation.h" 
 #include "AliNeutralMesonSelection.h"
+#include "AliAODCaloCluster.h"
+#include "AliAODEvent.h"
 #include "Riostream.h"
 #include "AliLog.h"
 
@@ -50,7 +55,7 @@ ClassImp(AliAnaGamma)
     fOutputContainer(0x0), 
     fAnaType(0),  fCalorimeter(0), fData(0x0), fKine(0x0), 
     fReader(0x0), fGammaDirect(0x0), fGammaCorrelation(0x0),
-    fNeutralMesonSelection(0x0)
+    fNeutralMesonSelection(0x0), fAODclusters(0x0), fNAODclusters(0)
 {
   //Default Ctor
 
@@ -64,6 +69,8 @@ ClassImp(AliAnaGamma)
   if(!fNeutralMesonSelection)
     fNeutralMesonSelection = new AliNeutralMesonSelection();
 
+  fAODclusters = 0;
+
   InitParameters();
   
 }
@@ -75,7 +82,8 @@ 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),
-  fNeutralMesonSelection(g.fNeutralMesonSelection)
+  fNeutralMesonSelection(g.fNeutralMesonSelection),  
+  fAODclusters(g. fAODclusters), fNAODclusters(g.fNAODclusters)
 {
   // cpy ctor
   
@@ -136,6 +144,10 @@ void AliAnaGamma::Init()
   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)
+    fGammaDirect->SetMC();//Only useful with AliGammaMCDataReader, by default kFALSE
+  
   if(fAnaType == kCorrelation){
     
     //Check if selected options are correct
@@ -226,7 +238,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;
   if(!fOutputContainer)
     AliFatal("Histograms not initialized");
 
@@ -235,6 +247,9 @@ Bool_t AliAnaGamma::ProcessEvent(Long64_t entry){
   TClonesArray * plCTS      = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
   TClonesArray * plEMCAL    = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter (EMCAL)
   TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured  Gamma calorimeter
+  TClonesArray * plPrimCTS      = new TClonesArray("TParticle",1000); // primary tracks
+  TClonesArray * plPrimEMCAL    = new TClonesArray("TParticle",1000);   // primary emcal clusters
+  TClonesArray * plPrimPHOS     = new TClonesArray("TParticle",1000);  // primary phos clusters
   TClonesArray * plParton   = new TClonesArray("TParticle",1000);  // All partons
   //Fill lists with photons, neutral particles and charged particles
   //look for the highest energy photon in the event inside fCalorimeter
@@ -242,24 +257,28 @@ Bool_t AliAnaGamma::ProcessEvent(Long64_t entry){
   //Fill particle lists 
   if(fReader->GetDataType() == AliGammaReader::kData){
     AliDebug(1,"Data analysis");
-    fReader->CreateParticleList(fData, NULL,plCTS,plEMCAL,plPHOS,NULL); 
+    fReader->CreateParticleList(fData, NULL,plCTS,plEMCAL,plPHOS,NULL,NULL,NULL); 
   }
   else if( fReader->GetDataType()== AliGammaReader::kMC){
     AliDebug(1,"Kinematics analysis");
-    fReader->CreateParticleList(fKine, NULL,plCTS,plEMCAL,plPHOS,plParton); 
+    fReader->CreateParticleList(fKine, NULL,plCTS,plEMCAL,plPHOS,plParton,NULL,NULL); 
   }
   else if(fReader->GetDataType() == AliGammaReader::kMCData) {
    AliDebug(1,"Data + Kinematics analysis");
-   fReader->CreateParticleList(fData, fKine,plCTS,plEMCAL,plPHOS,NULL); 
+   fReader->CreateParticleList(fData, fKine,plCTS,plEMCAL,plPHOS,plPrimCTS,plPrimEMCAL,plPrimPHOS); 
   }
   else
     AliError("Option not implemented");
-  
+
+  //Fill AOD with calorimeter
+  //Temporal solution, just for testing
+  //FillAODs(plPHOS,plEMCAL);  
+
   //Search highest energy prompt gamma in calorimeter
   if(fCalorimeter == "PHOS")
-    MakeAnalysis(plPHOS, plEMCAL, plCTS, plParton) ; 
+    MakeAnalysis(plPHOS, plEMCAL, plCTS, plParton, plPrimPHOS) ; 
   else if (fCalorimeter == "EMCAL")
-    MakeAnalysis(plEMCAL, plPHOS, plCTS,plParton) ; 
+    MakeAnalysis(plEMCAL, plPHOS, plCTS,plParton, plPrimEMCAL) ; 
   else
     AliFatal("Wrong calorimeter name");
 
@@ -267,18 +286,24 @@ Bool_t AliAnaGamma::ProcessEvent(Long64_t entry){
   plEMCAL->Clear() ;
   plPHOS->Clear() ;
   plParton->Clear() ;
+  plPrimCTS->Clear() ;
+  plPrimEMCAL->Clear() ;
+  plPrimPHOS->Clear() ;
 
   delete plCTS ;
   delete plPHOS ;
   delete plEMCAL ;
   delete plParton ;
+  delete plPrimCTS ;
+  delete plPrimPHOS ;
+  delete plPrimEMCAL ;
 
   return kTRUE;
 
 }
 
 //____________________________________________________________________________
-void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plParton)  {
+void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plParton, TClonesArray * plPrimCalo)  {
   
   TParticle * pGamma = new TParticle ;
   Bool_t isInCalo = kFALSE ;
@@ -302,7 +327,7 @@ void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClon
            
          default :
            {
-             fGammaDirect->GetPromptGamma(plCalo, plCTS,pGamma,isInCalo);
+             fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
              if(!isInCalo)
                AliDebug(1,"Prompt gamma not found");
            }
@@ -316,7 +341,7 @@ void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClon
       {
        AliDebug(1,"kCorrelation analysis");
        //Find prompt photon    
-       fGammaDirect->GetPromptGamma(plCalo, plCTS,pGamma,isInCalo);
+       fGammaDirect->GetPromptGamma(plCalo, plCTS,plPrimCalo, pGamma,isInCalo);
 
        if(isInCalo){//If prompt photon found, do correlation
          
@@ -360,3 +385,44 @@ void AliAnaGamma::MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClon
   delete pGamma ; 
   
 }
+
+//____________________________________________________
+void AliAnaGamma::AddCluster(AliAODCaloCluster p)
+{
+// Add new jet to the list
+  cout<<"AOD list pointer "<<fAODclusters<<" nAODclusters "<<fNAODclusters<<endl;
+  cout<<"list entries "<<fAODclusters->GetEntries()<<endl;
+  new ((*fAODclusters)[fNAODclusters++]) AliAODCaloCluster(p);
+}
+
+//___________________________________________________
+void AliAnaGamma::ConnectAOD(AliAODEvent* aod)
+{
+// Connect to the AOD
+  fAODclusters = aod->GetCaloClusters();
+}
+
+//____________________________________________________________________________
+void AliAnaGamma::FillAODs(TClonesArray * plPHOS, TClonesArray * plEMCAL){
+  //Fill AOD caloClusters
+  //Temporal method, just for testing AOD creation
+  
+  //Fill AOD with PHOS clusters
+  Int_t nphos =  plPHOS->GetEntries() ;
+  cout<<"PHOS entries "<<nphos<<endl;
+  for(Int_t ipr = 0;ipr < nphos ; ipr ++ ){
+    TParticle * particle = dynamic_cast<TParticle *>(plPHOS->At(ipr)) ;
+    Float_t pos []= {0,0,0};
+    AliAODCaloCluster phos(ipr,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kPHOSNeutral,0);
+    AddCluster(phos);
+  }
+  
+  //Fill AOD with EMCAL clusters
+  cout<<"EMCAL entries "<<plEMCAL->GetEntries()<<endl;
+  for(Int_t ipr = 0;ipr < plEMCAL->GetEntries() ; ipr ++ ){
+    TParticle * particle = dynamic_cast<TParticle *>(plEMCAL->At(ipr)) ;
+    Float_t pos []= {0,0,0};
+    AliAODCaloCluster emcal(ipr+nphos,0,0x0,particle->Energy(),pos,0x0,AliAODCluster::kEMCALClusterv1,0);
+    AddCluster(emcal);
+  }
+}
index db34234..235cdee 100644 (file)
@@ -7,6 +7,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
  *
 #include <TH2F.h>
 #include<TObject.h>
 #include <TTree.h>
+
+// --- AliRoot system ---
 #include <AliLog.h>
+#include "AliAODCaloCluster.h"
 
 class AliGammaReader ;
 class AliAnaGammaDirect ;
 class AliAnaGammaCorrelation ;
 class AliAnaGammaJetLeadCone ;
 class AliNeutralMesonSelection ;
+class AliAODEvent;
 
 // --- AliRoot
 class AliAnaGamma : public TObject {
@@ -44,21 +51,12 @@ public:
 
   enum anatype_t {kPrompt, kCorrelation};
 
-  //General methods
+  //Setter and getters
   TList * GetOutputContainer()      const {return fOutputContainer ; }
-  
-  void Init();
-  void InitParameters();
 
   Int_t GetAnalysisType(){  return fAnaType ; }
   void SetAnalysisType(Int_t ana ){  fAnaType = ana ; }
 
-  void Print(const Option_t * opt) const;
-
-  void MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray *plParton)  ;  
-  Bool_t ProcessEvent(Long64_t entry) ;
-  //TTree * MakeTreeG(TString name) ;
-
   TString GetCalorimeter() {return fCalorimeter ; }
   void SetCalorimeter(TString calo) {if (calo == "PHOS" || calo == "EMCAL") fCalorimeter = calo ;
     else AliFatal("Wrong calorimeter name") ; }
@@ -74,7 +72,22 @@ public:
   void SetGammaDirect(AliAnaGammaDirect * dg) { fGammaDirect = dg ; }
   void SetGammaCorrelation(AliAnaGammaCorrelation * gc) { fGammaCorrelation = gc ;}
   void SetNeutralMesonSelection(AliNeutralMesonSelection * nms) { fNeutralMesonSelection = nms ; }
-  
+
+  //AOD stuff  
+  void   AddCluster(AliAODCaloCluster p);
+  void   ConnectAOD(AliAODEvent* aod);
+  void   FillAODs(TClonesArray * plPHOS, TClonesArray * plEMCAL);
+
+  //Others
+  void Init();
+  void InitParameters();
+
+  void MakeAnalysis(TClonesArray * plCalo, TClonesArray * plNe, TClonesArray * plCTS, TClonesArray *plParton, TClonesArray * plPrimCalo)  ;  
+
+  void Print(const Option_t * opt) const;
+
+  Bool_t ProcessEvent(Long64_t entry) ;
+
  private:
   
   //General Data members
@@ -87,8 +100,10 @@ public:
   AliAnaGammaDirect *   fGammaDirect ; //! Pointer to prompt gamma algorithm 
   AliAnaGammaCorrelation *   fGammaCorrelation ; //! Pointer to gamma correlation algorithm
   AliNeutralMesonSelection *  fNeutralMesonSelection ; //! Pointer to pair selection for pi0 identification.
-  
-  ClassDef(AliAnaGamma,0)
+  TClonesArray* fAODclusters;        //! reconstructed jets
+  Int_t         fNAODclusters;       //! number of reconstructed jets
+
+  ClassDef(AliAnaGamma,1)
 } ;
  
 
index baf83e1..dd9f6fb 100644 (file)
@@ -7,6 +7,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
  *
@@ -90,7 +93,7 @@ public:
   Double_t   fRatioMaxCut ;    // Leading particle/gamma Ratio cut maximum (kLeadJetCone)
   Double_t   fRatioMinCut ;    // Leading particle/gamma Ratio cut minimum (kLeadJetCone)
   
-  ClassDef(AliAnaGammaCorrelation,0)
+  ClassDef(AliAnaGammaCorrelation,1)
 } ;
  
 
index 6611cff..9024ef9 100644 (file)
@@ -17,6 +17,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.6  2007/08/17 12:40:04  schutz
+ * New analysis classes by Gustavo Conesa
+ *
  * Revision 1.4.4.4  2007/07/26 10:32:09  schutz
  * new analysis classes in the the new analysis framework
  *
 #include <TParticle.h>
 #include <TH2.h>
 #include <TList.h>
-#include "AliAnaGammaDirect.h" 
 #include "Riostream.h"
+#include "TROOT.h"
+
+// --- AliRoot system --- 
+#include "AliAnaGammaDirect.h" 
 #include "AliLog.h"
-  
+
 ClassImp(AliAnaGammaDirect)
   
 //____________________________________________________________________________
@@ -48,10 +54,12 @@ ClassImp(AliAnaGammaDirect)
     TObject(),
     fMinGammaPt(0.),
     fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
-    fICMethod(0),fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0),  
+    fICMethod(0), fAnaMC(0),
+    fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0), 
+    fntuplePrompt(0),
     //kSeveralIC
     fNCones(0),fNPtThres(0), fConeSizes(),  fPtThresholds(), 
-    fhPtThresIsolated(), fhPtSumIsolated()
+    fhPtThresIsolated(), fhPtSumIsolated(), fntSeveralIC()
 {
   //default ctor
   
@@ -68,7 +76,10 @@ AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
   fPtThreshold(g.fPtThreshold),
   fPtSumThreshold(g.fPtSumThreshold), 
   fICMethod(g.fICMethod),
-  fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),fhEtaGamma(g.fhEtaGamma),  
+  fAnaMC( g.fAnaMC),
+  fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),
+  fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),    
+  fntuplePrompt(g.fntuplePrompt),
   //kSeveralIC
   fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(), 
   fhPtThresIsolated(), fhPtSumIsolated()
@@ -79,12 +90,15 @@ AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
   for(Int_t i = 0; i < fNCones ; i++){
     fConeSizes[i] =  g.fConeSizes[i];
     fhPtSumIsolated[i] = g.fhPtSumIsolated[i]; 
+    fntSeveralIC[i]= g.fntSeveralIC[i];  
     for(Int_t j = 0; j < fNPtThres ; j++)
       fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j]; 
   }
   
   for(Int_t i = 0; i < fNPtThres ; i++)
     fPtThresholds[i]=   g.fPtThresholds[i];
+
+
 }
 
 //_________________________________________________________________________
@@ -99,10 +113,15 @@ AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & sou
   fPtThreshold = source.fPtThreshold ;
   fPtSumThreshold = source.fPtSumThreshold ; 
   fICMethod = source.fICMethod ;
+  fAnaMC = source.fAnaMC ;
+
   fhNGamma = source.fhNGamma ; 
   fhPhiGamma = source.fhPhiGamma ;
   fhEtaGamma = source.fhEtaGamma ;
-  
+  fhConeSumPt = source.fhConeSumPt ;
+
+  fntuplePrompt = source.fntuplePrompt ;
+
   //kSeveralIC
   fNCones = source.fNCones ;
   fNPtThres = source.fNPtThres ; 
@@ -110,13 +129,14 @@ AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & sou
   for(Int_t i = 0; i < fNCones ; i++){
     fConeSizes[i] =  source.fConeSizes[i];
     fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
+    fntSeveralIC[i]= source.fntSeveralIC[i];  
     for(Int_t j = 0; j < fNPtThres ; j++)
       fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
   }
   
   for(Int_t i = 0; i < fNPtThres ; i++)
     fPtThresholds[i]=   source.fPtThresholds[i];
-  
+
   return *this;
   
 }
@@ -129,11 +149,14 @@ AliAnaGammaDirect::~AliAnaGammaDirect()
   delete fhNGamma    ;  
   delete fhPhiGamma  ; 
   delete fhEtaGamma   ;  
-  
+  delete fhConeSumPt ;
+  delete fntuplePrompt    ;  
+
   //kSeveralIC
   delete [] fhPtThresIsolated ;
   delete [] fhPtSumIsolated ;
-  
+  delete [] fntSeveralIC ;
+
 }
 
 //________________________________________________________________________
@@ -163,6 +186,16 @@ TList *  AliAnaGammaDirect::GetCreateOutputObjects()
   fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
   outputContainer->Add(fhEtaGamma) ;
 
+  fhConeSumPt  = new TH2F
+    ("ConePtSum","#Sigma p_{T}  in cone ",200,0,120,100,0,100);
+  fhConeSumPt->SetYTitle("#Sigma p_{T}");
+  fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
+  outputContainer->Add(fhConeSumPt) ;
+  
+  //NTUPLE
+  fntuplePrompt = new TNtuple("ntuplePromptGamma", "Tree of prompt #gamma", "ptcluster:phicluster:etacluster:ptprimary:phiprimary:etaprimary:pdgprimary:statusprimary");
+  outputContainer->Add(fntuplePrompt) ;
+
   if(fICMethod == kSeveralIC){
     char name[128];
     char title[128];
@@ -174,6 +207,11 @@ TList *  AliAnaGammaDirect::GetCreateOutputObjects()
       fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
       outputContainer->Add(fhPtSumIsolated[icone]) ; 
       
+      sprintf(name,"nt_Cone_%d",icone);
+      sprintf(title,"ntuple for cone size %d",icone);
+      fntSeveralIC[icone] = new TNtuple(name, title, "ptcand:phicand:etacand:ptsum:type:ncone0:ncone1:ncone2:ncone3:ncone4:ncone5");
+      outputContainer->Add(fntSeveralIC[icone]) ; 
+
       for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
        sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
        sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
@@ -184,21 +222,26 @@ TList *  AliAnaGammaDirect::GetCreateOutputObjects()
     }//ipt loop
   }
 
+  gROOT->cd();
+
   return outputContainer ;
 
 }
 
 //____________________________________________________________________________
-void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &found) const 
+void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plCTS, TClonesArray * plPrimCalo, TParticle *pGamma, Bool_t &found) const 
 {
   //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
 
   Double_t pt = 0;
+  Int_t n = 0;
   Int_t index = -1; 
-  for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
-    TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
+  Float_t coneptsum = 0 ;
+
+  for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
+    TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
 
-    if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
+    if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt) && (particle->GetPdgCode() == 22)){
       index = ipr ;
       pt = particle->Pt();
       pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
@@ -209,10 +252,9 @@ void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS,
   //Do Isolation?
   if( ( fICMethod == kPtIC  ||  fICMethod == kSumPtIC )  && found)
     {
-      Float_t coneptsum = 0 ;
       Bool_t  icPtThres   = kFALSE;
       Bool_t  icPtSum     = kFALSE;
-      MakeIsolationCut(plCTS,pl, pGamma, index, 
+      MakeIsolationCut(plCTS,plCalo, pGamma, index,n,
                       icPtThres, icPtSum,coneptsum);
       if(fICMethod == kPtIC) //Pt thres method
        found = icPtThres ;
@@ -221,12 +263,38 @@ void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS,
     }
   
   if(found){
-    AliDebug(1,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
+    AliDebug(1,Form("Cluster with p_{T} larger than the pt treshold %f found in calorimeter ", fMinGammaPt)) ;
     AliDebug(1,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
     //Fill prompt gamma histograms 
-    fhNGamma->Fill(pGamma->Pt());
-    fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
-    fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
+    Float_t ptcluster = pGamma->Pt();
+    Float_t phicluster = pGamma->Phi();
+    Float_t etacluster = pGamma->Eta();
+
+    fhNGamma   ->Fill(ptcluster);
+    fhPhiGamma ->Fill(ptcluster,phicluster);
+    fhEtaGamma ->Fill(ptcluster,etacluster);
+    fhConeSumPt->Fill(ptcluster,coneptsum);
+
+    Float_t ptprimary = 0 ;
+    Float_t phiprimary = 0 ;
+    Float_t etaprimary = 0 ;
+    Int_t pdgprimary = 0 ;
+    Int_t statusprimary = 0 ;
+
+    if(fAnaMC){
+      TParticle * primary = dynamic_cast<TParticle *>(plPrimCalo->At(index)) ;
+      ptprimary = primary->Pt();
+      phiprimary = primary->Phi();
+      etaprimary = primary->Eta();
+      pdgprimary =  TMath::Abs(primary->GetPdgCode()) ;
+      statusprimary = primary->GetStatusCode(); // = 2 means decay gamma!!!
+      AliDebug(1, Form("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f",ptcluster,pdgprimary,ptprimary));
+      //printf("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f \n",ptcluster,pdgprimary,ptprimary);
+    }
+
+    //Fill ntuple with cluster / MC data
+    fntuplePrompt->Fill(ptcluster,phicluster,etacluster,ptprimary,phiprimary, etaprimary,pdgprimary,statusprimary);
   }
   else
     AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
@@ -240,27 +308,27 @@ void AliAnaGammaDirect::InitParameters()
   fMinGammaPt  = 5. ;
 
   //Fill particle lists when PID is ok
-  fConeSize             = 0.2 ; 
-  fPtThreshold         = 2.0; 
+  fConeSize             = 0.4 ; 
+  fPtThreshold         = 1.; 
   fPtSumThreshold  = 1.; 
 
   fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
-
+  fAnaMC = kFALSE ;
  //-----------kSeveralIC-----------------
-  fNCones           = 4 ; 
-  fNPtThres         = 4 ; 
-  fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
-  fPtThresholds[0]=1.; fPtThresholds[1]=2.; fPtThresholds[2]=3.; fPtThresholds[3]=4.;
+  fNCones           = 5 ; 
+  fNPtThres         = 5 ; 
+  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.;
 
 }
 
 //__________________________________________________________________
 void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS, 
-                                          TClonesArray * plNe, 
-                                          TParticle * pCandidate, 
-                                          Int_t index, 
-                                          Bool_t  &icmpt,  Bool_t  &icms, 
-                                          Float_t &coneptsum) const 
+                                         TClonesArray * plNe, 
+                                         TParticle * pCandidate, 
+                                         Int_t index, Int_t & n,
+                                         Bool_t  &icmpt,  Bool_t  &icms, 
+                                         Float_t &coneptsum) const 
 {  
   //Search in cone around a candidate particle if it is isolated 
   Float_t phiC  = pCandidate->Phi() ;
@@ -269,10 +337,10 @@ void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
   Float_t eta   = -100.  ;
   Float_t phi    = -100.  ;
   Float_t rad   = -100 ;
-  Int_t    n        = 0 ;
   TParticle * particle  = new TParticle;
 
-  coneptsum = 0; 
+  n = 0 ;
+  coneptsum = 0.; 
   icmpt = kFALSE;
   icms  = kFALSE;
 
@@ -284,9 +352,8 @@ void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
     phi  = particle->Phi() ;
     
     //Check if there is any particle inside cone with pt larger than  fPtThreshold
-    rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
-                     TMath::Power((phi-phiC),2));
-    if(rad<fConeSize){
+    rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
+    if(rad < fConeSize){
       AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
       coneptsum+=pt;
       if(pt > fPtThreshold ) n++;
@@ -302,9 +369,8 @@ void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
       phi  = particle->Phi() ;
       
       //Check if there is any particle inside cone with pt larger than  fPtThreshold
-      rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
-                       TMath::Power((phi-phiC),2));
-      if(rad<fConeSize){
+      rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
+      if(rad < fConeSize){
        AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
        coneptsum+=pt;
        if(pt > fPtThreshold ) n++;
@@ -314,7 +380,7 @@ void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
   
   if(n == 0) 
     icmpt =  kTRUE ;
-  if(coneptsum < fPtSumThreshold)
+  if(coneptsum <= fPtSumThreshold)
     icms  =  kTRUE ;
 
 }
@@ -323,41 +389,63 @@ void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
 void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS) 
 {
   //Isolation Cut Analysis for both methods and different pt cuts and cones
-
   if (fICMethod != kSeveralIC)
     AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
+  Int_t type = 0;
+
+  //Search maximum energy photon in the event
+  Double_t ptC = 0;
+  Int_t index = -1; 
+  TParticle * pCandidate = new TParticle();
+  Bool_t found = kFALSE;
   
-  for(Int_t ipr = 0; ipr < plCalo->GetEntries() ; ipr ++ ){
-    TParticle * pCandidate = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
+  for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
+    TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
     
-    if(pCandidate->Pt() > fMinGammaPt){
-      
-      Bool_t  icPtThres   = kFALSE;
-      Bool_t  icPtSum     = kFALSE;
-      
-      Float_t ptC      = pCandidate->Pt() ;
-   
-      fhNGamma->Fill(ptC);
-      fhPhiGamma->Fill( ptC,pCandidate->Phi());
-      fhEtaGamma->Fill(ptC,pCandidate->Eta());
+    if((particle->Pt() > fMinGammaPt) && (particle->Pt() > ptC) && (particle->GetPdgCode() == 22)){
+      index = ipr ;
+      ptC = particle->Pt();
+      pCandidate = particle ;
+      found  = kTRUE;
+    }
+  }
+  
+  if(found){
     
-      for(Int_t icone = 0; icone<fNCones; icone++){
-       fConeSize=fConeSizes[icone] ;
-       Float_t coneptsum = 0 ;
-       for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
-         fPtThreshold=fPtThresholds[ipt] ;
-         MakeIsolationCut(plCTS,plCalo, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
-         AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
-                         pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
-
-         fhPtThresIsolated[icone][ipt]->Fill(ptC); 
-       }//pt thresh loop
-       fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
-      }//cone size loop
-    }//min pt candidate
-  }//candidate loop
+    fhNGamma->Fill(ptC);
+    fhPhiGamma->Fill(ptC,pCandidate->Phi());
+    fhEtaGamma->Fill(ptC,pCandidate->Eta());
+    
+    Int_t ncone[fNCones][fNPtThres];
+    Bool_t  icPtThres   = kFALSE;
+    Bool_t  icPtSum     = kFALSE;
+    
+    for(Int_t icone = 0; icone<fNCones; icone++){
+      fConeSize=fConeSizes[icone] ;
+      Float_t coneptsum = 0 ;
+      for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
+       ncone[icone][ipt]=0;
+       fPtThreshold=fPtThresholds[ipt] ;
+       MakeIsolationCut(plCTS,plCalo, pCandidate, index,  
+                        ncone[icone][ipt]=0, icPtThres, icPtSum,coneptsum);
+       AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
+                       pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
+       if(ptC >15 && ptC < 25 && (icPtThres || icPtSum) && ipt ==0){
+         printf("R %0.1f, ptthres %1.1f, ptsum %1.1f, Candidate pt %2.2f,  eta %2.2f, phi %2.2f, pt in cone %2.2f, Isolated? ICPt %d, ICSum %d\n",
+                fConeSize,  fPtThreshold, fPtSumThreshold, ptC, pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), coneptsum, icPtThres, icPtSum);
+         //cout<<"mother label "<<pCandidate->GetMother(0)<<endl;
+       }
+       
+       fhPtThresIsolated[icone][ipt]->Fill(ptC); 
+      }//pt thresh loop
+      fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
+      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
+
 }
 
+//__________________________________________________________________
 void AliAnaGammaDirect::Print(const Option_t * opt) const
 {
   
index f99bdf8..4985f5d 100644 (file)
@@ -7,6 +7,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.5  2007/08/17 12:40:04  schutz
+ * New analysis classes by Gustavo Conesa
+ *
  * Revision 1.4.4.3  2007/07/26 10:32:09  schutz
  * new analysis classes in the the new analysis framework
  *
@@ -27,6 +30,7 @@
 #include <TClonesArray.h> 
 #include "TObject.h" 
 #include <TH2F.h>
+#include <TNtuple.h>
 
 class TList ;
 
@@ -46,13 +50,14 @@ public:
   Float_t     GetPtThreshold()      const {return fPtThreshold ; }
   Float_t     GetPtSumThres()     const {return fPtSumThreshold ; }
   Int_t        GetICMethod()          const {return fICMethod ; }
+  Bool_t     IsMC() const {return fAnaMC ; };
 
   TList *  GetCreateOutputObjects();
-  void GetPromptGamma(TClonesArray * plNe, TClonesArray * plCTS, TParticle * pGamma, Bool_t &Is)  const;
+  void GetPromptGamma(TClonesArray * plNe, TClonesArray * plCTS, TClonesArray * plPrimNe,  TParticle * pGamma, Bool_t &Is)  const;
   
   void MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS); 
   void MakeIsolationCut(TClonesArray * plCTS, TClonesArray * plNe, 
-                       TParticle *pCandidate, Int_t index, 
+                       TParticle *pCandidate, Int_t index, Int_t &n,
                        Bool_t &imcpt, Bool_t &icms, Float_t &ptsum) const ;  
   
   void Print(const Option_t * opt)const;
@@ -62,7 +67,8 @@ public:
   void SetPtThreshold(Float_t pt)        {fPtThreshold = pt; };
   void SetPtSumThreshold(Float_t pt) {fPtSumThreshold = pt; };
   void SetICMethod(Int_t i )          {fICMethod = i ; }
-  
+  void SetMC()    {fAnaMC = kTRUE ; }
+
   Int_t    GetNCones()                  const {return fNCones ; }
   Int_t    GetNPtThresholds()                const {return fNPtThres ; }
   Float_t GetConeSizes(Int_t i)      const {return fConeSizes[i] ; }
@@ -87,11 +93,16 @@ public:
                                            // kPtIC: Pt threshold method
                                            // kSumPtIC: Cone pt sum method
                                            // kSeveralIC: Analysis for several cuts
+  Bool_t fAnaMC ; //Set in case of using MCData reader 
+
   //Histograms  
   TH1F * fhNGamma    ;  //Number of (isolated) gamma identified
-  TH2F * fhPhiGamma    ; // Phi of identified gamma
-  TH2F * fhEtaGamma    ; // eta of identified gamma
-  
+  TH2F * fhPhiGamma  ; // Phi of identified gamma
+  TH2F * fhEtaGamma  ; // eta of identified gamma
+  TH2F * fhConeSumPt ; // Sum Pt in the cone
+
+  TNtuple *    fntuplePrompt ; //List of found prompt photons, pt, eta and phi. Also primary information.
+
   //Prompt photon analysis data members for multiple cones and pt thresholds kIsolationCut
   Int_t         fNCones   ; //Number of cone sizes to test
   Int_t         fNPtThres ; //Number of ptThres to test
@@ -100,8 +111,9 @@ public:
   
   TH1F* fhPtThresIsolated[20][20]; // Isolated gamma with pt threshold 
   TH2F* fhPtSumIsolated[20] ;  //  Isolated gamma with threshold on cone pt sume
+  TNtuple *    fntSeveralIC[20] ; //ntuple 
 
-  ClassDef(AliAnaGammaDirect,0)
+  ClassDef(AliAnaGammaDirect,1)
 } ;
  
 
index fc9f2c4..c8d5fbc 100644 (file)
@@ -17,6 +17,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.5  2007/08/17 12:40:04  schutz
+ * New analysis classes by Gustavo Conesa
+ *
  * Revision 1.4.4.2  2007/07/26 10:32:09  schutz
  * new analysis classes in the the new analysis framework
  *
@@ -45,7 +48,8 @@ ClassImp(AliAnaGammaHadron)
     AliAnaGammaCorrelation(),
     fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0), 
     fhDeltaPhiGammaCharged(0),  fhDeltaPhiGammaNeutral(0), 
-    fhDeltaEtaGammaCharged(0), fhDeltaEtaGammaNeutral(0), 
+    fhDeltaEtaGammaCharged(0), fhDeltaEtaGammaNeutral(0),
+    fhDeltaPhiChargedPt(0),  
     fhCorrelationGammaNeutral(0), fhCorrelationGammaCharged(0)
 {
   //Default Ctor
@@ -64,6 +68,7 @@ AliAnaGammaHadron::AliAnaGammaHadron(const AliAnaGammaHadron & g) :
   fhDeltaPhiGammaNeutral(g.fhDeltaPhiGammaNeutral), 
   fhDeltaEtaGammaCharged(g.fhDeltaEtaGammaCharged), 
   fhDeltaEtaGammaNeutral(g.fhDeltaEtaGammaNeutral), 
+  fhDeltaPhiChargedPt(g.fhDeltaPhiChargedPt), 
   fhCorrelationGammaNeutral(g.fhCorrelationGammaNeutral), 
   fhCorrelationGammaCharged(g.fhCorrelationGammaCharged)
 {
@@ -85,6 +90,7 @@ AliAnaGammaHadron & AliAnaGammaHadron::operator = (const AliAnaGammaHadron & sou
   fhDeltaPhiGammaNeutral = source.fhDeltaPhiGammaNeutral ; 
   fhDeltaEtaGammaCharged = source.fhDeltaEtaGammaCharged ; 
   fhDeltaEtaGammaNeutral = source.fhDeltaEtaGammaNeutral ; 
+  fhDeltaPhiChargedPt = source.fhDeltaPhiChargedPt ;
 
   fhCorrelationGammaNeutral = source.fhCorrelationGammaNeutral ; 
   fhCorrelationGammaCharged = source.fhCorrelationGammaCharged ; 
@@ -105,6 +111,7 @@ AliAnaGammaHadron::~AliAnaGammaHadron()
   delete fhDeltaPhiGammaNeutral   ; 
   delete fhDeltaEtaGammaCharged  ; 
   delete fhDeltaEtaGammaNeutral  ; 
+  delete fhDeltaPhiChargedPt  ;
 
   delete fhCorrelationGammaNeutral  ; 
   delete fhCorrelationGammaCharged  ;
@@ -138,7 +145,13 @@ TList *  AliAnaGammaHadron::GetCreateOutputObjects()
      200,0,120,200,0,6.4); 
   fhDeltaPhiGammaCharged->SetYTitle("#Delta #phi");
   fhDeltaPhiGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
-  
+
+  fhDeltaPhiChargedPt  = new TH2F
+    ("DeltaPhiChargedPt","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #pi}",
+     200,0,120,200,0,6.4);
+  fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
+  fhDeltaPhiChargedPt->SetXTitle("p_{T #pi} (GeV/c)");
+
   fhDeltaEtaGammaCharged  = new TH2F
     ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
      200,0,120,200,-2,2); 
@@ -156,6 +169,7 @@ TList *  AliAnaGammaHadron::GetCreateOutputObjects()
   outputContainer->Add(fhDeltaPhiGammaCharged) ; 
   outputContainer->Add(fhDeltaEtaGammaCharged) ;
   outputContainer->Add(fhCorrelationGammaCharged) ;
+  outputContainer->Add(fhDeltaPhiChargedPt) ;
 
   if(!AreJetsOnlyInCTS()){
     //---- kHadron and kJetLeadCone ----
@@ -269,6 +283,8 @@ void  AliAnaGammaHadron::MakeGammaChargedCorrelation(TParticle * pGamma, TClones
     fhPhiCharged->Fill(ptg,phi);
     fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-particle->Eta());
     fhDeltaPhiGammaCharged->Fill(ptg,phig-phi);
+    fhDeltaPhiChargedPt->Fill(pt,phig-phi);
+
     //Selection within angular and energy limits
     if(((phig-phi)> GetDeltaPhiMinCut()) && ((phig-phi)<GetDeltaPhiMaxCut()) && pt > GetMinPtHadron()){
       AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
index ce5358e..5f8d567 100644 (file)
@@ -7,6 +7,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.4  2007/08/17 12:40:04  schutz
+ * New analysis classes by Gustavo Conesa
+ *
  * Revision 1.3.4.2  2007/07/26 10:32:09  schutz
  * new analysis classes in the the new analysis framework
  *
@@ -49,11 +52,12 @@ public:
   TH2F * fhDeltaPhiGammaNeutral   ; 
   TH2F * fhDeltaEtaGammaCharged  ; 
   TH2F * fhDeltaEtaGammaNeutral  ; 
+  TH2F * fhDeltaPhiChargedPt  ;
 
   TH2F * fhCorrelationGammaNeutral  ; //Neutral hadron correlation histogram 
   TH2F * fhCorrelationGammaCharged  ; //Charged hadron correlation histogram
   
-  ClassDef(AliAnaGammaHadron,0)
+  ClassDef(AliAnaGammaHadron,1)
 } ;
  
 
index fcf6657..435eafb 100644 (file)
@@ -7,6 +7,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
  *
@@ -43,7 +46,7 @@ class AliAnaGammaJetFinder : public AliAnaGammaCorrelation {
        TH2F * fhDeltaPtJet;
        TH2F * fhPtRatJet;
        
-       ClassDef(AliAnaGammaJetFinder,0)
+       ClassDef(AliAnaGammaJetFinder,1)
  } ;
 
 
index 7b77efa..4f0f268 100644 (file)
@@ -7,6 +7,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
  *
@@ -182,7 +185,7 @@ public:
   TH2F * fhBkgPtDists[5][5];
   
 
-  ClassDef(AliAnaGammaJetLeadCone,0)
+  ClassDef(AliAnaGammaJetLeadCone,1)
 } ;
  
 
index 37185fa..8642ea9 100644 (file)
@@ -7,6 +7,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.3  2007/09/26 11:07:19  schutz
+ * Update classes for the new analysis framwork
+ *
  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  * new analysis classes in the the new analysis framework
  *
@@ -43,7 +46,7 @@ class AliAnaGammaParton : public AliAnaGammaCorrelation {
        TH2F * fhDeltaPtParton;
        TH2F * fhPtRatParton;
        
-       ClassDef(AliAnaGammaParton,0)
+       ClassDef(AliAnaGammaParton,1)
  } ;
 
 
index 72db8ee..9a09834 100644 (file)
@@ -127,7 +127,7 @@ void AliAnaGammaPhos::CreateOutputObjects()
   AliAODHandler* handler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());  
   fTreeA = handler->GetTree() ; 
   fAOD   = handler->GetAOD();
-  fAODPhotons = fAOD->GetClusters() ; 
+  fAODPhotons = fAOD->GetCaloClusters() ; 
   
 
   OpenFile(1) ; 
@@ -192,7 +192,7 @@ void AliAnaGammaPhos::Exec(Option_t *)
       fhPHOSPos->Fill( pos[0], pos[1], pos[2] ) ;
       fhPHOSDigits->Fill(entry, caloCluster->GetNumberOfDigits() ) ;
       numberOfDigitsInPhos += caloCluster->GetNumberOfDigits() ;
-      Float_t * pid = caloCluster->GetPid() ;
+      Double_t * pid = caloCluster->GetPid() ;
       if(pid[AliPID::kPhoton] > GetPhotonId() ) {
        phosVector[fPhotonsInPhos] = new TVector3(pos[0],pos[1],pos[2]) ;
        phosPhotonsEnergy[fPhotonsInPhos]=caloCluster->E() ;
index 8942716..2963e8b 100644 (file)
@@ -59,6 +59,6 @@ private:
   TH1D    * fhPHOSInvariantMass ;
   TH1I    * fhPHOSDigitsEvent ;
    
-  ClassDef(AliAnaGammaPhos, 0); // a PHOS photon analysis task 
+  ClassDef(AliAnaGammaPhos, 1); // a PHOS photon analysis task 
 };
 #endif // ALIANAGAMMAPHOS_H
index 0769adc..f4a2f60 100644 (file)
 #include "AliAnaScale.h" 
 #include "AliAnalysisManager.h"
 #include "AliLog.h"
+#include "Riostream.h"
 
 //______________________________________________________________________________
 AliAnaScale::AliAnaScale() : 
   fDebug(0),
   fScale(1.0),
   fInputList(0x0), 
-  fOutputList(0x0), 
-  fhInPHOSEnergy(0),
-  fhOuPHOSEnergy(0)
+  fOutputList(0x0) 
 {
   //Default constructor
 }
@@ -48,9 +47,7 @@ AliAnaScale::AliAnaScale(const char *name) :
   fDebug(0),
   fScale(1.0), 
   fInputList(0x0), 
-  fOutputList(0x0), 
-  fhInPHOSEnergy(0),
-  fhOuPHOSEnergy(0)
+  fOutputList(0x0) 
 {
   // Constructor.
   // Called only after the event loop
@@ -76,7 +73,6 @@ void AliAnaScale::ConnectInputData(const Option_t*)
     
   AliInfo(Form("*** Initialization of %s", GetName())) ; 
   fInputList     = dynamic_cast<TList*>(GetInputData(0)) ;  
-  fhInPHOSEnergy = dynamic_cast<TH1D*>(fInputList->At(2));
 }
 //________________________________________________________________________
 void AliAnaScale::CreateOutputObjects()
@@ -91,15 +87,23 @@ void AliAnaScale::Exec(Option_t *)
 {
   // Do the Scaling
     
-  fhOuPHOSEnergy = static_cast<TH1D*>(fhInPHOSEnergy->Clone("PHOSEnergyScaled")) ;
-  
-  // create output container
-  
   fOutputList = new TList() ; 
   fOutputList->SetName(GetName()) ; 
-
-  fOutputList->AddAt(fhOuPHOSEnergy,          0) ; 
-  fhOuPHOSEnergy->Scale(fScale) ; 
+  TIter next(fInputList) ;     
+  TObject * h ; 
+  while ( (h = next()) ) { 
+    if(h){
+      if ( strcmp(h->ClassName(),"TNtuple") ) {
+      char name[20] ; 
+      sprintf(name, "%sScaled", h->GetName()) ; 
+      TH1 * hout = dynamic_cast<TH1*> (h->Clone(name)) ; 
+      hout->Scale(fScale) ;  
+      fOutputList->Add(hout) ; 
+      } 
+      else  fOutputList->Add(h) ; 
+    }
+  }
+  cout<<"end"<<endl;
   PostData(0, fOutputList);
 }
 
@@ -117,24 +121,5 @@ void AliAnaScale::Terminate(Option_t *)
 {
   // Processing when the event loop is ended
   
-  AliInfo(Form(" *** %s Report:", GetName())) ; 
-  printf("        PHOS Energy Integral In         : %5.3e \n", fhInPHOSEnergy->Integral() ) ;
-  printf("        PHOS Energy Integral Ou         : %5.3e \n", fhOuPHOSEnergy->Integral() ) ;
-
-  TCanvas  * cPHOS = new TCanvas("cPHOS", "PHOS ESD Test", 400, 10, 600, 700) ;
-
-  cPHOS->Divide(2, 1);
-  cPHOS->cd(1) ; 
-  if ( fhInPHOSEnergy->GetMaximum() > 0. ) 
-    gPad->SetLogy();
-  fhInPHOSEnergy->SetAxisRange(0, 25.);
-  fhInPHOSEnergy->SetLineColor(2);
-  fhInPHOSEnergy->Draw();
 
-  cPHOS->cd(2) ; 
-  if ( fhOuPHOSEnergy->GetMaximum() > 0. ) 
-    gPad->SetLogy();
-  fhOuPHOSEnergy->SetAxisRange(0, 25.);
-  fhOuPHOSEnergy->SetLineColor(2);
-  fhOuPHOSEnergy->Draw();
 }
index ac9b2e2..b09e1a7 100644 (file)
@@ -38,9 +38,7 @@ private:
   // Histograms
   TList   * fInputList ;  //! input data list
   TList   * fOutputList ; //! output data list
-  TH1D    * fhInPHOSEnergy ;
-  TH1D    * fhOuPHOSEnergy ;
    
-  ClassDef(AliAnaScale, 0); // a post event loop scaling 
+  ClassDef(AliAnaScale, 1); // a post event loop scaling 
 };
 #endif // ALIANASCALE_H
index ba14f06..d680967 100644 (file)
@@ -12,7 +12,8 @@
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
+
+// root
 #include <TROOT.h>
 #include <TSystem.h>
 #include <TInterpreter.h>
 #include <TFile.h>
 #include <Riostream.h>
 
+// analysis
 #include "AliAnalysisTaskGamma.h"
 #include "AliAnalysisManager.h"
 #include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
 #include "AliAnaGamma.h"
 #include "AliGammaReader.h"
 #include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliAODHandler.h"
 #include "AliStack.h"
 #include "AliLog.h"
 
@@ -37,7 +42,8 @@ AliAnalysisTaskGamma::AliAnalysisTaskGamma():
     fAna(0x0),
     fChain(0x0),
     fESD(0x0),
-    fTreeG(0x0),//Not used for the moment
+    fAOD(0x0),
+    fTreeG(0x0),
     fOutputContainer(0x0)
 {
   // Default constructor
@@ -49,41 +55,48 @@ AliAnalysisTaskGamma::AliAnalysisTaskGamma(const char* name):
     fAna(0x0),
     fChain(0x0),
     fESD(0x0),
-    fTreeG(0x0),//Not used for the moment
+    fAOD(0x0),
+    fTreeG(0x0),
     fOutputContainer(0x0)
 {
   // Default constructor
-  
-  Init();
  
   DefineInput (0, TChain::Class());
-  //DefineOutput(0, TTree::Class());// to create AODs, to be done
-  DefineOutput(0, TList::Class());
+  DefineOutput(0, TTree::Class());
+  DefineOutput(1, TList::Class());
 
 }
 
 //_____________________________________________________
 AliAnalysisTaskGamma::~AliAnalysisTaskGamma() 
 {
-
   // Remove all pointers
-  fOutputContainer->Clear() ; 
-  delete fOutputContainer ;
+  if(fOutputContainer){
+    fOutputContainer->Clear() ; 
+    delete fOutputContainer ;
+  }
+  
+  if(fTreeG) delete fTreeG ; 
 
-  delete fTreeG ; //Not used for the moment
 }
 
 //_____________________________________________________
 void AliAnalysisTaskGamma::CreateOutputObjects()
 {
   // Create the output container
-  //OpenFile(0);
-  //fTreeG = new TTree ; // fAna->MakeTreeG("TreeG");// to create AODs, to be done
   
+  //AODs
   OpenFile(0);
-  fOutputContainer = fAna->GetOutputContainer();
+  AliAODHandler* handler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+  fAOD   = handler->GetAOD();
+  fTreeG = handler->GetTree();
+  fAna->ConnectAOD(fAOD);
 
+  //Histograms container
+  OpenFile(1);
+  fOutputContainer = fAna->GetOutputContainer();
+  
 }
 
 //_____________________________________________________
@@ -91,7 +104,7 @@ void AliAnalysisTaskGamma::Init()
 {
   // Initialization
   AliDebug(1,"Begin");
-   
+  
   // Call configuration file
   gROOT->LoadMacro("ConfigGammaAnalysis.C");
   fAna = (AliAnaGamma*) gInterpreter->ProcessLine("ConfigGammaAnalysis()");
@@ -146,10 +159,9 @@ void AliAnalysisTaskGamma::Exec(Option_t */*option*/)
     AliMCEventHandler*    mctruth = (AliMCEventHandler*) 
       ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
     
-    if(mctruth){
-      stack = mctruth->Stack();
-      //printf("AliAnalysisTaskJets: Number of tracks %5d\n", stack->GetNtrack());
-    }
+    if(mctruth)
+      stack = mctruth->MCEvent()->Stack();
+    
   }
   
   //Get Event
@@ -162,6 +174,7 @@ void AliAnalysisTaskGamma::Exec(Option_t */*option*/)
     AliError("fESD is not connected to the input!") ; 
     return ; 
   } 
+  
   fAna->SetData(fESD);
   
   //In case of montecarlo analysis, pass the stack also.
@@ -171,8 +184,8 @@ void AliAnalysisTaskGamma::Exec(Option_t */*option*/)
   //Process event
   fAna->ProcessEvent(ientry);
   
-  //PostData(0, fTreeG); // Create AODs, to be done.
-  PostData(0, fOutputContainer);
+  PostData(0, fTreeG); 
+  PostData(1, fOutputContainer);
   
 }
 
index bbe3042..7fe978b 100644 (file)
@@ -7,6 +7,7 @@
 #include "AliAnalysisTask.h"
 class AliAnaGamma;
 class AliESDEvent;
+class AliAODEvent;
 class TChain;
 class TList;
 
@@ -21,7 +22,7 @@ class AliAnalysisTaskGamma : public AliAnalysisTask
     virtual void ConnectInputData(Option_t *option = "");
     virtual void CreateOutputObjects();
     virtual void Init();
-    //virtual void LocalInit() {Init();}
+    virtual void LocalInit() {Init();}
     virtual void Exec(Option_t *option);
     virtual void Terminate(Option_t *option);
    
@@ -30,10 +31,11 @@ class AliAnalysisTaskGamma : public AliAnalysisTask
     AliAnaGamma* fAna; //  Pointer to the jet finder 
     TChain*       fChain;     //! chained files
     AliESDEvent*       fESD;       //! ESD
+    AliAODEvent*       fAOD;       //! AOD
     TTree*        fTreeG;     //  tree of prompt gamma, does nothing for the moment 
     TList * fOutputContainer ; // Histogram container
 
-    ClassDef(AliAnalysisTaskGamma, 0); // Analysis task for standard gamma correlation analysis
+    ClassDef(AliAnalysisTaskGamma, 1); // Analysis task for standard gamma correlation analysis
 };
  
 #endif //ALIANALYSISTASKGAMMA_H
index 42c0f07..3ba3e68 100644 (file)
@@ -18,6 +18,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
  *
@@ -48,27 +51,18 @@ ClassImp(AliGammaDataReader)
 
 //____________________________________________________________________________
 AliGammaDataReader::AliGammaDataReader() : 
-  AliGammaReader(),  
-  fEMCALPID(0),fPHOSPID(0),
-  fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.), fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.)  
+  AliGammaReader()
 {
   //Default Ctor
   
   //Initialize parameters
   fDataType=kData;
-  InitParameters();
   
 }
 
 //____________________________________________________________________________
 AliGammaDataReader::AliGammaDataReader(const AliGammaDataReader & g) :   
-  AliGammaReader(g),
-  fEMCALPID(g.fEMCALPID), 
-  fPHOSPID(g.fPHOSPID),
-  fEMCALPhotonWeight(g.fEMCALPhotonWeight), 
-  fEMCALPi0Weight(g.fEMCALPi0Weight), 
-  fPHOSPhotonWeight(g.fPHOSPhotonWeight),
-  fPHOSPi0Weight(g.fPHOSPi0Weight)
+  AliGammaReader(g)
 {
   // cpy ctor
 }
@@ -80,13 +74,6 @@ AliGammaDataReader & AliGammaDataReader::operator = (const AliGammaDataReader &
 
   if(&source == this) return *this;
 
-  fEMCALPID = source.fEMCALPID ;
-  fPHOSPID = source.fPHOSPID ;
-  fEMCALPhotonWeight = source. fEMCALPhotonWeight ;
-  fEMCALPi0Weight = source.fEMCALPi0Weight ;
-  fPHOSPhotonWeight = source.fPHOSPhotonWeight ;
-  fPHOSPi0Weight = source.fPHOSPi0Weight ;
-
   return *this;
 
 }
@@ -95,15 +82,19 @@ AliGammaDataReader & AliGammaDataReader::operator = (const AliGammaDataReader &
 void AliGammaDataReader::CreateParticleList(TObject * data, TObject *,
                                            TClonesArray * plCTS, 
                                            TClonesArray * plEMCAL,  
-                                           TClonesArray * plPHOS, TClonesArray*){
+                                           TClonesArray * plPHOS, 
+                                           TClonesArray * , 
+                                           TClonesArray * ,  
+                                           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.
 
   AliESDEvent* esd = (AliESDEvent*) data;
 
   Int_t npar  = 0 ;
-  Float_t *pid = new Float_t[AliPID::kSPECIESN];  
+  Double_t *pid = new Double_t[AliPID::kSPECIESN];  
    AliDebug(3,"Fill particle lists");
   
   //Get vertex for momentum calculation  
@@ -141,21 +132,52 @@ void AliGammaDataReader::CreateParticleList(TObject * data, TObject *,
        
        pid=clus->GetPid();     
        Int_t pdg = 22;
-       
-       if(fPHOSPID){
-         AliDebug(4, Form("PID: photon %f, pi0 %f, electron %f",pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron]));
-         if( pid[AliPID::kPhoton] > fPHOSPhotonWeight) pdg=22;
-         else if( pid[AliPID::kPi0] > fPHOSPi0Weight) pdg=111;
-         else pdg = 0;
+
+       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.
        }
        
-       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(3,Form("PHOS added: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
-       
-       new((*plPHOS)[indexPH++])   TParticle(*particle) ;
-       
+       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()));        
+
       }//pt, eta, phi cut
       else     AliDebug(4,"Particle not added");
     }//not charged
@@ -202,18 +224,13 @@ void AliGammaDataReader::CreateParticleList(TObject * data, TObject *,
   //################ EMCAL ##############
   
   Int_t indexEM  = plEMCAL->GetEntries() ; 
-  //Int_t begem = esd->GetFirstEMCALCluster();  
-  //Int_t endem = esd->GetFirstEMCALCluster() + 
-  //esd->GetNumberOfEMCALClusters() ;  
-  
-  //AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
   
   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::kClusterv1 && clus->GetTrackMatched()==-1 ){
+    if(clustertype == AliESDCaloCluster::kEMCALClusterv1 && clus->GetTrackMatched()==-1 ){
       
       TLorentzVector momentum ;
       clus->GetMomentum(momentum, v); 
@@ -224,57 +241,49 @@ void AliGammaDataReader::CreateParticleList(TObject * data, TObject *,
       
        pid=clus->GetPid();     
        Int_t pdg = 22;
-       if(fEMCALPID){
-         AliDebug(4, Form("PID: photon %f, pi0 %f, electron %f",pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron]));
-         if( pid[AliPID::kPhoton] > fEMCALPhotonWeight) pdg=22;
-         else if( pid[AliPID::kPi0] > fEMCALPi0Weight) pdg=111;
-         else pdg = 0 ;
-       }
 
+       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 ;
+       }
        
-       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(3,Form("EMCAL added: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
+       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((*plEMCAL)[indexEM++])   TParticle(*particle) ;
-  
       }//pt, phi, eta cut
       else     AliDebug(4,"Particle not added");
     }//not charged, not pseudocluster
   }//Cluster loop
   
-  AliDebug(3,"Particle lists filled");
-  
-}
-
-  //____________________________________________________________________________
-void AliGammaDataReader::InitParameters()
-{
-  //Initialize the parameters of the analysis.
-
-  //Fill particle lists when PID is ok
-  fEMCALPID = kFALSE;
-  fPHOSPID = kFALSE;
-  fEMCALPhotonWeight = 0.5 ;
-  fEMCALPi0Weight = 0.5 ;
-  fPHOSPhotonWeight = 0.8 ;
+  AliDebug(3,Form("Particle lists filled, PHOS clusters %d, EMCAL clusters %d, CTS tracks %d", 
+                 indexPH,indexEM,indexCh));
 
 }
 
-
-void AliGammaDataReader::Print(const Option_t * opt) const
-{
-
-  //Print some relevant parameters set for the analysis
-  if(! opt)
-    return;
-
-  Info("Print", "%s %s", GetName(), GetTitle() ) ;
-  printf("PHOS PID on?               =     %d\n",  fPHOSPID) ; 
-  printf("EMCAL PID  on?         =     %d\n",  fEMCALPID) ;
-  printf("PHOS PID weight , photon  %f\n",  fPHOSPhotonWeight) ; 
-  printf("EMCAL PID weight, photon %f, pi0 %f\n",   fEMCALPhotonWeight,  fEMCALPi0Weight) ; 
-
-}
index 1fa8838..f4a62f5 100644 (file)
@@ -7,6 +7,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.2  2007/08/21 08:38:20  schutz
+ * Missing heared file added
+ *
  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  * new analysis classes in the the new analysis framework
  *
@@ -36,39 +39,12 @@ public:
   AliGammaDataReader & operator = (const AliGammaDataReader & g) ;//cpy assignment
   virtual ~AliGammaDataReader() {;} //virtual dtor
 
-  void InitParameters();
-
-  Bool_t   IsEMCALPIDOn() const {return fEMCALPID ; }
-  Bool_t   IsPHOSPIDOn() const {return fPHOSPID ; }
-  Float_t  GetEMCALPhotonWeight() { return  fEMCALPhotonWeight  ; }
-  Float_t  GetEMCALPi0Weight()    {  return fEMCALPi0Weight  ; }
-  Float_t  GetPHOSPhotonWeight()  {  return fPHOSPhotonWeight  ; }
-  Float_t  GetPHOSPi0Weight()  {  return fPHOSPi0Weight  ; }
-  
-  void Print(const Option_t * opt)const;
-  
-  void SetEMCALPIDOn(Bool_t pid){ fEMCALPID= pid ; }
-  void SetPHOSPIDOn(Bool_t pid){ fPHOSPID= pid ; }
-  void SetEMCALPhotonWeight(Float_t  w){  fEMCALPhotonWeight = w ; }
-  void SetEMCALPi0Weight(Float_t  w){  fEMCALPi0Weight = w ; }
-  void SetPHOSPhotonWeight(Float_t  w){  fPHOSPhotonWeight = w ; }
-  void SetPHOSPi0Weight(Float_t  w){  fPHOSPi0Weight = w ; }
-
   void CreateParticleList(TObject * esd, TObject *,TClonesArray * plCh, 
-                         TClonesArray * plEMCAL, TClonesArray * plPHOS, TClonesArray *);
-
+                         TClonesArray * plEMCAL, TClonesArray * plPHOS, TClonesArray *,TClonesArray *, TClonesArray *);
   
  private:
 
-  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 
-  Float_t      fEMCALPi0Weight;  //Bayesian PID weight for pi0 in EMCAL 
-  Float_t      fPHOSPhotonWeight; //Bayesian PID weight for photons in PHOS 
-  Float_t      fPHOSPi0Weight; //Bayesian PID weight for pi0 in PHOS 
-
-  ClassDef(AliGammaDataReader,0)
+  ClassDef(AliGammaDataReader,1)
 } ;
  
 
index 8c11c70..474ebd0 100644 (file)
@@ -18,6 +18,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
  *
@@ -49,27 +52,18 @@ ClassImp(AliGammaMCDataReader)
 
 //____________________________________________________________________________
   AliGammaMCDataReader::AliGammaMCDataReader() : 
-    AliGammaReader(),  
-    fEMCALPID(0),fPHOSPID(0),
-    fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.), fPHOSPhotonWeight(0.),  fPHOSPi0Weight(0.) 
+    AliGammaReader()
 {
   //Default Ctor
   
   //Initialize parameters
   fDataType=kMCData;
-  InitParameters();
   
 }
 
 //____________________________________________________________________________
 AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) :   
-  AliGammaReader(g),
-  fEMCALPID(g.fEMCALPID), 
-  fPHOSPID(g.fPHOSPID),
-  fEMCALPhotonWeight(g.fEMCALPhotonWeight), 
-  fEMCALPi0Weight(g.fEMCALPi0Weight), 
-  fPHOSPhotonWeight(g.fPHOSPhotonWeight),
-  fPHOSPi0Weight(g.fPHOSPi0Weight)
+  AliGammaReader(g)
 {
   // cpy ctor
 }
@@ -81,13 +75,7 @@ AliGammaMCDataReader & AliGammaMCDataReader::operator = (const AliGammaMCDataRea
   
   if(&source == this) return *this;
   
-  fEMCALPID = source.fEMCALPID ;
-  fPHOSPID = source.fPHOSPID ;
-  fEMCALPhotonWeight = source. fEMCALPhotonWeight ;
-  fEMCALPi0Weight = source.fEMCALPi0Weight ;
-  fPHOSPhotonWeight = source.fPHOSPhotonWeight ;
-  fPHOSPi0Weight = source.fPHOSPi0Weight ;
-  
+
   return *this;
   
 }
@@ -96,7 +84,10 @@ AliGammaMCDataReader & AliGammaMCDataReader::operator = (const AliGammaMCDataRea
 void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
                                              TClonesArray * plCTS, 
                                              TClonesArray * plEMCAL,  
-                                             TClonesArray * plPHOS, TClonesArray *){
+                                             TClonesArray * plPHOS,   
+                                             TClonesArray * plPrimCTS, 
+                                             TClonesArray * plPrimEMCAL, 
+                                             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 
@@ -105,7 +96,7 @@ void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
   AliStack* stack = (AliStack*) kine;
   
   Int_t npar  = 0 ;
-  Float_t *pid = new Float_t[AliPID::kSPECIESN];  
+  Double_t *pid = new Double_t[AliPID::kSPECIESN];  
   AliDebug(3,"Fill particle lists");
   
   //Get vertex for momentum calculation  
@@ -114,43 +105,91 @@ void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
   
   //########### PHOS ##############
   
-  Int_t begphos = esd->GetFirstPHOSCluster();  
-  Int_t endphos = esd->GetFirstPHOSCluster() + 
-    esd->GetNumberOfPHOSClusters() ;  
+  Int_t begphos = 0;  
+  Int_t endphos = esd->GetNumberOfCaloClusters() ;  
   Int_t indexPH = plPHOS->GetEntries() ;
-  AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
+  Int_t begem = 0;
+  Int_t endem = endphos;
+  
+  AliDebug(3,Form("First Calorimeter particle %d, last particle %d", begphos,endphos));
   
   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);
-      pid=clus->GetPid();      
-      Int_t pdg = 22;
+      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(fPHOSPID){
-       if( pid[AliPID::kPhoton] > fPHOSPhotonWeight) pdg=22;
-       if( pid[AliPID::kPi0] > fPHOSPi0Weight) pdg=111;
-       else pdg = 0 ;
-      }
-      
-      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);
+       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(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array
 
-      AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
-      
-      //Look  if parent is prompt photon
-      Int_t label = clus->GetLabel();
-      if(label>=0){
-       TParticle * pmother = stack->Particle(label);
-       Int_t imother = pmother->GetFirstMother(); 
-       if(imother == 6 || imother == 7)
-         pmother->SetFirstMother(22);
-      }
-      
-      new((*plPHOS)[indexPH++])   TParticle(*particle) ;
-    
+         //###############
+         //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()));
+         //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++;
+         
+       }
+       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, eta, phi cut
     }//not charged
   }//cluster loop
 
@@ -182,91 +221,209 @@ void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
       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) ;
-      new((*plCTS)[indexCh++])       TParticle(*particle) ;    
+
+      //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() ; 
-  Int_t begem = esd->GetFirstEMCALCluster();  
-  Int_t endem = esd->GetFirstEMCALCluster() + 
-    esd->GetNumberOfEMCALClusters() ;  
   
-  AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
   
   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::kClusterv1 && clus->GetTrackMatched()==-1 ){
+    if(clustertype == AliESDCaloCluster::kEMCALClusterv1 && clus->GetTrackMatched()==-1 ){
       
       TLorentzVector momentum ;
-      clus->GetMomentum(momentum, v);            
+      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(fEMCALPID){
-       if( pid[AliPID::kPhoton] > fEMCALPhotonWeight) pdg=22;
-       else if( pid[AliPID::kPi0] > fEMCALPi0Weight) pdg=111;
-       else pdg = 0;
+      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 ;
       }
       
-      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 clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
-
-      //Look  if parent is prompt photon
-      Int_t label = clus->GetLabel();
-      if(label>=0){
-       TParticle * pmother = stack->Particle(label);
-       Int_t imother = pmother->GetFirstMother(); 
-       if(imother == 6 || imother == 7)
-         pmother->SetFirstMother(22);
+      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()));
       
-      new((*plEMCAL)[indexEM++])   TParticle(*particle) ;
+      }//pt, phi, eta cut
+      else     AliDebug(4,"Particle not added");
       
     }//not charged, not pseudocluster
   }//Cluster loop
   
-  AliDebug(3,"Particle lists filled");
+  AliDebug(3,Form("Particle lists filled, PHOS clusters %d, EMCAL clusters %d, CTS tracks %d", 
+                 indexPH,indexEM,indexCh));
   
 }
 
-  //____________________________________________________________________________
-void AliGammaMCDataReader::InitParameters()
+
+TParticle * AliGammaMCDataReader:: GetMotherParticle(Int_t label, AliStack *stack, TString calo,  TLorentzVector momentum)
 {
-  //Initialize the parameters of the analysis.
+  //Gets the primary particle and do some checks:
+  //Check if primary is inside calorimeter and look the mother outsie
+  //Check if mother is a decay photon, in which case check if decay was overlapped
+  
+  Float_t minangle = 0;
+  Float_t ipdist = 0;
+  TParticle * pmother = new TParticle();
 
-  //Fill particle lists when PID is ok
-  fEMCALPID = kFALSE;
-  fPHOSPID = kFALSE;
-  fEMCALPhotonWeight = 0.5 ;
-  fEMCALPi0Weight = 0.5 ;
-  fPHOSPhotonWeight = 0.8 ;
-  fPHOSPi0Weight = 0.5 ;
+  if(calo == "PHOS"){
+    ipdist= fPHOSIPDistance;
+    minangle = fPHOSMinAngle; 
+  }
+  else if (calo == "EMCAL"){
+    ipdist = fEMCALIPDistance;
+    minangle = fEMCALMinAngle;
+  }
 
-}
 
+  if(label>=0){
+    pmother = stack->Particle(label);
+    Int_t mostatus = pmother->GetStatusCode();
+    Int_t mopdg    = TMath::Abs(pmother->GetPdgCode());
+    
+    //Check if mother is a conversion inside the calorimeter
+    //In such case, return the mother before the calorimeter.
+    //First approximation.
+    Float_t vy = TMath::Abs(pmother->Vy()) ;
 
-void AliGammaMCDataReader::Print(const Option_t * opt) const
-{
+    if( mostatus == 0 && vy >= ipdist){
 
-  //Print some relevant parameters set for the analysis
-  if(! opt)
-    return;
-  
-  Info("Print", "%s %s", GetName(), GetTitle() ) ;
-  printf("PHOS PID on?               =     %d\n",  fPHOSPID) ; 
-  printf("EMCAL PID  on?         =     %d\n",  fEMCALPID) ;
-  printf("PHOS PID weight , photon  %f, pi0 %f\n\n",  fPHOSPhotonWeight,  fPHOSPi0Weight) ; 
-  printf("EMCAL PID weight, photon %f, pi0 %f\n",   fEMCALPhotonWeight,  fEMCALPi0Weight) ; 
+      //cout<<"Calo: "<<calo<<" vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
+      //  <<pmother->GetName()<<endl;
+
+      while( vy >= ipdist){//inside calorimeter
+       AliDebug(4,Form("Conversion inside calorimeter, mother vertex %0.2f, IP distance %0.2f", vy, ipdist));
+       pmother =  stack->Particle(pmother->GetMother(0));
+       vy = TMath::Abs(pmother->Vy()) ;
+       //cout<<" label "<<label<<" Mother: "<<pmother->GetName()<<" E "<<pmother->Energy()<<" Status "<<pmother->GetStatusCode()<<"  and vertex "<<vy<<endl;
+       mostatus = pmother->GetStatusCode();
+       mopdg    = TMath::Abs(pmother->GetPdgCode());
+      }//while vertex is inside calorimeter
+      //cout<<"Calo: "<<calo<<" final vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
+      //  <<pmother->GetName()<<endl;
+    }//check status and vertex
+
+    AliDebug(4,Form("%s, ESD E %2.2f, PID %d,  mother: E %2.2f, label %d, status %d,  vertex %3.2f, mother 2 %d, grandmother %d \n",
+                   calo.Data(),momentum.E(),pmother->Energy(), label, pmother->GetPdgCode(),
+                   pmother->GetStatusCode(), vy, pmother->GetMother(0), stack->GetPrimary(label)));
+    
+    //Check Decay photons
+    if(mopdg == 22){
+      
+      //his mother was a pi0?
+      TParticle * pmotherpi0 =  stack->Particle(pmother->GetMother(0));
+      if( pmotherpi0->GetPdgCode() == 111){
+
+       AliDebug(4,Form(" %s: E cluster %2.2f, E gamma %2.2f, Mother Pi0, E %0.2f, status %d, daughters %d\n",
+                       calo.Data(), momentum.E(),pmother->Energy(),pmotherpi0->Energy(), 
+                       pmotherpi0->GetStatusCode(), pmotherpi0->GetNDaughters()));
+       
+       if(pmotherpi0->GetNDaughters() == 1) mostatus = 2 ; //signal that this photon can only be decay, not overlapped.
+       else if(pmotherpi0->GetNDaughters() == 2){
+         
+         TParticle * pd1 = stack->Particle(pmotherpi0->GetDaughter(0));
+         TParticle * pd2 = stack->Particle(pmotherpi0->GetDaughter(1));
+         //if(pmotherpi0->Energy()> 10 ){
+//       cout<<"Two "<<calo<<" daugthers, pi0 label "<<pmother->GetMother(0)<<" E :"<<pmotherpi0->Energy()<<endl;
+//       cout<<" 1) label "<<pmotherpi0->GetDaughter(0)<<" pdg "<<pd1->GetPdgCode()<<" E  "<<pd1->Energy()
+//           <<" phi "<<pd1->Phi()*TMath::RadToDeg()<<" eta "<<pd1->Eta()
+//           <<" mother label "<<pd1->GetMother(0)<<" n daug "<<pd1->GetNDaughters() <<endl;
+//         cout<<" 2) label "<<pmotherpi0->GetDaughter(1)<<" pdg "<<pd2->GetPdgCode()<<" E  "<<pd2->Energy()
+//             <<" phi "<<pd2->Phi()*TMath::RadToDeg()<<" eta "<<pd2->Eta()<<" mother label "
+//             <<pd2->GetMother(0)<<" n daug "<<pd2->GetNDaughters() <<endl;
+           //}
+         if(pd1->GetPdgCode() == 22 && pd2->GetPdgCode() == 22){
+           TLorentzVector lv1 , lv2 ;
+           pd1->Momentum(lv1);
+           pd2->Momentum(lv2);
+           Double_t angle = lv1.Angle(lv2.Vect());
+//         if(pmotherpi0->Energy()> 10 )
+//           cout<<"angle*ipdist "<<angle*ipdist<<" mindist "<< minangle <<endl;
+           if (angle < minangle){
+             //There is overlapping, pass the pi0 as mother
+             
+             AliDebug(4,Form(">>>> %s cluster with E %2.2f is a overlapped pi0, E %2.2f, angle %2.4f < anglemin %2.4f\n",
+                             calo.Data(), momentum.E(), pmotherpi0->Energy(), angle, minangle));           
+            
+             pmother = pmotherpi0 ;
+             
+           }
+           else mostatus = 2 ; // gamma decay not overlapped
+         }// daughters are photons
+         else mostatus = 2; // daughters are e-gamma or e-e, no overlapped, or charged cluster
+       }//2 daughters
+       else AliDebug(4,Form("pi0 has %d daughters",pmotherpi0->GetNDaughters()));
+      }//pi0 decay photon?
+    }//photon 
+
+    pmother->SetStatusCode(mostatus); // if status = 2, decay photon.
+
+  }//label >=0
+  else AliInfo(Form("Negative Kinematic label of PHOS cluster:  %d",label));
+
+  return pmother ;
 
 }
index dcd97a9..f6c1c1d 100644 (file)
@@ -7,6 +7,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
  *
@@ -36,39 +39,15 @@ public:
   AliGammaMCDataReader & operator = (const AliGammaMCDataReader & g) ;//cpy assignment
   virtual ~AliGammaMCDataReader() {;} //virtual dtor
 
-  void InitParameters();
-
-  Bool_t   IsEMCALPIDOn() const {return fEMCALPID ; }
-  Bool_t   IsPHOSPIDOn() const {return fPHOSPID ; }
-  Float_t  GetEMCALPhotonWeight() { return  fEMCALPhotonWeight  ; }
-  Float_t  GetEMCALPi0Weight()    {  return fEMCALPi0Weight  ; }
-  Float_t  GetPHOSPhotonWeight()  {  return fPHOSPhotonWeight  ; }
-  Float_t  GetPHOSPi0Weight()  {  return fPHOSPi0Weight  ; }
-  
-  void Print(const Option_t * opt)const;
-  
-  void SetEMCALPIDOn(Bool_t pid){ fEMCALPID= pid ; }
-  void SetPHOSPIDOn(Bool_t pid){ fPHOSPID= pid ; }
-  void SetEMCALPhotonWeight(Float_t  w){  fEMCALPhotonWeight = w ; }
-  void SetEMCALPi0Weight(Float_t  w){  fEMCALPi0Weight = w ; }
-  void SetPHOSPhotonWeight(Float_t  w){  fPHOSPhotonWeight = w ; }
-  void SetPHOSPi0Weight(Float_t  w){  fPHOSPi0Weight = w ; }
-
-  void CreateParticleList(TObject * esd, TObject * stack, TClonesArray * plCh, 
-                         TClonesArray * plEMCAL, TClonesArray * plPHOS, TClonesArray *);
-
-  
- private:
-
-  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 
-  Float_t      fEMCALPi0Weight;  //Bayesian PID weight for pi0 in EMCAL 
-  Float_t      fPHOSPhotonWeight; //Bayesian PID weight for photons in PHOS 
-  Float_t      fPHOSPi0Weight; //Bayesian PID weight for pi0 in PHOS 
-
-  ClassDef(AliGammaMCDataReader,0)
+  void CreateParticleList(TObject * esd, TObject * stack, 
+                         TClonesArray * plCh, TClonesArray * plEMCAL, TClonesArray * plPHOS,
+                         TClonesArray * plPrimCh, TClonesArray * plPrimEMCAL, TClonesArray * plPrimPHOS);
+
+  TParticle *  GetMotherParticle(Int_t label, AliStack *stack,  TString calo, TLorentzVector momentum) ;
+
+  private:
+
+  ClassDef(AliGammaMCDataReader,1)
 } ;
  
 
index a593f2e..8195314 100644 (file)
@@ -9,7 +9,7 @@
  * without fee, provided that the above copyright notice appears in all   *
  * copies and that both the copyright notice and this permission notice   *
  * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
+ * about the suitability of this software for any purpose. It is         *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 /* $Id$ */
@@ -17,9 +17,6 @@
 /* History of cvs commits:
  *
  * $Log$
- * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
- * new analysis classes in the the new analysis framework
- *
  *
  */
 
@@ -49,26 +46,19 @@ ClassImp(AliGammaMCReader)
 
 //____________________________________________________________________________
 AliGammaMCReader::AliGammaMCReader() : 
-  AliGammaReader(),
-  fEMCALIPDistance(0.),   fPHOSIPDistance(0.), 
-  fEMCALMinDistance(0.),   fPHOSMinDistance(0.), 
-  fDecayPi0(kFALSE)
+  AliGammaReader(), fDecayPi0(0), fCheckOverlapping(kFALSE)
 {
   //Ctor
   
   //Initialize parameters
-  fDataType = kMC;
   InitParameters();
-  
+  fDataType = kMC;  
   
 }
 
 //____________________________________________________________________________
 AliGammaMCReader::AliGammaMCReader(const AliGammaMCReader & g) :   
-  AliGammaReader(g),
-  fEMCALIPDistance(g.fEMCALIPDistance),  fPHOSIPDistance(g.fPHOSIPDistance),  
-  fEMCALMinDistance(g.fEMCALMinDistance),  fPHOSMinDistance(g.fPHOSMinDistance), 
-  fDecayPi0(g.fDecayPi0)
+  AliGammaReader(g), fDecayPi0(g.fDecayPi0), fCheckOverlapping(g.fCheckOverlapping)
 {
   // cpy ctor
 }
@@ -80,24 +70,102 @@ AliGammaMCReader & AliGammaMCReader::operator = (const AliGammaMCReader & source
 
   if(&source == this) return *this;
 
-  fEMCALIPDistance = source.fEMCALIPDistance; 
-  fPHOSIPDistance = source.fPHOSIPDistance; 
-  fEMCALMinDistance = source.fEMCALMinDistance; 
-  fPHOSMinDistance = source.fPHOSMinDistance; 
-  fDecayPi0 = source.fDecayPi0;
+  fDecayPi0 = source.fDecayPi0; 
+  fCheckOverlapping = source.fCheckOverlapping;
 
   return *this;
 
 }
 
+//___________________________________________________________________________
+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.
+      FillListWithDecayGammaOrPi0(pmother, pdaug0, pdaug1, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
+    
+  }//mother is a pi0 with 2 daughters
+  else{
+    if(IsInPHOS(particle->Phi(),particle->Eta()))
+      new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
+    else if(IsInEMCAL(particle->Phi(),particle->Eta()))
+      new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
+  }
+
+}
+
+//___________________________________________________________________________
+ void AliGammaMCReader::CaseGeantDecay(TParticle * particle, AliStack * stack,
+                                  TClonesArray * plEMCAL, Int_t &indexEMCAL,
+                                  TClonesArray * plPHOS, Int_t &indexPHOS){
+   //Find decay gamma from pi0, decayed by GEANT and put them in the list.
+   
+   Int_t ndaug = particle->GetNDaughters() ;
+   TParticle * d1 = new TParticle();
+   TParticle * d2 = new TParticle();
+   
+   if(ndaug > 0){//At least 1 daugther
+     d1 = stack->Particle(particle->GetDaughter(0));
+     if (ndaug > 1 ) 
+       d2 = stack->Particle(particle->GetDaughter(1));
+     
+     FillListWithDecayGammaOrPi0(particle, d1, d2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
+     
+   }// ndaugh > 0
+ }
+
+//___________________________________________________________________________
+void AliGammaMCReader::CasePi0Decay(TParticle * pPi0, 
+                                   TClonesArray * plEMCAL, Int_t &indexEMCAL,
+                                   TClonesArray * plPHOS, Int_t &indexPHOS){
+  
+  //Decays pi0, see if aperture angle is small and then add the pi0 or the 2 gamma
+  
+  TLorentzVector lvPi0, lvGamma1, lvGamma2 ;
+  //Double_t angle = 0;
+  
+  lvPi0.SetPxPyPzE(pPi0->Px(),pPi0->Py(),pPi0->Pz(),pPi0->Energy());
+
+  //Decay
+  MakePi0Decay(lvPi0,lvGamma1,lvGamma2);//,angle);
+  
+  //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
+  TParticle * pPhoton1 = new TParticle(22,1,0,0,0,0,lvGamma1.Px(),lvGamma1.Py(),
+                                     lvGamma1.Pz(),lvGamma1.E(),0,0,0,0);   
+  TParticle * pPhoton2 = new TParticle(22,1,0,0,0,0,lvGamma2.Px(),lvGamma2.Py(),
+                                     lvGamma2.Pz(),lvGamma2.E(),0,0,0,0);
+
+  FillListWithDecayGammaOrPi0(pPi0, pPhoton1, pPhoton2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
+}
+
+//_______________________________________________________________
+void AliGammaMCReader::InitParameters()
+{
+  
+  //Initialize the parameters of the analysis.
+
+  fDecayPi0 = kGeantDecay;
+  fCheckOverlapping = kTRUE ;
+}
 
 //____________________________________________________________________________
 void AliGammaMCReader::CreateParticleList(TObject * data, TObject *, 
-                                        TClonesArray * plCh, 
-                                        TClonesArray * plEMCAL,  
+                                         TClonesArray * plCh, 
+                                         TClonesArray * plEMCAL,  
                                          TClonesArray * plPHOS,
-                                         TClonesArray *plParton)
+                                         TClonesArray *plParton,TClonesArray *,TClonesArray *)
 {
   //Create list of particles from EMCAL, PHOS and CTS. 
   AliStack * stack = (AliStack *) data ;
@@ -111,13 +179,14 @@ void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
   for (iParticle=0 ; iParticle <  stack->GetNprimary() ; iParticle++) {
     TParticle * particle = stack->Particle(iParticle); 
     
+
     //Keep partons
     if(particle->GetStatusCode() == 21 && iParticle>=2){//All partons, not nucleus
       new((*plParton)[indexParton++])  TParticle(*particle) ;
     }
 
     //Keep Stable particles 
-    if((particle->GetStatusCode() == 0) && (particle->Pt() > 0)){
+    if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
       
       charge = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode())->Charge();
       
@@ -134,10 +203,17 @@ void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
              TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos
        
        if(particle->GetPdgCode()!=111){
-         if(IsInPHOS(particle->Phi(),particle->Eta()))
-           new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
-         else if(IsInEMCAL(particle->Phi(),particle->Eta()))
-           new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
+         //In case that we force PYTHIA to decay pi0, and we want to check the overlapping of 
+         // the decay gamma.
+         if(particle->GetPdgCode() == 22 && fDecayPi0==kDecayGamma){
+           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) ;
+         }
        }//no pi0
        else{
          if(fDecayPi0 == kNoDecay){//keep the pi0 do not decay
@@ -147,15 +223,64 @@ void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
              new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
          }
          else if(fDecayPi0 == kDecay)
-           MakePi0Decay(particle,plEMCAL,indexEMCAL,plPHOS, indexPHOS);
+           CasePi0Decay(particle,plEMCAL,indexEMCAL,plPHOS, indexPHOS);
          else if(fDecayPi0 == kGeantDecay)
-           SetGeantDecay(particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS);
+           CaseGeantDecay(particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS);
        }//pi0  
       }//neutral particle
     }//stable particle
   }//particle loop
 }
 
+
+//___________________________________________________________________________
+void AliGammaMCReader::FillListWithDecayGammaOrPi0(TParticle * pPi0, TParticle * pdaug0, TParticle * pdaug1,
+                                  TClonesArray * plEMCAL, Int_t &indexEMCAL,
+                                  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 (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;
+      }
+    }
+  }//fCheckOverlapping
+  
+  //Fill with gammas if not overlapp
+  if(!overlap){
+    if(pdaug0->GetPdgCode() == 22 || TMath::Abs(pdaug0->GetPdgCode() ) == 11 ){
+      if(IsInEMCAL(pdaug0->Phi(),pdaug0->Eta()) &&  pdaug0->Pt() > fNeutralPtCut)
+       new((*plEMCAL)[indexEMCAL++])       TParticle(*pdaug0) ;
+      else if(IsInPHOS(pdaug0->Phi(),pdaug0->Eta()) &&  pdaug0->Pt() > fNeutralPtCut)
+       new((*plPHOS)[indexPHOS++])       TParticle(*pdaug0) ;
+    }
+    
+    if(pdaug1->GetPdgCode() == 22 || TMath::Abs(pdaug1->GetPdgCode() ) == 11 ){
+      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){
   //Check if particle is in EMCAL acceptance
@@ -180,113 +305,9 @@ Bool_t  AliGammaMCReader::IsInPHOS(Double_t phi, Double_t eta){
   return kFALSE ;
 }
 
-//___________________________________________________________________________
-void AliGammaMCReader::SetGeantDecay(TParticle * particle, AliStack * stack,
-                                    TClonesArray * plEMCAL, Int_t &indexEMCAL,
-                                    TClonesArray * plPHOS, Int_t &indexPHOS){
-  //Find decay gamma from pi0 and put them in the list.
-  
-  Int_t ndaug = particle->GetNDaughters() ;
-  if(ndaug<=2 && ndaug >0){//At least 1 daugther
-    TParticle * d1 = stack->Particle(particle->GetDaughter(0));
-    if(d1->GetPdgCode()==22){
-      if(IsInEMCAL(d1->Phi(),d1->Eta()))
-       new((*plEMCAL)[indexEMCAL++])       TParticle(*d1) ;
-      else if(IsInPHOS(d1->Phi(),d1->Eta()))
-       new((*plPHOS)[indexPHOS++])       TParticle(*d1) ;
-    }
-    
-    if(ndaug>1){//second daugther if present
-      TParticle * d2 = stack->Particle(particle->GetDaughter(1));
-      if(IsInEMCAL(d2->Phi(),d2->Eta()))
-       new((*plEMCAL)[indexEMCAL++])       TParticle(*d2) ;
-      else if(IsInPHOS(d2->Phi(),d2->Eta()))
-       new((*plPHOS)[indexPHOS++])       TParticle(*d2) ;
-    }
-  }
-}
-
-//___________________________________________________________________________
-void AliGammaMCReader::MakePi0Decay(TParticle * particle, 
-                                   TClonesArray * plEMCAL, Int_t &indexEMCAL,
-                                   TClonesArray * plPHOS, Int_t &indexPHOS){
-  
-  //Decays pi0, see if aperture angle is small and then add the pi0 or the 2 gamma
-  
-  TLorentzVector pPi0, pGamma1, pGamma2 ;
-  Double_t angle = 0, cellDistance = 0.;
-  Bool_t checkPhoton = kTRUE;
-  
-  pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
-  
-  //Decay
-  Pi0Decay(pPi0,pGamma1,pGamma2,angle);
-  
-  //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
-  if(IsInPHOS(particle->Phi(), particle->Eta())){
-    cellDistance = angle*fPHOSIPDistance;
-    if (cellDistance < fPHOSMinDistance){
-      new((*plPHOS)[indexPHOS++])       TParticle(*particle) ;
-      checkPhoton = kFALSE;
-    }
-  } 
-  else if(IsInEMCAL(particle->Phi(), particle->Eta())){
-    cellDistance = angle*fEMCALIPDistance;
-    if (cellDistance < fEMCALMinDistance) {
-      new((*plEMCAL)[indexEMCAL++])       TParticle(*particle) ;
-      checkPhoton = kFALSE;
-    }
-  } 
-  else checkPhoton = kTRUE ;
-  
-  if (checkPhoton) {
-    //Gamma Not overlapped
-    TParticle * photon1 = new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
-                                       pGamma1.Pz(),pGamma1.E(),0,0,0,0);    
-    if(photon1->Pt() > fNeutralPtCut) {
-      
-      if(IsInPHOS(photon1->Phi(), photon1->Eta()))
-       new((*plPHOS)[indexPHOS++])       TParticle(*photon1) ;
-      
-      else if(IsInEMCAL(photon1->Phi(), photon1->Eta()))
-       new((*plEMCAL)[indexEMCAL++])       TParticle(*photon1) ;
-      
-    }// photon 1 of pi0 in acceptance
-    
-    
-    TParticle * photon2 = new TParticle(22,1,0,0,0,0,pGamma2.Px(),pGamma2.Py(),
-                                       pGamma2.Pz(),pGamma2.E(),0,0,0,0);
-    
-    if(photon2->Pt() > fNeutralPtCut) {
-      
-      if(IsInPHOS(photon2->Phi(), photon2->Eta()))
-       new((*plPHOS)[indexPHOS++])       TParticle(*photon2) ;
-      
-      else if(IsInEMCAL(photon2->Phi(), photon2->Eta()))
-       new((*plEMCAL)[indexEMCAL++])       TParticle(*photon2) ;
-      
-    }// photon 2 of pi0 in acceptance
-  }//Not overlapped gamma
-}
-
-
-//_______________________________________________________________
-void AliGammaMCReader::InitParameters()
-{
-  
-  //Initialize the parameters of the analysis.
-  
-  //Fill particle lists when PID is ok
-  fEMCALMinDistance    = 10. ;  
-  fPHOSMinDistance    = 3.6 ;
-  fEMCALIPDistance    = 450. ;//cm  
-  fPHOSIPDistance    = 460. ;//cm
-  fDecayPi0 = kGeantDecay;
-}
-
 //____________________________________________________________________________
-void AliGammaMCReader::Pi0Decay(TLorentzVector &p0, TLorentzVector &p1, 
-                               TLorentzVector &p2, Double_t &angle) {
+void AliGammaMCReader::MakePi0Decay(TLorentzVector &p0, TLorentzVector &p1, 
+                               TLorentzVector &p2){//, Double_t &angle) {
   // Perform isotropic decay pi0 -> 2 photons
   // p0 is pi0 4-momentum (inut)
   // p1 and p2 are photon 4-momenta (output)
@@ -321,7 +342,7 @@ void AliGammaMCReader::Pi0Decay(TLorentzVector &p0, TLorentzVector &p1,
   p2.Boost(b);
   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
   //cout<<"angle"<<endl;
-  angle = p1.Angle(p2.Vect());
+  //angle = p1.Angle(p2.Vect());
   //cout<<angle<<endl;
 }
 
@@ -335,10 +356,11 @@ void AliGammaMCReader::Print(const Option_t * opt) const
   
   Info("Print", "%s %s", GetName(), GetTitle() ) ;
   
-  printf("IP distance to PHOS         : %f\n", fPHOSIPDistance) ;
-  printf("IP distance to EMCAL         : %f\n", fEMCALIPDistance) ;
-  printf("Min gamma decay distance in PHOS         : %f\n", fPHOSMinDistance) ;
-  printf("Min gamma decay distance in EMCAL         : %f\n", fEMCALMinDistance) ;
   printf("Decay Pi0?          : %d\n", fDecayPi0) ;
-  
+  printf("Check Overlapping?          : %d\n", fCheckOverlapping) ;
+
 }
+
+
+
+
index b1ab3b6..cfbade3 100644 (file)
@@ -7,9 +7,6 @@
 /* History of cvs commits:
  *
  * $Log$
- * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
- * new analysis classes in the the new analysis framework
- *
  *
  */
 
@@ -37,7 +34,7 @@ public:
   AliGammaMCReader & operator = (const AliGammaMCReader & g) ;//cpy assignment
   virtual ~AliGammaMCReader() {;} //virtual dtor
 
-  enum decay_t {kNoDecay, kGeantDecay, kDecay};
+  enum decay_t {kNoDecay, kGeantDecay, kDecay, kDecayGamma};
  
   void InitParameters();
 
@@ -45,44 +42,45 @@ public:
   Bool_t  IsInPHOS(Double_t phi, Double_t eta) ;
 
   Int_t    GetDecayPi0Flag() const {return fDecayPi0 ; }
-  Float_t  GetEMCALIPDistance()  {  return fEMCALIPDistance ; }
-  Float_t  GetPHOSIPDistance()  {  return fPHOSIPDistance ; }
-  Float_t  GetEMCALMinDistance()  {  return fEMCALMinDistance ; }
-  Float_t  GetPHOSMinDistance()  {  return fPHOSMinDistance ; }
 
   void Print(const Option_t * opt)const;
   
   void SetDecayPi0Flag(Int_t d){ fDecayPi0 = d ; }
-  void SetEMCALIPDistance(Float_t  d){  fEMCALIPDistance = d ; }
-  void SetPHOSIPDistance(Float_t  d){  fPHOSIPDistance = d ; }
-  void SetEMCALMinDistance(Float_t  d){  fEMCALMinDistance = d ; }
-  void SetPHOSMinDistance(Float_t  d){  fPHOSMinDistance = d ; }
 
+  void SetCheckOverlapping(Bool_t check){fCheckOverlapping = check ;}
+  Bool_t IsCheckOverlappingOn() {return fCheckOverlapping ;}
+
+  private:
+  
+  void CaseDecayGamma(Int_t index, TParticle * particle, AliStack * stack,
+                     TClonesArray * plEMCAL, Int_t &indexEMCAL,
+                     TClonesArray * plPHOS, Int_t &indexPHOS);
+  
+  void CaseGeantDecay(TParticle * particle, AliStack * stack,
+                     TClonesArray * plEMCAL, Int_t &indexEMCAL,
+                     TClonesArray * plPHOS, Int_t &indexPHOS);
+  
+  void CasePi0Decay(TParticle * particle, TClonesArray * plEMCAL, Int_t &indexEMCAL,
+                   TClonesArray * plPHOS, Int_t &indexPHOS);
+  
   void CreateParticleList(TObject * stack, TObject * ,
                          TClonesArray * plCh, TClonesArray * plEMCAL, 
-                         TClonesArray * plPHOS, TClonesArray * plParton);
-  void MakePi0Decay(TParticle * particle, TClonesArray * plEMCAL, Int_t &indexEMCAL,
-                   TClonesArray * plPHOS, Int_t &indexPHOS);
+                         TClonesArray * plPHOS, TClonesArray * plParton,TClonesArray *,TClonesArray *);
+  
+  void FillListWithDecayGammaOrPi0(TParticle * pPi0, TParticle * pdaug0, TParticle * pdaug1,
+                                  TClonesArray * plEMCAL, Int_t &indexEMCAL,
+                                  TClonesArray * plPHOS, Int_t &indexPHOS);  
+  void MakePi0Decay(TLorentzVector &p0, TLorentzVector &p1, TLorentzVector &p2);//, Double_t &angle);
+  
+  
+  private:
   
-  void Pi0Decay(TLorentzVector &p0, TLorentzVector &p1, TLorentzVector &p2, 
-               Double_t &angle);
-
-  void SetGeantDecay(TParticle * particle, AliStack * stack,
-                    TClonesArray * plEMCAL, Int_t &indexEMCAL,
-                    TClonesArray * plPHOS, Int_t &indexPHOS);
- private:
-
-  Float_t      fEMCALIPDistance; //Calorimeter IP distance.
-  Float_t      fPHOSIPDistance; //Calorimeter IP distance
-  Float_t      fEMCALMinDistance; //Gamma decay minimum aperture.
-  Float_t      fPHOSMinDistance; //Gamma decay minimum aperture.
-
   Int_t      fDecayPi0; //Decay Pi0.
-
-  ClassDef(AliGammaMCReader,0)
+  Bool_t   fCheckOverlapping; // if True, check if gammas from decay overlapp in calorimeters.
+  
+  ClassDef(AliGammaMCReader,1)
 } ;
+
 
 #endif //ALIGAMMAMCREADER_H
 
index 209e570..8a42521 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
  *
@@ -43,8 +46,15 @@ ClassImp(AliGammaReader)
 AliGammaReader::AliGammaReader() : 
   TObject(), fDataType(0),
   fCTSEtaCut(0.), fEMCALEtaCut(0.), fPHOSEtaCut(0.),
-  fNeutralPtCut(0.),
-  fChargedPtCut(0.)
+  fNeutralPtCut(0.), fChargedPtCut(0.),
+  fEMCALIPDistance(0.),   fPHOSIPDistance(0.), 
+  fEMCALMinAngle(0.),   fPHOSMinAngle(0.), 
+  fEMCALPID(0),fPHOSPID(0),
+  fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),  fEMCALElectronWeight(0.),  
+  fEMCALChargeWeight(0.),fEMCALNeutralWeight(0.),
+  fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),  fPHOSElectronWeight(0.), 
+  fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.) ,
+  fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0) 
 {
   //Ctor
 
@@ -61,8 +71,24 @@ AliGammaReader::AliGammaReader() :
 AliGammaReader::AliGammaReader(const AliGammaReader & g) :   
   TObject(g), fDataType(g.fDataType),
   fCTSEtaCut(g.fCTSEtaCut),  fEMCALEtaCut(g.fEMCALEtaCut),  fPHOSEtaCut(g.fPHOSEtaCut),
-  fNeutralPtCut(g.fNeutralPtCut),
-  fChargedPtCut(g.fChargedPtCut)
+  fNeutralPtCut(g.fNeutralPtCut), fChargedPtCut(g.fChargedPtCut),
+  fEMCALIPDistance(g.fEMCALIPDistance),  fPHOSIPDistance(g.fPHOSIPDistance),  
+  fEMCALMinAngle(g.fEMCALMinAngle),  fPHOSMinAngle(g.fPHOSMinAngle), 
+  fEMCALPID(g.fEMCALPID), 
+  fPHOSPID(g.fPHOSPID),
+  fEMCALPhotonWeight(g.fEMCALPhotonWeight), 
+  fEMCALPi0Weight(g.fEMCALPi0Weight), 
+  fEMCALElectronWeight(g.fEMCALElectronWeight), 
+  fEMCALChargeWeight(g.fEMCALChargeWeight), 
+  fEMCALNeutralWeight(g.fEMCALNeutralWeight), 
+  fPHOSPhotonWeight(g.fPHOSPhotonWeight),
+  fPHOSPi0Weight(g.fPHOSPi0Weight),
+  fPHOSElectronWeight(g.fPHOSElectronWeight), 
+  fPHOSChargeWeight(g.fPHOSChargeWeight),
+  fPHOSNeutralWeight(g.fPHOSNeutralWeight),
+  fPHOSWeightFormula(g.fPHOSWeightFormula), 
+  fPHOSPhotonWeightFormula(g.fPHOSPhotonWeightFormula), 
+  fPHOSPi0WeightFormula(g.fPHOSPi0WeightFormula) 
 {
   // cpy ctor
 
@@ -91,6 +117,31 @@ AliGammaReader & AliGammaReader::operator = (const AliGammaReader & source)
   fPhiPHOSCut[0]=source.fPhiPHOSCut[0];
   fPhiPHOSCut[1]=source.fPhiPHOSCut[1];
 
+  fEMCALIPDistance = source.fEMCALIPDistance; 
+  fPHOSIPDistance = source.fPHOSIPDistance; 
+  fEMCALMinAngle = source.fEMCALMinAngle; 
+  fPHOSMinAngle = source.fPHOSMinAngle; 
+
+  fEMCALPID = source.fEMCALPID ;
+  fPHOSPID = source.fPHOSPID ;
+
+  fEMCALPhotonWeight = source. fEMCALPhotonWeight ;
+  fEMCALPi0Weight = source.fEMCALPi0Weight ;
+  fEMCALElectronWeight = source.fEMCALElectronWeight; 
+  fEMCALChargeWeight = source.fEMCALChargeWeight;
+  fEMCALNeutralWeight = source.fEMCALNeutralWeight;
+
+  fPHOSPhotonWeight = source.fPHOSPhotonWeight ;
+  fPHOSPi0Weight = source.fPHOSPi0Weight ;
+  fPHOSElectronWeight = source.fPHOSElectronWeight; 
+  fPHOSChargeWeight = source.fPHOSChargeWeight;
+  fPHOSNeutralWeight = source.fPHOSNeutralWeight;
+
+  fPHOSWeightFormula       = source.fPHOSWeightFormula; 
+  fPHOSPhotonWeightFormula = source.fPHOSPhotonWeightFormula; 
+  fPHOSPi0WeightFormula    = source.fPHOSPi0WeightFormula;
+
+
   return *this;
 
 }
@@ -111,6 +162,34 @@ void AliGammaReader::InitParameters()
   fNeutralPtCut   = 0.5 ;
   fChargedPtCut   = 0.5 ;
 
+  fEMCALMinAngle    =   2.5 * TMath::DegToRad() ;  
+  fPHOSMinAngle    = 0.45 * TMath::DegToRad() ; //3.6 ;
+  fEMCALIPDistance    = 450. ;//cm  
+  fPHOSIPDistance    = 445. ;//cm 460 (EMCA) - 15 (CPV)
+
+  //pid, only for ESD data
+  fEMCALPID = kFALSE;
+  fPHOSPID = kFALSE;
+
+  fEMCALPhotonWeight = 0.8 ;
+  fEMCALPi0Weight = 0.5 ;
+  fEMCALElectronWeight = 0.8 ;
+  fEMCALChargeWeight = 0.5 ;
+  fEMCALNeutralWeight = 0.5 ;
+
+  fPHOSPhotonWeight = 0.75 ;
+  fPHOSPi0Weight = 0.8 ;
+  fPHOSElectronWeight = 0.5 ;
+  fPHOSChargeWeight = 0.5 ;
+  fPHOSNeutralWeight = 0.5 ;
+
+  //Formula to set the PID weight threshold for photon or pi0
+  fPHOSWeightFormula = kTRUE;
+  fPHOSPhotonWeightFormula = 
+    new TFormula("photonWeight","0.98*(x<40)+ 0.68*(x>=100)+(x>=40 && x<100)*(0.98+x*(6e-3)-x*x*(2e-04)+x*x*x*(1.1e-06))");
+  fPHOSPi0WeightFormula = 
+    new TFormula("pi0Weight","0.98*(x<65)+ 0.915*(x>=100)+(x>=65 && x-x*(1.95e-3)-x*x*(4.31e-05)+x*x*x*(3.61e-07))");
+
 }
 
 
@@ -132,6 +211,32 @@ void AliGammaReader::Print(const Option_t * opt) const
   printf("pT neutral cut           : %f GeV/c\n", fNeutralPtCut) ;
   printf("pT charged cut           : %f GeV/c\n", fChargedPtCut) ;
 
+  if(fDataType == kMC || fDataType == kMCData){
+    printf("IP distance to PHOS         : %f\n", fPHOSIPDistance) ;
+    printf("IP distance to EMCAL         : %f\n", fEMCALIPDistance) ;
+    printf("Min gamma decay aperture angle in PHOS         : %f\n", fPHOSMinAngle) ;
+    printf("Min gamma decay aperture angle in EMCAL         : %f\n", fEMCALMinAngle) ;
+  }
+  
+  if(fDataType != kMC){
+    printf("PHOS PID on?               =     %d\n",  fPHOSPID) ; 
+    printf("EMCAL PID  on?         =     %d\n",  fEMCALPID) ;
+    printf("PHOS PID weight , photon %f, pi0 %f, e %f, charge %f, neutral %f \n",  
+          fPHOSPhotonWeight,  fPHOSPi0Weight, 
+          fPHOSElectronWeight,  fPHOSChargeWeight,   fPHOSNeutralWeight) ; 
+    printf("EMCAL PID weight, photon %f, pi0 %f, e %f, charge %f, neutral %f\n",   
+          fEMCALPhotonWeight,  fEMCALPi0Weight, 
+          fEMCALElectronWeight,  fEMCALChargeWeight,  fEMCALNeutralWeight) ; 
+    
+    printf("PHOS Parametrized weight on?               =     %d\n",  fPHOSWeightFormula) ; 
+    if(fPHOSWeightFormula){
+      printf(">>>>>>>>>>> Photon weight formula<<<<<<<<<<<<\n");
+      fPHOSPhotonWeightFormula->Print();
+      printf(">>>>>>>>>>> Pi0    weight formula<<<<<<<<<<<<\n");
+      fPHOSPhotonWeightFormula->Print();
+    }
+  }
+
 } 
 
 
index 6ae0f4c..08fcb66 100644 (file)
@@ -7,6 +7,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
  *
@@ -19,6 +22,7 @@
 
 // --- ROOT system ---
 #include <TParticle.h> 
+#include <TFormula.h>
 #include <TClonesArray.h> 
 #include "AliStack.h"
 #include "TObject.h" 
@@ -36,6 +40,18 @@ public:
   AliGammaReader & operator = (const AliGammaReader & g) ;//cpy assignment
   virtual ~AliGammaReader() {;} //virtual dtor
 
+  enum PidType {
+    kPhoton = 22,
+    kPi0 = 111,
+    kEta = 221, 
+    kElectron = 11, 
+    kEleCon = 13, 
+    kNeutralHadron = 2112, 
+    kChargedHadron = 211, 
+    kNeutralUnknown = 130, 
+    kChargedUnknown=321
+  };
+
   enum datatype_t {kData, kMC, kMCData};
   
   void InitParameters();
@@ -51,6 +67,28 @@ public:
   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 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 void Print(const Option_t * opt)const;
   
   virtual void SetCTSEtaCut(Float_t eta){ fCTSEtaCut= eta ; }
@@ -63,9 +101,31 @@ public:
   virtual void SetNeutralPtCut(Float_t  pt){  fNeutralPtCut = pt ; }
   virtual void SetChargedPtCut(Float_t  pt){  fChargedPtCut = pt ; }
 
+ virtual void SetEMCALIPDistance(Float_t  d){  fEMCALIPDistance = d ; }
+  virtual void SetPHOSIPDistance(Float_t  d){  fPHOSIPDistance = d ; }
+  virtual void SetEMCALMinAngle(Float_t  d){  fEMCALMinAngle = d ; }
+  virtual void SetPHOSMinAngle(Float_t  d){  fPHOSMinAngle = d ; }
+  
+  virtual void SetEMCALPIDOn(Bool_t pid){ fEMCALPID= pid ; }
+  virtual void SetPHOSPIDOn(Bool_t pid){ fPHOSPID= pid ; }
+  virtual void SetEMCALPhotonWeight(Float_t  w){  fEMCALPhotonWeight = w ; }
+  virtual void SetEMCALPi0Weight(Float_t  w){  fEMCALPi0Weight = w ; }
+  virtual void SetEMCALElectronWeight(Float_t  w){  fEMCALElectronWeight = w ; }
+  virtual void SetEMCALChargeWeight(Float_t  w){  fEMCALChargeWeight = w ; }
+  virtual void SetEMCALNeutralWeight(Float_t  w){  fEMCALNeutralWeight = w ; }
+  virtual void SetPHOSPhotonWeight(Float_t  w){  fPHOSPhotonWeight = w ; }
+  virtual void SetPHOSPi0Weight(Float_t  w){  fPHOSPi0Weight = w ; }
+  virtual void SetPHOSElectronWeight(Float_t  w){  fPHOSElectronWeight = w ; }
+  virtual void SetPHOSChargeWeight(Float_t  w){  fPHOSChargeWeight = w ; }
+  virtual void SetPHOSNeutralWeight(Float_t  w){  fPHOSNeutralWeight = w ; }
+
+  virtual void UsePHOSPIDWeightFormula(Bool_t par)  { fPHOSWeightFormula  = par; } 
+  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 * plCh, TClonesArray * plEMCAL, TClonesArray * plPHOS,     
+                                 TClonesArray * parton, TClonesArray * plPrimEMCAL, TClonesArray * plPrimPHOS) {;}
  protected:
   Int_t        fDataType ;
   Float_t      fCTSEtaCut ;//CTS  pseudorapidity acceptance
@@ -76,9 +136,31 @@ public:
   Float_t      fNeutralPtCut; //
   Float_t      fChargedPtCut;  // 
 
-   ClassDef(AliGammaReader,0)
+  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.
+
+  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 
+  Float_t      fEMCALPi0Weight;  //Bayesian PID weight for pi0 in EMCAL 
+  Float_t      fEMCALElectronWeight; //Bayesian PID weight for electrons in EMCAL 
+  Float_t      fEMCALChargeWeight;  //Bayesian PID weight for charged hadrons in EMCAL 
+  Float_t      fEMCALNeutralWeight;  //Bayesian PID weight for neutral hadrons in EMCAL 
+  Float_t      fPHOSPhotonWeight; //Bayesian PID weight for photons in PHOS 
+  Float_t      fPHOSPi0Weight; //Bayesian PID weight for pi0 in PHOS 
+  Float_t      fPHOSElectronWeight; //Bayesian PID weight for electrons in PHOS 
+  Float_t      fPHOSChargeWeight; //Bayesian PID weight for charged hadrons in PHOS 
+  Float_t      fPHOSNeutralWeight; //Bayesian PID weight for neutral hadrons in PHOS 
+
+  Bool_t  fPHOSWeightFormula ; //Use parametrized weight threshold, function of energy
+  TFormula * fPHOSPhotonWeightFormula ; //Formula for photon weight
+  TFormula * fPHOSPi0WeightFormula ; //Formula for pi0 weight
+  
+  ClassDef(AliGammaReader,1)
 } ;
+
 
 #endif //ALIGAMMAREADER_H
 
index b27f41a..73421fe 100644 (file)
@@ -7,6 +7,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
  *
@@ -98,7 +101,7 @@ class AliNeutralMesonSelection : public TObject {
   TH2F * fhInvMassPairOpeningAngleCut  ; 
   TH2F * fhInvMassPairAllCut   ; 
   
-  ClassDef(AliNeutralMesonSelection,0)
+  ClassDef(AliNeutralMesonSelection,1)
     
     } ;
 
index adc893c..13d3e35 100644 (file)
@@ -13,8 +13,8 @@ AliAnaGamma*  ConfigGammaAnalysis()
     //----------------------------------------------------------
 
     // -----Option 1------ Data, ESDs
-    AliGammaDataReader *reader = new AliGammaDataReader();
-//     //AliGammaMCDataReader *reader = new AliGammaMCDataReader(); //Copy of AliGammaDataReader for the moment
+    //AliGammaDataReader *reader = new AliGammaDataReader();
+    //       AliGammaMCDataReader *reader = new AliGammaMCDataReader(); //Copy of AliGammaDataReader + Montecarlo Information
 //     //Set detectors acceptance (in real and montecarlo data)
 //     reader->SetCTSEtaCut(1); reader->SetEMCALEtaCut(1); reader->SetPHOSEtaCut(0.2);
 //     reader->SetPhiEMCALCut(40*TMath::DegToRad(), 200*TMath::DegToRad());     
@@ -22,26 +22,33 @@ AliAnaGamma*  ConfigGammaAnalysis()
 //     //Set minimum pt for particles in analysis  (in real and montecarlo data)
 //     reader->SetNeutralPtCut(0.4); reader->SetChargedPtCut(0.4); 
 //     //pid of measured particles
-//     reader->SetEMCALPIDOn(kFALSE); reader->SetPHOSPIDOn(kFALSE); //No pid, accept all particles 
+    //reader->SetEMCALPIDOn(kTRUE); reader->SetPHOSPIDOn(kTRUE); //No pid, accept all particles 
 //     //if previous kTrue
 //     reader->SetPHOSPhotonWeight(0.7);    reader->SetPHOSPi0Weight(0.7); 
 //     reader->SetEMCALPhotonWeight(0.7);    reader->SetEMCALPi0Weight(0.7);
-   
+    //reader->UsePHOSPIDWeightFormula(kTRUE);
+    //TFormula * photonF = new TFormula("photonWeight","0.98*(x<40)+ 0.68*(x>=100)+(x>=40 && x<100)*(0.98+x*(6e-3)-x*x*(2e-04)+x*x*x*(1.1e-06))");
+    //TFormula * pi0F = new TFormula("pi0Weight","0.98*(x<65)+ 0.915*(x>=100)+(x>=65 && x-x*(1.95e-3)-x*x*(4.31e-05)+x*x*x*(3.61e-07))");
+    //reader->SetPHOSPhotonWeightFormula(photonF);
+    //reader->SetPHOSPi0WeightFormula(pi0F);
+
      // -----Option 2------ Kinematics
-    //    AliGammaMCReader *reader = new AliGammaMCReader();
+       AliGammaMCReader *reader = new AliGammaMCReader();
 //     //Set detectors acceptance (in real and montecarlo data)
 //     reader->SetCTSEtaCut(1); reader->SetEMCALEtaCut(1); reader->SetPHOSEtaCut(0.2);
 //     reader->SetPhiEMCALCut(40*TMath::DegToRad(), 200*TMath::DegToRad());     
 //     reader->SetPhiPHOSCut(200*TMath::DegToRad(), 350*TMath::DegToRad()); 
 //     //Set minimum pt for particles in analysis  (in real and montecarlo data)
 //     reader->SetNeutralPtCut(0.4); reader->SetChargedPtCut(0.4); 
-//     reader->SetDecayPi0Flag(AliGammaMCReader::kNoDecay) ; //Options
+     reader->SetDecayPi0Flag(AliGammaMCReader::kGeantDecay) ; //Options
 //     //kNoDecay :Do not decay pi0, keep them in the list as they are
 //     //kGeantDecay: Look for gamma decayed by GEANT
 //     //kDecay: Decay pi0 by hand (geant was not used)
+//     //kDecayGamma: Pi0 is decayed by PYTHIA, pi0 is not final check if photons overlapp
 //     //parameters to study if decay is overlapped:
 //     reader->SetEMCALIPDistance(450.); reader->SetPHOSIPDistance(460.);
 //     reader->SetEMCALMinDistance(3.6);    reader->SetPHOSMinDistance(11.); //Miimum overlapp distance
+       reader->SetCheckOverlapping(kTRUE);
  
     //----------------------------------------------------------
     //----------------------------------------------------
@@ -65,16 +72,16 @@ AliAnaGamma*  ConfigGammaAnalysis()
     AliAnaGammaDirect *gd = new AliAnaGammaDirect();
     gd->SetMinGammaPt(1.);
     //not used in option kIsolationCut
-    gd->SetConeSize(0.3); gd->SetPtThreshold(1.); gd->SetPtSumThreshold(10.);
-    gd->SetICMethod(AliAnaGammaDirect::kNoIC) ;//Options:
+    //gd->SetConeSize(0.5); gd->SetPtThreshold(0.); gd->SetPtSumThreshold(0.);
+    gd->SetICMethod(AliAnaGammaDirect::kPtIC) ;//Options:
           //kNoIC: Accept all photons, no isolation used
           //kPtIC: IC with cut on pT
           //kSumPtIC: IC with cut on pT sum in cone
           //kSeveralIC: Not allowed in kCorrelation analysis
 //     //for option kSeveralIC:
-//     gd->SetNCones(3); gd->SetNPtThresholds(3);
-//     gd->SetConeSizes(0,0.2); gd->SetConeSizes(1,0.3); gd->SetConeSizes(2,0.4); 
-//     gd->SetPtThresholds(0,0.5); gd->SetPtThresholds(1,1.); gd->SetPtThresholds(2,2);
+//      gd->SetNCones(2); gd->SetNPtThresholds(3);
+//      gd->SetConeSizes(0,0.2); gd->SetConeSizes(1,0.3); //gd->SetConeSizes(2,0.4); 
+//      gd->SetPtThresholds(0,0); gd->SetPtThresholds(1,1.); gd->SetPtThresholds(2,2);
  
     //============================
     //Second, select the correlation algoritm
@@ -86,13 +93,13 @@ AliAnaGamma*  ConfigGammaAnalysis()
     //No associated setters and getters for the moment.
 
     //--- Option 2 ---
-    //       AliAnaGammaHadron *gc = new AliAnaGammaHadron();
+    AliAnaGammaHadron *gc = new AliAnaGammaHadron();
 //     gc->SetDeltaPhiCutRange(1,4); //Correlation in delta phi (gamma-particle), radians
 //     gc->SetMinPtHadron(5.);
 //     gc->SetJetsOnlyInCTS(kFALSE); // Don't consider particles in opposite calorimeter
     
     //--- Option 3 ---
-    AliAnaGammaJetLeadCone *gc = new AliAnaGammaJetLeadCone();
+    //    AliAnaGammaJetLeadCone *gc = new AliAnaGammaJetLeadCone();
 //     gc->SetDeltaPhiCutRange(2.8,3.5); //Correlation with leading particle delta phi (gamma-particle), radians
 //     gc->SetRatioCutRange(0.1,1.2);
 //     gc->SetJetsOnlyInCTS(kFALSE); // Don't consider particles in opposite calorimeter
@@ -106,7 +113,7 @@ AliAnaGamma*  ConfigGammaAnalysis()
     //correlated with the opposite prompt photon. Need to play with this class
     AliNeutralMesonSelection *nms = new AliNeutralMesonSelection();
     nms->SetInvMassCutRange(0.1,0.17);
-    nms->KeepNeutralMesonSelectionHistos(kFALSE); //keep in file several histograms or not.
+    nms->KeepNeutralMesonSelectionHistos(kTRUE); //keep in file several histograms or not.
     // and other parameters
 
     //---------------------------------------------------------------------
@@ -114,7 +121,7 @@ AliAnaGamma*  ConfigGammaAnalysis()
     //---------------------------------------------------------------------
     ana = new AliAnaGamma();
     ana->SetReader(reader);//pointer to reader
-    ana->SetAnalysisType(AliAnaGamma::kPrompt); //set kPrompt, kCorrelation
+    ana->SetAnalysisType(AliAnaGamma::kCorrelation); //set kPrompt, kCorrelation
     ana->SetGammaDirect(gd);//pointer to direct photon algorithm
     ana->SetGammaCorrelation(gc);//pointer to correlation algorithm
     ana->SetNeutralMesonSelection(nms); //pi0 pair selection
index bbaf6c6..a02c647 100644 (file)
@@ -22,9 +22,9 @@ Bool_t LoadLib( const char* pararchivename)
     gROOT->ProcessLine(processline) ;
     gSystem->ChangeDirectory(cdir) ; 
     sprintf(processline, ".! mv /tmp/%s .", parpar) ;
-    gROOT->ProcessLine(processline) ;  
+    //gROOT->ProcessLine(processline) ;        
     sprintf(processline,".! tar xvzf %s",parpar);
-    gROOT->ProcessLine(processline);
+    //gROOT->ProcessLine(processline);
    }
    gSystem->ChangeDirectory(pararchivename);
    
@@ -62,11 +62,10 @@ Bool_t LoadLib( const char* pararchivename)
 void ana(const Int_t kEvent=10)  
 { 
 
-  //AliLog::SetGlobalDebugLevel(5);
+  //AliLog::SetGlobalDebugLevel(1);
   gSystem->AddIncludePath("-I\"$ALICE_ROOT/include\"");
   gSystem->Load("libANALYSIS.so");
-  gSystem->Load("libANALYSISRL.so");
-
   if (! gIsAnalysisLoaded ) {
     //    LoadLib("ESD") ; 
     //   LoadLib("AOD") ;
@@ -77,10 +76,10 @@ void ana(const Int_t kEvent=10)
   // create the analysis goodies object
   AliAnalysisGoodies * ag = new AliAnalysisGoodies() ; 
   
-  AliAnalysisTaskGamma * task = new AliAnalysisTaskGamma ("GammaCorrelations");
+  AliAnalysisTaskGamma * task = new AliAnalysisTaskGamma ("Gamma");
   ag->ConnectInput(task, TChain::Class(), 0) ; 
-  ag->ConnectOuput(task, TList::Class(), 0) ;  
-  //ag->ConnectOuput(task, TTree::Class(), 0, "AOD") ;  
+  ag->ConnectOuput(task, TTree::Class(), 0, "AOD") ;  
+  ag->ConnectOuput(task, TList::Class(), 1) ;  
    
   // get the data to analyze