suggestion by Chiara to add TestPreprocessor macro to svn
[u/mrichter/AliRoot.git] / EMCAL / macros / TestESDCaloClusterAndCell.C
CommitLineData
34f81413 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
34TChain * 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//======================================================================
68TChain * 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//=====================================================================
90void 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