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