]>
Commit | Line | Data |
---|---|---|
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 | ||
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 |