7be13e68 |
1 | #include "AliHBTReaderESD.h" |
7e35966a |
2 | //____________________________________________________________________ |
3 | ////////////////////////////////////////////////////////////////////// |
4 | // // |
5 | // class AliHBTReaderESD // |
6 | // // |
7 | // reader for ALICE Event Summary Data (ESD). // |
8 | // // |
9 | // Piotr.Skowronski@cern.ch // |
10 | // // |
11 | ////////////////////////////////////////////////////////////////////// |
7be13e68 |
12 | |
13 | #include <TPDGCode.h> |
14 | #include <TString.h> |
15 | #include <TObjString.h> |
16 | #include <TTree.h> |
17 | #include <TFile.h> |
18 | #include <TParticle.h> |
19 | |
bed069a4 |
20 | #include <AliRunLoader.h> |
21 | #include <AliStack.h> |
7be13e68 |
22 | #include <AliESDtrack.h> |
23 | #include <AliESD.h> |
24 | |
25 | #include "AliHBTRun.h" |
26 | #include "AliHBTEvent.h" |
27 | #include "AliHBTParticle.h" |
28 | #include "AliHBTParticleCut.h" |
29 | |
30 | ClassImp(AliHBTReaderESD) |
31 | |
bed069a4 |
32 | AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename): |
7be13e68 |
33 | fESDFileName(esdfilename), |
bed069a4 |
34 | fGAlFileName(galfilename), |
35 | fFile(0x0), |
36 | fRunLoader(0x0), |
37 | fReadParticles(kFALSE) |
7be13e68 |
38 | { |
39 | //cosntructor |
40 | if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES)) |
41 | Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra."); |
42 | } |
43 | /********************************************************************/ |
44 | |
bed069a4 |
45 | AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename): |
7be13e68 |
46 | AliHBTReader(dirs), |
7be13e68 |
47 | fESDFileName(esdfilename), |
bed069a4 |
48 | fGAlFileName(galfilename), |
49 | fFile(0x0), |
50 | fRunLoader(0x0), |
51 | fReadParticles(kFALSE) |
7be13e68 |
52 | { |
53 | //cosntructor |
54 | if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES)) |
55 | Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra."); |
56 | } |
57 | /********************************************************************/ |
58 | |
59 | AliHBTReaderESD::~AliHBTReaderESD() |
60 | { |
61 | //desctructor |
bed069a4 |
62 | delete fRunLoader; |
7e35966a |
63 | delete fFile; |
7be13e68 |
64 | } |
bed069a4 |
65 | /**********************************************************/ |
66 | Int_t AliHBTReaderESD::ReadNext() |
7be13e68 |
67 | { |
bed069a4 |
68 | //reads next event from fFile |
69 | |
70 | //****** Tentative particle type "concentrations" |
71 | static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05}; |
7be13e68 |
72 | |
73 | Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track |
74 | Double_t w[kNSpecies]; |
75 | Double_t mom[3];//momentum |
76 | Double_t pos[3];//position |
7be13e68 |
77 | |
83a9f895 |
78 | if (AliHBTParticle::GetDebug()) |
bed069a4 |
79 | Info("ReadNext","Entered"); |
80 | |
7be13e68 |
81 | TDatabasePDG* pdgdb = TDatabasePDG::Instance(); |
82 | if (pdgdb == 0x0) |
83 | { |
bed069a4 |
84 | Error("Next","Can not get PDG Database Instance."); |
7be13e68 |
85 | return 1; |
86 | } |
bed069a4 |
87 | |
88 | if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent(); |
89 | if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent(); |
90 | |
91 | fParticlesEvent->Reset(); |
92 | fTracksEvent->Reset(); |
93 | |
7be13e68 |
94 | do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./" |
95 | { |
bed069a4 |
96 | if (fFile == 0x0) |
7be13e68 |
97 | { |
bed069a4 |
98 | fFile = OpenFile(fCurrentDir);//rl is opened here |
99 | if (fFile == 0x0) |
7be13e68 |
100 | { |
bed069a4 |
101 | Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir); |
102 | fCurrentDir++; |
103 | continue; |
104 | } |
105 | fCurrentEvent = 0; |
106 | } |
107 | |
108 | //try to read |
109 | TString esdname; |
110 | esdname+=fCurrentEvent; |
111 | AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname)); |
112 | if (esd == 0x0) |
113 | { |
83a9f895 |
114 | if (AliHBTParticle::GetDebug() > 2 ) |
bed069a4 |
115 | { |
116 | Info("ReadNext","Can not find ESD object named %s.",esdname.Data()); |
117 | } |
118 | fCurrentDir++; |
119 | delete fFile;//we have to assume there is no more ESD objects in the fFile |
120 | delete fRunLoader; |
121 | fFile = 0x0; |
7e35966a |
122 | fRunLoader = 0x0; |
bed069a4 |
123 | continue; |
124 | } |
125 | AliStack* stack = 0x0; |
126 | if (fReadParticles && fRunLoader) |
127 | { |
128 | fRunLoader->GetEvent(fCurrentEvent); |
129 | stack = fRunLoader->Stack(); |
130 | } |
131 | Info("ReadNext","Reading Event %d",fCurrentEvent); |
7be13e68 |
132 | |
bed069a4 |
133 | Int_t ntr = esd->GetNumberOfTracks(); |
134 | Info("ReadNext","Found %d tracks.",ntr); |
135 | for (Int_t i = 0;i<ntr; i++) |
136 | { |
137 | AliESDtrack *esdtrack = esd->GetTrack(i); |
138 | if (esdtrack == 0x0) |
139 | { |
140 | Error("Next","Can not get track %d", i); |
141 | continue; |
142 | } |
7e35966a |
143 | |
aafd28f6 |
144 | //if (esdtrack->HasVertexParameters() == kFALSE) |
145 | if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE) |
7e35966a |
146 | { |
83a9f895 |
147 | if (AliHBTParticle::GetDebug() > 2) |
7e35966a |
148 | Info("ReadNext","Particle skipped: Data at vertex not available."); |
149 | continue; |
150 | } |
151 | |
bed069a4 |
152 | if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) |
153 | { |
83a9f895 |
154 | if (AliHBTParticle::GetDebug() > 2) |
bed069a4 |
155 | Info("ReadNext","Particle skipped: PID BIT is not set."); |
156 | continue; |
157 | } |
7e35966a |
158 | |
bed069a4 |
159 | esdtrack->GetESDpid(pidtable); |
aafd28f6 |
160 | //esdtrack->GetVertexPxPyPz(mom[0],mom[1],mom[2]); |
161 | esdtrack->GetPxPyPz(mom); |
162 | //esdtrack->GetVertexXYZ(pos[0],pos[1],pos[2]); |
163 | esdtrack->GetXYZ(pos); |
bed069a4 |
164 | |
3c69895b |
165 | Double_t extx; |
166 | Double_t extp[5]; |
167 | esdtrack->GetExternalParameters(extx,extp); |
168 | |
169 | Int_t charge = (extp[4] > 0)?-1:1;//if curvature=-charg/Pt is positive charge is negative |
170 | |
bed069a4 |
171 | //Particle from kinematics |
172 | AliHBTParticle* particle = 0; |
173 | Bool_t keeppart = kFALSE; |
174 | if ( fReadParticles && stack ) |
175 | { |
176 | if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track |
177 | TParticle *p = stack->Particle(esdtrack->GetLabel()); |
178 | if (p==0x0) |
179 | { |
180 | Error("ReadNext","Can not find track with such label."); |
181 | continue; |
182 | } |
3c69895b |
183 | // if(p->GetPdgCode()<0) charge = -1; |
bed069a4 |
184 | particle = new AliHBTParticle(*p,i); |
3c69895b |
185 | |
bed069a4 |
186 | } |
187 | |
188 | //Here we apply Bayes' formula |
189 | Double_t rc=0.; |
190 | for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s]; |
191 | if (rc==0.0) |
192 | { |
83a9f895 |
193 | if (AliHBTParticle::GetDebug() > 2) |
bed069a4 |
194 | Info("ReadNext","Particle rejected since total bayessian PID probab. is zero."); |
195 | continue; |
196 | } |
7be13e68 |
197 | |
bed069a4 |
198 | for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc; |
7be13e68 |
199 | |
83a9f895 |
200 | if (AliHBTParticle::GetDebug() > 4) |
bed069a4 |
201 | { |
202 | Info("ReadNext","###########################################################################"); |
203 | Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]); |
204 | Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]); |
205 | TString msg("Pid list got from track:"); |
206 | for (Int_t s = 0;s<kNSpecies;s++) |
207 | { |
208 | msg+="\n "; |
209 | msg+=s; |
210 | msg+="("; |
3c69895b |
211 | msg+=charge*GetSpeciesPdgCode((ESpecies)s); |
bed069a4 |
212 | msg+="): "; |
213 | msg+=w[s]; |
214 | msg+=" ("; |
215 | msg+=pidtable[s]; |
216 | msg+=")"; |
217 | } |
218 | Info("ReadNext","%s",msg.Data()); |
83a9f895 |
219 | }//if (AliHBTParticle::GetDebug()>4) |
7be13e68 |
220 | |
bed069a4 |
221 | for (Int_t s = 0; s<kNSpecies; s++) |
222 | { |
dd82cadc |
223 | Float_t pp = w[s]; |
224 | if (pp == 0.0) continue; |
7be13e68 |
225 | |
3c69895b |
226 | Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s); |
bed069a4 |
227 | if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type |
7be13e68 |
228 | |
bed069a4 |
229 | Double_t mass = pdgdb->GetParticle(pdgcode)->Mass(); |
230 | Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track |
7be13e68 |
231 | |
bed069a4 |
232 | AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, |
233 | mom[0], mom[1], mom[2], tEtot, |
234 | pos[0], pos[1], pos[2], 0.); |
235 | //copy probabilitis of other species (if not zero) |
236 | for (Int_t k = 0; k<kNSpecies; k++) |
237 | { |
238 | if (k == s) continue; |
239 | if (w[k] == 0.0) continue; |
3c69895b |
240 | track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]); |
bed069a4 |
241 | } |
242 | if(Pass(track))//check if meets all criteria of any of our cuts |
243 | //if it does not delete it and take next good track |
244 | { |
245 | delete track; |
246 | continue; |
247 | } |
248 | |
249 | fTracksEvent->AddParticle(track); |
250 | if (particle) fParticlesEvent->AddParticle(particle); |
251 | keeppart = kTRUE; |
252 | |
83a9f895 |
253 | if (AliHBTParticle::GetDebug() > 4 ) |
bed069a4 |
254 | { |
255 | Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode); |
256 | track->Print(); |
257 | } |
258 | }//for (Int_t s = 0; s<kNSpecies; s++) |
259 | |
260 | if (keeppart == kFALSE) delete particle;//particle was not stored in event |
261 | |
262 | }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks |
263 | |
264 | Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).", |
dd82cadc |
265 | fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(), |
bed069a4 |
266 | fNEventsRead,fCurrentEvent,fCurrentDir); |
7be13e68 |
267 | |
bed069a4 |
268 | fCurrentEvent++; |
269 | fNEventsRead++; |
270 | delete esd; |
271 | return 0;//success -> read one event |
272 | }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array |
273 | |
274 | return 1; //no more directories to read |
275 | } |
276 | /**********************************************************/ |
277 | void AliHBTReaderESD::Rewind() |
278 | { |
7e35966a |
279 | //rewinds reading |
bed069a4 |
280 | delete fFile; |
281 | fFile = 0; |
282 | delete fRunLoader; |
283 | fCurrentDir = 0; |
284 | fNEventsRead = 0; |
285 | fCurrentEvent++; |
bed069a4 |
286 | } |
7be13e68 |
287 | /**********************************************************/ |
bed069a4 |
288 | |
7be13e68 |
289 | TFile* AliHBTReaderESD::OpenFile(Int_t n) |
290 | { |
bed069a4 |
291 | //opens fFile with kine tree |
7be13e68 |
292 | |
293 | const TString& dirname = GetDirName(n); |
294 | if (dirname == "") |
295 | { |
296 | Error("OpenFiles","Can not get directory name"); |
297 | return 0x0; |
298 | } |
299 | TString filename = dirname +"/"+ fESDFileName; |
300 | TFile *ret = TFile::Open(filename.Data()); |
301 | |
302 | if ( ret == 0x0) |
303 | { |
bed069a4 |
304 | Error("OpenFiles","Can't open fFile %s",filename.Data()); |
7be13e68 |
305 | return 0x0; |
306 | } |
307 | if (!ret->IsOpen()) |
308 | { |
bed069a4 |
309 | Error("OpenFiles","Can't open fFile %s",filename.Data()); |
7be13e68 |
310 | return 0x0; |
311 | } |
312 | |
bed069a4 |
313 | if (fReadParticles) |
314 | { |
315 | fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName); |
316 | if (fRunLoader == 0x0) |
317 | { |
318 | Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data()); |
319 | delete ret; |
320 | return 0x0; |
321 | } |
322 | fRunLoader->LoadHeader(); |
323 | if (fRunLoader->LoadKinematics()) |
324 | { |
325 | Error("Next","Error occured while loading kinematics."); |
326 | delete fRunLoader; |
327 | delete ret; |
328 | return 0x0; |
329 | } |
330 | } |
7be13e68 |
331 | return ret; |
332 | } |
333 | /**********************************************************/ |
334 | |
335 | Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron |
336 | { |
337 | //returns pdg code from the PID index |
338 | //ask jura about charge |
339 | switch (spec) |
340 | { |
341 | case kESDElectron: |
342 | return kPositron; |
343 | break; |
344 | case kESDMuon: |
345 | return kMuonPlus; |
346 | break; |
347 | case kESDPion: |
348 | return kPiPlus; |
349 | break; |
350 | case kESDKaon: |
351 | return kKPlus; |
352 | break; |
353 | case kESDProton: |
354 | return kProton; |
355 | break; |
356 | default: |
357 | ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec); |
358 | break; |
359 | } |
360 | return 0; |
361 | } |