correct previous non intended commit, but leaving some new things
[u/mrichter/AliRoot.git] / EMCAL / macros / TestESDCaloClusterAndCell.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2
3 //Root include files
4 #include <Riostream.h>
5 #include <TFile.h>
6 #include <TChain.h>
7 #include <TParticle.h>
8 #include <TNtuple.h>
9 #include <TCanvas.h>
10 #include <TObjArray.h>
11 #include <TSystem.h>
12 #include <TString.h> 
13 #include <TH1F.h> 
14 #include <TVector.h> 
15 #include <TParticle.h> 
16 #include <TRefArray.h>
17 #include <TArrayS.h>
18
19 //AliRoot include files
20 #include "AliRunLoader.h"
21 #include "AliStack.h"
22 #include "AliESDEvent.h"
23 #include "AliESDVertex.h"
24 #include "AliESDCaloCluster.h"
25 #include "AliESDCaloCells.h"
26 #include "AliPID.h"
27 #include "AliLog.h" 
28
29 #endif
30
31 //example macro to extract information from AliESDCaloClusters and AliESDCaloCells.
32 // Author: Gustavo.Conesa.Balbastre@cern.ch (INFN-LNF)
33
34 TChain * AliReadESDfromdisk(const UInt_t eventsToRead, 
35                                    const TString dirName, 
36                                    const TString esdTreeName, 
37                                    const char *  pattern) 
38 {
39   // Reads ESDs from Disk
40  TChain *  rv = 0; 
41   
42   // create a TChain of all the files 
43   TChain * cESDTree = new TChain(esdTreeName) ; 
44
45   // read from the directory file until the require number of events are collected
46   void * from = gSystem->OpenDirectory(dirName) ;
47    if (!from) 
48      rv = 0 ;   
49   else{ // reading file names from directory
50     const char * subdir ; 
51     // search all subdirectories witch matching pattern
52     while( (subdir = gSystem->GetDirEntry(from))  && 
53            (cESDTree->GetEntries() < eventsToRead)) {
54       if ( strstr(subdir, pattern) != 0 ) { 
55         char file[200] ; 
56         sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir);     
57         cESDTree->Add(file) ;
58       }
59     } // while file
60   
61     rv = cESDTree ; 
62     
63   } // reading file names from directory
64   return rv ; 
65 }
66
67 //======================================================================
68 TChain * AliReadESD(const UInt_t eventsToRead,
69                   const TString dirName, 
70                   const TString esdTreeName, 
71                   const char *  pattern)  
72 {
73   // Read AliESDs files and return a Chain of events
74  
75   if ( dirName == "" ) 
76     return 0 ; 
77   if ( esdTreeName == "" ) 
78     return AliReadESDfromdisk(eventsToRead, dirName,"","") ;//Last 2 arguments are not necessary but pdsf compiler complains "","'
79   else if ( strcmp(pattern, "") == 0 )
80     return AliReadESDfromdisk(eventsToRead, dirName, esdTreeName,"") ;//Last argument is not necessary but pdsf compiler complains "","'
81   else 
82     return AliReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ;        
83 }
84
85 //=====================================================================
86 //  Do:
87 //  .L TestESDCaloClusterAndCell.C++
88 //  TestESDCaloClusterAndCell(calo, number of events to process)
89 //=====================================================================
90 void TestESDCaloClusterAndCell(TString det = "EMCAL", const UInt_t eventsToProcess = 5, 
91                         TString dirName = "./", 
92                         const TString esdTreeName = "esdTree", 
93                         const char *  pattern = ".") 
94
95
96   if(det !="EMCAL" && det !="PHOS" )
97  {
98    cout << "Wrong detector name "<<det<<endl; 
99    return;
100  }
101
102   //Create chain of esd trees
103   //By default the root files are in the same directory 
104   TChain * t = AliReadESD(eventsToProcess, dirName,esdTreeName,pattern) ; 
105  
106   
107   // ESD
108   AliESDEvent * esd = new AliESDEvent(); 
109   esd->ReadFromTree(t); // This checks also if the branch address is already set
110
111   
112   //Define few variables to be used in macro
113   TString alirunName = "" ; 
114
115   UInt_t event ;
116   Float_t pos[3] ; 
117
118   for (event = 0; event < eventsToProcess; event++) {//event loop
119     //AliInfo( Form("Event %d \n",event) );  
120     Int_t nbytes = t->GetEntry(event); // store event in esd
121     if ( nbytes == 0 ) //If nothing in ESD 
122       break ; 
123     
124     // Check that name of file is correct
125     if (alirunName != t->GetFile()->GetName()) {        
126       alirunName = t->GetFile()->GetName() ; 
127       alirunName.ReplaceAll("galice.root", "AliESDs.root") ;
128     }
129     
130     AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),  "read");
131     rl->LoadKinematics();  
132     rl->GetEvent(event);
133     AliStack *sta=rl->Stack();
134
135     //get reconstructed vertex position     
136     Double_t vertex_position[3] ; 
137     esd->GetVertex()->GetXYZ(vertex_position) ; 
138     
139     //cout<<"Event >>>>>>>>>>>  "<<event<<" vertex >>> "<<vertex_position[0]<<" "<<vertex_position[1]<<" "<<vertex_position[2]<<endl;
140     
141     //Get CaloCells
142
143     AliESDCaloCells &cells= *(esd->GetEMCALCells());
144     if(det == "PHOS")  AliESDCaloCells &cells = *(esd->GetPHOSCells());
145     Int_t ncell = cells.GetNumberOfCells() ;  
146     Int_t type = cells.GetType();
147     cout<<">>> Event "<<event<<" ncells "<<ncell<< " type "<<type<<endl;
148     // Uncomment to see the full list of digit amplitudes and times.
149 //     for (Int_t icell=  0; icell <  ncell; icell++) {
150 //       cout<<"icell "<<icell<<
151 //      " ID "<<cells.GetCellNumber(icell)<<
152 //      " amp "<<cells.GetAmplitude(icell)<<
153 //      " time "<<cells.GetTime(icell)<<endl;
154
155 //     }// cell loop
156
157     //Get the CaloClusters
158
159     //select EMCAL clusters only 
160     TRefArray * caloClusters  = new TRefArray();
161
162     if(det == "EMCAL") esd->GetEMCALClusters(caloClusters);
163     else if(det == "PHOS") esd->GetPHOSClusters(caloClusters);
164     
165     Int_t nclus = caloClusters->GetEntries() ;  
166    
167     cout<<"Event: "<<event<<"; Number of clusters "<<nclus<<endl;   
168     for (Int_t iclus =  0; iclus <  nclus; iclus++) {//////////////EMCAL cluster loop
169       AliESDCaloCluster * clus = (AliESDCaloCluster *) caloClusters->At(iclus) ; // retrieve cluster from esd
170       //cout<<">>>> Cluster >>>"<<iclus<<endl;
171         
172       //Get the cluster info
173       
174       Float_t energy   = clus->E() ;
175       //Float_t disp     = clus->GetClusterDisp() ;
176       Int_t iprim = clus->GetLabel();
177       
178       clus->GetPosition(pos) ;
179       TVector3 vpos(pos[0],pos[1],pos[2]) ;
180       TLorentzVector p ;
181       clus->GetMomentum(p,vertex_position);
182       
183       Int_t mult       = clus->GetNCells() ;
184       UShort_t * index = clus->GetCellsAbsId() ;
185       Double_t * fraction = clus->GetCellsAmplitudeFraction() ;
186       //Print cluster info
187       cout<<"Cluster "<<iclus <<"; digits mult "<<mult<<"; type "<<(Int_t )clus->GetClusterType()
188           <<"; Energy "<<energy
189           <<"; Phi "<<vpos.Phi()*180/TMath::Pi()<<"; Eta "<<vpos.Eta()
190           <<"; label "<<iprim<<endl;      
191       //Print primary info
192       if(iprim>=0 && iprim < sta->GetNtrack()){
193         TParticle * particle = sta->Particle(clus->GetLabel());
194         //Print primary values
195         cout<<" Primary: "<<particle->GetName()<< "; Energy "<<particle->Energy()<<endl;    
196       }
197       else if( iprim >= sta->GetNtrack()) cout <<"PROBLEM, label is too large : "<<iprim<<" >= particles in stack "<< sta->GetNtrack() <<endl;
198       else cout<<"Negative label!!!  : "<<iprim<<endl;
199       
200       //Get CaloCells of cluster and print       
201       for(Int_t i = 0; i < mult ; i++){
202         Int_t absId =   index[i]; // or clus->GetCellNumber(i) ;
203         Double_t ampFract =  fraction[i];
204         Float_t amp = cells.GetCellAmplitude(absId) ;
205         Float_t time = cells.GetCellTime(absId);
206         cout<<"AbsID : "<< absId << "; Amplitude "<< amp << "; Fraction "<<ampFract<<"; Time " <<time<<endl;
207       }
208       
209     }// clusters
210     
211   }// event loop
212   
213 }
214
215