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