]>
Commit | Line | Data |
---|---|---|
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 | |
bed069a4 | 78 | if (AliHBTParticle::fgDebug) |
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 | { | |
114 | if (AliHBTParticle::fgDebug > 2 ) | |
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 | |
bfb09ece | 144 | if (esdtrack->HasVertexParameters() == kFALSE) |
7e35966a | 145 | { |
146 | if (AliHBTParticle::fgDebug > 2) | |
147 | Info("ReadNext","Particle skipped: Data at vertex not available."); | |
148 | continue; | |
149 | } | |
150 | ||
bed069a4 | 151 | if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) |
152 | { | |
153 | if (AliHBTParticle::fgDebug > 2) | |
154 | Info("ReadNext","Particle skipped: PID BIT is not set."); | |
155 | continue; | |
156 | } | |
7e35966a | 157 | |
bed069a4 | 158 | esdtrack->GetESDpid(pidtable); |
bfb09ece | 159 | esdtrack->GetVertexPxPyPz(mom[0],mom[1],mom[2]); |
160 | esdtrack->GetVertexXYZ(pos[0],pos[1],pos[2]); | |
bed069a4 | 161 | |
3c69895b | 162 | Double_t extx; |
163 | Double_t extp[5]; | |
164 | esdtrack->GetExternalParameters(extx,extp); | |
165 | ||
166 | Int_t charge = (extp[4] > 0)?-1:1;//if curvature=-charg/Pt is positive charge is negative | |
167 | ||
bed069a4 | 168 | //Particle from kinematics |
169 | AliHBTParticle* particle = 0; | |
170 | Bool_t keeppart = kFALSE; | |
171 | if ( fReadParticles && stack ) | |
172 | { | |
173 | if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track | |
174 | TParticle *p = stack->Particle(esdtrack->GetLabel()); | |
175 | if (p==0x0) | |
176 | { | |
177 | Error("ReadNext","Can not find track with such label."); | |
178 | continue; | |
179 | } | |
3c69895b | 180 | // if(p->GetPdgCode()<0) charge = -1; |
bed069a4 | 181 | particle = new AliHBTParticle(*p,i); |
3c69895b | 182 | |
bed069a4 | 183 | } |
184 | ||
185 | //Here we apply Bayes' formula | |
186 | Double_t rc=0.; | |
187 | for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s]; | |
188 | if (rc==0.0) | |
189 | { | |
190 | if (AliHBTParticle::fgDebug > 2) | |
191 | Info("ReadNext","Particle rejected since total bayessian PID probab. is zero."); | |
192 | continue; | |
193 | } | |
7be13e68 | 194 | |
bed069a4 | 195 | for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc; |
7be13e68 | 196 | |
bed069a4 | 197 | if (AliHBTParticle::fgDebug > 4) |
198 | { | |
199 | Info("ReadNext","###########################################################################"); | |
200 | Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]); | |
201 | Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]); | |
202 | TString msg("Pid list got from track:"); | |
203 | for (Int_t s = 0;s<kNSpecies;s++) | |
204 | { | |
205 | msg+="\n "; | |
206 | msg+=s; | |
207 | msg+="("; | |
3c69895b | 208 | msg+=charge*GetSpeciesPdgCode((ESpecies)s); |
bed069a4 | 209 | msg+="): "; |
210 | msg+=w[s]; | |
211 | msg+=" ("; | |
212 | msg+=pidtable[s]; | |
213 | msg+=")"; | |
214 | } | |
215 | Info("ReadNext","%s",msg.Data()); | |
216 | }//if (AliHBTParticle::fgDebug>4) | |
7be13e68 | 217 | |
bed069a4 | 218 | for (Int_t s = 0; s<kNSpecies; s++) |
219 | { | |
dd82cadc | 220 | Float_t pp = w[s]; |
221 | if (pp == 0.0) continue; | |
7be13e68 | 222 | |
3c69895b | 223 | Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s); |
bed069a4 | 224 | if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type |
7be13e68 | 225 | |
bed069a4 | 226 | Double_t mass = pdgdb->GetParticle(pdgcode)->Mass(); |
227 | 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 | 228 | |
bed069a4 | 229 | AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, |
230 | mom[0], mom[1], mom[2], tEtot, | |
231 | pos[0], pos[1], pos[2], 0.); | |
232 | //copy probabilitis of other species (if not zero) | |
233 | for (Int_t k = 0; k<kNSpecies; k++) | |
234 | { | |
235 | if (k == s) continue; | |
236 | if (w[k] == 0.0) continue; | |
3c69895b | 237 | track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]); |
bed069a4 | 238 | } |
239 | if(Pass(track))//check if meets all criteria of any of our cuts | |
240 | //if it does not delete it and take next good track | |
241 | { | |
242 | delete track; | |
243 | continue; | |
244 | } | |
245 | ||
246 | fTracksEvent->AddParticle(track); | |
247 | if (particle) fParticlesEvent->AddParticle(particle); | |
248 | keeppart = kTRUE; | |
249 | ||
250 | if (AliHBTParticle::fgDebug > 4 ) | |
251 | { | |
252 | Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode); | |
253 | track->Print(); | |
254 | } | |
255 | }//for (Int_t s = 0; s<kNSpecies; s++) | |
256 | ||
257 | if (keeppart == kFALSE) delete particle;//particle was not stored in event | |
258 | ||
259 | }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks | |
260 | ||
261 | Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).", | |
dd82cadc | 262 | fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(), |
bed069a4 | 263 | fNEventsRead,fCurrentEvent,fCurrentDir); |
7be13e68 | 264 | |
bed069a4 | 265 | fCurrentEvent++; |
266 | fNEventsRead++; | |
267 | delete esd; | |
268 | return 0;//success -> read one event | |
269 | }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array | |
270 | ||
271 | return 1; //no more directories to read | |
272 | } | |
273 | /**********************************************************/ | |
274 | void AliHBTReaderESD::Rewind() | |
275 | { | |
7e35966a | 276 | //rewinds reading |
bed069a4 | 277 | delete fFile; |
278 | fFile = 0; | |
279 | delete fRunLoader; | |
280 | fCurrentDir = 0; | |
281 | fNEventsRead = 0; | |
282 | fCurrentEvent++; | |
bed069a4 | 283 | } |
7be13e68 | 284 | /**********************************************************/ |
bed069a4 | 285 | |
7be13e68 | 286 | TFile* AliHBTReaderESD::OpenFile(Int_t n) |
287 | { | |
bed069a4 | 288 | //opens fFile with kine tree |
7be13e68 | 289 | |
290 | const TString& dirname = GetDirName(n); | |
291 | if (dirname == "") | |
292 | { | |
293 | Error("OpenFiles","Can not get directory name"); | |
294 | return 0x0; | |
295 | } | |
296 | TString filename = dirname +"/"+ fESDFileName; | |
297 | TFile *ret = TFile::Open(filename.Data()); | |
298 | ||
299 | if ( ret == 0x0) | |
300 | { | |
bed069a4 | 301 | Error("OpenFiles","Can't open fFile %s",filename.Data()); |
7be13e68 | 302 | return 0x0; |
303 | } | |
304 | if (!ret->IsOpen()) | |
305 | { | |
bed069a4 | 306 | Error("OpenFiles","Can't open fFile %s",filename.Data()); |
7be13e68 | 307 | return 0x0; |
308 | } | |
309 | ||
bed069a4 | 310 | if (fReadParticles) |
311 | { | |
312 | fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName); | |
313 | if (fRunLoader == 0x0) | |
314 | { | |
315 | Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data()); | |
316 | delete ret; | |
317 | return 0x0; | |
318 | } | |
319 | fRunLoader->LoadHeader(); | |
320 | if (fRunLoader->LoadKinematics()) | |
321 | { | |
322 | Error("Next","Error occured while loading kinematics."); | |
323 | delete fRunLoader; | |
324 | delete ret; | |
325 | return 0x0; | |
326 | } | |
327 | } | |
7be13e68 | 328 | return ret; |
329 | } | |
330 | /**********************************************************/ | |
331 | ||
332 | Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron | |
333 | { | |
334 | //returns pdg code from the PID index | |
335 | //ask jura about charge | |
336 | switch (spec) | |
337 | { | |
338 | case kESDElectron: | |
339 | return kPositron; | |
340 | break; | |
341 | case kESDMuon: | |
342 | return kMuonPlus; | |
343 | break; | |
344 | case kESDPion: | |
345 | return kPiPlus; | |
346 | break; | |
347 | case kESDKaon: | |
348 | return kKPlus; | |
349 | break; | |
350 | case kESDProton: | |
351 | return kProton; | |
352 | break; | |
353 | default: | |
354 | ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec); | |
355 | break; | |
356 | } | |
357 | return 0; | |
358 | } |