Corrected event counter.
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderESD.cxx
CommitLineData
d7c6ab14 1// $Id$
2
3//____________________________________________________________________
4//////////////////////////////////////////////////////////////////////
5// //
6// class AliHBTReaderESD //
7// //
8// Reader for ALICE Event Summary Data (ESD). //
9// //
10// Piotr.Skowronski@cern.ch //
11// //
12//////////////////////////////////////////////////////////////////////
13
d9a296b1 14#include <Riostream.h>
d7c6ab14 15#include <TMath.h>
16#include <TPDGCode.h>
17#include <TString.h>
18#include <TObjString.h>
19#include <TTree.h>
20#include <TFile.h>
21#include <TKey.h>
22#include <TH1.h>
23#include <AliESD.h>
24#include <AliESDtrack.h>
04a02430 25#include <AliKalmanTrack.h>
d7c6ab14 26#include <AliJetEventParticles.h>
27#include "AliJetParticlesReaderESD.h"
28
29ClassImp(AliJetParticlesReaderESD)
30
301a24f1 31AliJetParticlesReaderESD::AliJetParticlesReaderESD(Bool_t constrained,
32 const Char_t* esdfilename) :
d7c6ab14 33 AliJetParticlesReader(),
301a24f1 34 fConstrained(constrained),
d7c6ab14 35 fESDFileName(esdfilename),
5d60c8f9 36 fESD(0),
d7c6ab14 37 fFile(0),
5d60c8f9 38 fTree(0),
d7c6ab14 39 fKeyIterator(0),
40 fPassFlag(AliESDtrack::kTPCrefit)
41{
42 //constructor
43}
44
45/********************************************************************/
46
47AliJetParticlesReaderESD::AliJetParticlesReaderESD(
301a24f1 48 Bool_t constrained,
d7c6ab14 49 TObjArray* dirs,
50 const Char_t* esdfilename) :
51 AliJetParticlesReader(dirs),
301a24f1 52 fConstrained(constrained),
d7c6ab14 53 fESDFileName(esdfilename),
5d60c8f9 54 fESD(0),
d7c6ab14 55 fFile(0),
5d60c8f9 56 fTree(0),
d7c6ab14 57 fKeyIterator(0),
58 fPassFlag(AliESDtrack::kTPCrefit)
59{
60 //constructor
61}
62
63/********************************************************************/
64
65AliJetParticlesReaderESD::~AliJetParticlesReaderESD()
66{
67 //desctructor
5d60c8f9 68 if(fESD) delete fESD;
69 if(fTree) delete fTree;
d7c6ab14 70 if(fKeyIterator) delete fKeyIterator;
04a02430 71 if(fFile) delete fFile;
d7c6ab14 72}
73
d7c6ab14 74/**********************************************************/
75
76Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
77{
78 //Reads one ESD
79
80 if (esd == 0)
81 {
82 Error("ReadESD","ESD is NULL");
83 return kFALSE;
84 }
85
86 //TDatabasePDG* pdgdb = TDatabasePDG::Instance();
87 //if (pdgdb == 0)
88 //{
89 // Error("ReadESD","Can not get PDG Database Instance.");
90 // return kFALSE;
91 // }
92
04a02430 93 Float_t mf = esd->GetMagneticField();
94 if (mf <= 0.0)
95 {
96 Error("ReadESD","Magnetic Field is 0. Skipping to next event.");
97 return kFALSE;
98 }
99 AliKalmanTrack::SetMagneticField(mf/10.);
d7c6ab14 100
a86edc4e 101 Info("ReadESD","Reading Event %d",fCurrentDir*1000+fCurrentEvent);
d7c6ab14 102 if((!fOwner) || (fEventParticles==0))
103 fEventParticles = new AliJetEventParticles();
104
04a02430 105 const Int_t kntr = esd->GetNumberOfTracks();
106 Info("ReadESD","Found %d tracks.",kntr);
107 fEventParticles->Reset(kntr);
108
109 TString headdesc="";
110 headdesc+="Run ";
111 headdesc+=esd->GetRunNumber();
112 headdesc+=": Ev ";
113 headdesc+=esd->GetEventNumber();
114 fEventParticles->SetHeader(headdesc);
115
d7c6ab14 116 Double_t vertexpos[3];//vertex position
117 const AliESDVertex* kvertex = esd->GetVertex();
118 if (kvertex == 0)
119 {
120 Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
121 vertexpos[0] = 0.0;
122 vertexpos[1] = 0.0;
123 vertexpos[2] = 0.0;
124 }
125 else
126 {
127 kvertex->GetXYZ(vertexpos);
128 }
129 fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
130
04a02430 131 //loop over tracks
d7c6ab14 132 for (Int_t i = 0;i<kntr; i++)
133 {
bad753b1 134
d7c6ab14 135 const AliESDtrack *kesdtrack = esd->GetTrack(i);
136 if (kesdtrack == 0)
137 {
04a02430 138 Error("ReadESD","Can't get track %d", i);
d7c6ab14 139 continue;
140 }
45cf33d2 141
142 ULong_t cmptest=(kesdtrack->GetStatus() & fPassFlag);
143 if (cmptest!=fPassFlag)
d7c6ab14 144 {
70d6bc80 145 //cout << i << " "; PrintESDtrack(kesdtrack); cout << endl;
45cf33d2 146 Info("ReadNext","Particle %d skipped: %u.",i,kesdtrack->GetStatus());
d7c6ab14 147 continue;
148 }
149
150 Double_t mom[3]; //momentum
301a24f1 151 Double_t xyz[3]; //position
152 if (fConstrained) {
153 if (kesdtrack->GetConstrainedChi2() > 25) continue;
154 kesdtrack->GetConstrainedPxPyPz(mom);
155 kesdtrack->GetConstrainedXYZ(xyz);
156 } else {
157 kesdtrack->GetPxPyPz(mom);
158 kesdtrack->GetXYZ(xyz);
159 }
d7c6ab14 160 const Float_t kmass=kesdtrack->GetMass();
161 const Float_t kp2=mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2];
162 const Float_t ketot=TMath::Sqrt(kmass*kmass+kp2);
163 const Float_t kpt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
164 const Float_t kp=TMath::Sqrt(kp2);
165 const Float_t keta=0.5*TMath::Log((kp+mom[2]+1e-30)/(kp-mom[2]+1e-30));
166 const Float_t kphi=TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
301a24f1 167 //Double_t dx = xyz[0]-vertexpos[0];
168 //Double_t dy = xyz[1]-vertexpos[1];
169 //Float_t dca = TMath::Sqrt(dx*dx + dy*dy);
170 //Float_t dz = xyz[2]-vertexpos[2];
171 UInt_t index[6];
172 const Int_t kncl=kesdtrack->GetITSclusters(index)
173 +kesdtrack->GetTPCclusters(NULL)
174 +kesdtrack->GetTRDclusters(NULL);
d7c6ab14 175 if(IsAcceptedParticle(kpt,kphi,keta))
301a24f1 176 fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kncl,kpt,kphi,keta);
d7c6ab14 177
178 } // loop over tracks
179
180 return kTRUE;
181}
182
183/**********************************************************/
184
185void AliJetParticlesReaderESD::Rewind()
186{
187 //rewinds reading
188
189 if(fFile) delete fFile;
190 if(fKeyIterator) delete fKeyIterator;
191 fKeyIterator = 0;
192 fFile = 0;
193 fCurrentDir = 0;
194 fNEventsRead = 0;
195}
196
197/**********************************************************/
198
d9a296b1 199Int_t AliJetParticlesReaderESD::ReadNext()
200{
201 //reads next event from fFile
202
203 do // is OK even if 0 dirs specified,
204 { // in that case we try to read from "./"
205 if (fFile == 0)
206 {
207 fFile = OpenFile(fCurrentDir);
208 if (fFile == 0)
209 {
210 Error("ReadNext","Can't get fFile for dir no. %d",fCurrentDir);
211 fCurrentDir++;
212 continue;
213 }
214
215 fCurrentEvent = 0;
216 //fFile->GetListOfKeys()->Print();
217
218 if(fTree) delete fTree;
219 fTree = dynamic_cast<TTree*>(fFile->Get("esdTree"));
220 if(fTree)
221 fTree->SetBranchAddress("ESD",&fESD);
222 else
223 fKeyIterator = new TIter(fFile->GetListOfKeys());
224 }
225
226 if(fTree)
227 {
228 if(AliKalmanTrack::GetConvConst()<=0.)
229 AliKalmanTrack::SetMagneticField(0.5);
230 if(fCurrentEvent>=fTree->GetEntries())
231 {
232 fCurrentDir++;
233 delete fTree;
234 fTree = 0;
235 delete fFile;
236 fFile = 0;
237 continue;
238 }
239 fTree->GetEvent(fCurrentEvent);
240 }
241 else
242 { // "old" way via ESD objects stored in root file
243 TKey* key = (TKey*)fKeyIterator->Next();
244 if (key == 0)
245 {
246 fCurrentDir++;
247 delete fKeyIterator;
248 fKeyIterator = 0;
249 delete fFile; //we have to assume there are no more ESD objects in the fFile
250 fFile = 0;
251 continue;
252 }
253 TString esdname = "ESD";
254 esdname+=fCurrentEvent;
255 if(fESD) delete fESD;
256 fESD = dynamic_cast<AliESD*>(fFile->Get(esdname));
257 if (fESD == 0)
258 {
259 Info("ReadNext","Can't find AliESD object named %s",esdname.Data());
260 fCurrentDir++;
261 delete fKeyIterator;
262 fKeyIterator = 0;
263 delete fFile;//we have to assume there is no more ESD objects in the fFile
264 fFile = 0;
265 continue;
266 }
267 }
268 ReadESD(fESD);
269 fCurrentEvent++;
270 fNEventsRead++;
271 return kTRUE;//success -> read one event
272 } while(fCurrentDir < GetNumberOfDirs());
273 //end of loop over directories specified in fDirs Obj Array
274 return kFALSE; //no more directories to read
275}
276
277/**********************************************************/
278
d7c6ab14 279TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
280{
281 //opens fFile with kine tree
282
283 const TString& dirname = GetDirName(n);
284 if (dirname == "")
285 {
04a02430 286 Error("OpenFiles","Can't get directory name");
d7c6ab14 287 return 0;
288 }
289 TString filename = dirname +"/"+ fESDFileName;
290 TFile *ret = TFile::Open(filename.Data());
291 if (ret == 0)
292 {
293 Error("OpenFiles","Can't open fFile %s",filename.Data());
294 return 0;
295 }
296 if (!ret->IsOpen())
297 {
298 Error("OpenFiles","Can't open fFile %s",filename.Data());
299 return 0;
300 }
301
302 return ret;
303}
45cf33d2 304
305/**********************************************************/
306
307void AliJetParticlesReaderESD::PrintESDtrack(const AliESDtrack *kesdtrack) const
308{
309 ULong_t status=kesdtrack->GetStatus();
310 cout << hex << status << dec << ": ";
311 if((status & AliESDtrack::kITSin) == AliESDtrack::kITSin) cout << "ITSin ";
312 if((status & AliESDtrack::kITSout) == AliESDtrack::kITSout) cout << "ITSout ";
313 if((status & AliESDtrack::kITSrefit) == AliESDtrack::kITSrefit) cout << "ITSrefit ";
314 if((status & AliESDtrack::kITSpid) == AliESDtrack::kITSpid) cout << "ITSpid ";
315
316 if((status & AliESDtrack::kTPCin) == AliESDtrack::kTPCin) cout << "TPCin ";
317 if((status & AliESDtrack::kTPCout) == AliESDtrack::kTPCout) cout << "TPCout ";
318 if((status & AliESDtrack::kTPCrefit) == AliESDtrack::kTPCrefit) cout << "TPCrefit ";
319 if((status & AliESDtrack::kTPCpid) == AliESDtrack::kTPCpid) cout << "TPCpid ";
320
321 if((status & AliESDtrack::kTRDin) == AliESDtrack::kTRDin) cout << "TRDin ";
322 if((status & AliESDtrack::kTRDout) == AliESDtrack::kTRDout) cout << "TRDout ";
323 if((status & AliESDtrack::kTRDrefit) == AliESDtrack::kTRDrefit) cout << "TRDrefit ";
324 if((status & AliESDtrack::kTRDpid) == AliESDtrack::kTRDpid) cout << "TRDpid ";
325 if((status & AliESDtrack::kTRDbackup) == AliESDtrack::kTRDbackup) cout << "TRDbackup ";
326 if((status & AliESDtrack::kTRDStop) == AliESDtrack::kTRDStop) cout << "TRDstop ";
327
328 if((status & AliESDtrack::kTOFin) == AliESDtrack::kTOFin) cout << "TOFin ";
329 if((status & AliESDtrack::kTOFout) == AliESDtrack::kTOFout) cout << "TOFout ";
330 if((status & AliESDtrack::kTOFrefit) == AliESDtrack::kTOFrefit) cout << "TOFrefit ";
331 if((status & AliESDtrack::kTOFpid) == AliESDtrack::kTOFpid) cout << "TOFpid ";
332
333 if((status & AliESDtrack::kPHOSpid) == AliESDtrack::kPHOSpid) cout << "PHOSpid ";
334 if((status & AliESDtrack::kRICHpid) == AliESDtrack::kRICHpid) cout << "RICHpid ";
335 if((status & AliESDtrack::kEMCALpid) == AliESDtrack::kEMCALpid) cout << "EMCALpid ";
336 if((status & AliESDtrack::kESDpid) == AliESDtrack::kESDpid) cout << "ESDpid ";
337 if((status & AliESDtrack::kTIME) == AliESDtrack::kTIME) cout << "TIME ";
338 cout << endl;
339 /*
340 cout << kesdtrack->GetConstrainedChi2()
341 << " " << kesdtrack->GetITSchi2()
342 << " " << kesdtrack->GetTPCchi2()
343 << " " << kesdtrack->GetTRDchi2()<< endl;
344 */
345}