7be13e68 |
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 <AliESDtrack.h> |
11 | #include <AliESD.h> |
12 | |
13 | #include "AliHBTRun.h" |
14 | #include "AliHBTEvent.h" |
15 | #include "AliHBTParticle.h" |
16 | #include "AliHBTParticleCut.h" |
17 | |
18 | ClassImp(AliHBTReaderESD) |
19 | |
20 | AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename): |
21 | fParticles(new AliHBTRun()), |
22 | fTracks(new AliHBTRun()), |
23 | fESDFileName(esdfilename), |
24 | fIsRead(kFALSE) |
25 | { |
26 | //cosntructor |
27 | if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES)) |
28 | Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra."); |
29 | } |
30 | /********************************************************************/ |
31 | |
32 | AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename): |
33 | AliHBTReader(dirs), |
34 | fParticles(new AliHBTRun()), |
35 | fTracks(new AliHBTRun()), |
36 | fESDFileName(esdfilename), |
37 | fIsRead(kFALSE) |
38 | { |
39 | //cosntructor |
40 | if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES)) |
41 | Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra."); |
42 | } |
43 | /********************************************************************/ |
44 | |
45 | AliHBTReaderESD::~AliHBTReaderESD() |
46 | { |
47 | //desctructor |
48 | delete fParticles; |
49 | delete fTracks; |
50 | } |
51 | /********************************************************************/ |
52 | |
53 | AliHBTEvent* AliHBTReaderESD::GetParticleEvent(Int_t n) |
54 | { |
55 | //returns Nth event with simulated particles |
56 | if (!fIsRead) |
57 | if(Read(fParticles,fTracks)) |
58 | { |
59 | Error("GetParticleEvent","Error in reading"); |
60 | return 0x0; |
61 | } |
62 | return fParticles->GetEvent(n); |
63 | } |
64 | /********************************************************************/ |
65 | |
66 | AliHBTEvent* AliHBTReaderESD::GetTrackEvent(Int_t n) |
67 | { |
68 | //returns Nth event with reconstructed tracks |
69 | if (!fIsRead) |
70 | if(Read(fParticles,fTracks)) |
71 | { |
72 | Error("GetTrackEvent","Error in reading"); |
73 | return 0x0; |
74 | } |
75 | return fTracks->GetEvent(n); |
76 | } |
77 | /********************************************************************/ |
78 | |
79 | Int_t AliHBTReaderESD::GetNumberOfPartEvents() |
80 | { |
81 | //returns number of events of particles |
82 | if (!fIsRead) |
83 | if ( Read(fParticles,fTracks)) |
84 | { |
85 | Error("GetNumberOfPartEvents","Error in reading"); |
86 | return 0; |
87 | } |
88 | return fParticles->GetNumberOfEvents(); |
89 | } |
90 | |
91 | /********************************************************************/ |
92 | Int_t AliHBTReaderESD::GetNumberOfTrackEvents() |
93 | { |
94 | //returns number of events of tracks |
95 | if (!fIsRead) |
96 | if(Read(fParticles,fTracks)) |
97 | { |
98 | Error("GetNumberOfTrackEvents","Error in reading"); |
99 | return 0; |
100 | } |
101 | return fTracks->GetNumberOfEvents(); |
102 | } |
103 | /********************************************************************/ |
104 | |
105 | |
106 | Int_t AliHBTReaderESD::Read(AliHBTRun* particles, AliHBTRun *tracks) |
107 | { |
108 | //reads data and puts put to the particles and tracks objects |
109 | //reurns 0 if everything is OK |
110 | // |
111 | Info("Read",""); |
112 | Int_t totalNevents = 0; |
113 | Int_t currentdir = 0; //number of current directory name is fDirs array |
114 | |
115 | Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track |
116 | Double_t w[kNSpecies]; |
117 | Double_t mom[3];//momentum |
118 | Double_t pos[3];//position |
119 | //****** Tentative particle type "concentrations" |
120 | const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05}; |
121 | |
122 | TDatabasePDG* pdgdb = TDatabasePDG::Instance(); |
123 | if (pdgdb == 0x0) |
124 | { |
125 | Error("Read","Can not get PDG Database Instance."); |
126 | return 1; |
127 | } |
128 | if (!particles) //check if an object is instatiated |
129 | { |
130 | Error("Read"," particles object must instatiated before passing it to the reader"); |
131 | return 1; |
132 | } |
133 | if (!tracks) //check if an object is instatiated |
134 | { |
135 | Error("Read"," tracks object must instatiated before passing it to the reader"); |
136 | return 1; |
137 | } |
138 | particles->Reset();//clear runs == delete all old events |
139 | tracks->Reset(); |
140 | |
141 | Int_t Ndirs; |
142 | if (fDirs) //if array with directories is supplied by user |
143 | { |
144 | Ndirs = fDirs->GetEntries(); //get the number if directories |
145 | } |
146 | else |
147 | { |
148 | Ndirs = 0; //if the array is not supplied read only from current directory |
149 | } |
150 | |
151 | do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./" |
152 | { |
153 | TFile* file = OpenFile(currentdir); |
154 | if (file == 0x0) |
155 | { |
156 | Error("Read","Cannot get File for dir no. %d",currentdir); |
157 | currentdir++; |
158 | continue; |
159 | } |
160 | |
161 | for(Int_t currentEvent = 0; ;currentEvent++)//loop over all events |
162 | { |
163 | TString esdname; |
164 | esdname+=currentEvent; |
165 | AliESD* esd = dynamic_cast<AliESD*>(file->Get(esdname)); |
166 | if (esd == 0x0) |
167 | { |
168 | if (AliHBTParticle::fgDebug > 2 ) |
169 | { |
170 | Info("Read","Can not find ESD object named %s.",esdname.Data()); |
171 | } |
172 | currentdir++; |
173 | break;//we have to assume there is no more ESD objects in the file |
174 | } |
175 | |
176 | Info("Read","Reading Event %d",currentEvent); |
177 | |
178 | Int_t ntr = esd->GetNumberOfTracks(); |
179 | Info("Read","Found %d tracks.",ntr); |
180 | for (Int_t i = 0;i<ntr; i++) |
181 | { |
182 | AliESDtrack *esdtrack = esd->GetTrack(i); |
183 | if (esdtrack == 0x0) |
184 | { |
185 | Error("Read","Can not get track %d", i); |
186 | continue; |
187 | } |
188 | if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) |
189 | { |
190 | if (AliHBTParticle::fgDebug > 2) |
191 | Info("Read","Particle skipped: PID BIT is not set."); |
192 | continue; |
193 | } |
194 | |
195 | esdtrack->GetESDpid(pidtable); |
196 | esdtrack->GetPxPyPz(mom); |
197 | esdtrack->GetXYZ(pos); |
198 | Double_t rc=0.; |
199 | |
200 | //Here we apply Bayes' formula |
201 | for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s]; |
202 | if (rc==0.0) |
203 | { |
204 | if (AliHBTParticle::fgDebug > 2) |
205 | Info("Read","Particle rejected since total bayessian PID probab. is zero."); |
206 | continue; |
207 | } |
208 | |
209 | for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc; |
210 | |
211 | if (AliHBTParticle::fgDebug > 4) |
212 | { |
213 | Info("Read","###########################################################################"); |
214 | Info("Read","Momentum: %f %f %f",mom[0],mom[1],mom[2]); |
215 | Info("Read","Position: %f %f %f",pos[0],pos[1],pos[2]); |
216 | Info("Read","Radius: %f ,rc: %f",TMath::Hypot(pos[0],pos[1]),rc); |
217 | TString msg("Pid list got from track:"); |
218 | for (Int_t s = 0;s<kNSpecies;s++) |
219 | { |
220 | msg+="\n "; |
221 | msg+=s; |
222 | msg+="("; |
223 | msg+=GetSpeciesPdgCode((ESpecies)s); |
224 | msg+="): "; |
225 | msg+=w[s]; |
226 | msg+=" ("; |
227 | msg+=pidtable[s]; |
228 | msg+=")"; |
229 | } |
230 | Info("Read","%s",msg.Data()); |
231 | }//if (AliHBTParticle::fgDebug>4) |
232 | |
233 | for (Int_t s = 0; s<kNSpecies; s++) |
234 | { |
235 | if (w[s] == 0.0) continue; |
236 | |
237 | Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s); |
238 | if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type |
239 | |
240 | Double_t mass = pdgdb->GetParticle(pdgcode)->Mass(); |
241 | Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track |
242 | |
243 | AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, |
244 | mom[0], mom[1], mom[2], tEtot, |
245 | pos[0], pos[1], pos[2], 0.); |
246 | //copy probabilitis of other species (if not zero) |
247 | for (Int_t k = 0; k<kNSpecies; k++) |
248 | { |
249 | if (k == s) continue; |
250 | if (w[k] == 0.0) continue; |
251 | track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]); |
252 | } |
253 | if(Pass(track))//check if meets all criteria of any of our cuts |
254 | //if it does not delete it and take next good track |
255 | { |
256 | delete track; |
257 | continue; |
258 | } |
259 | tracks->AddParticle(totalNevents,track); |
260 | if (AliHBTParticle::fgDebug > 4 ) |
261 | { |
262 | Info("Read","\n\nAdding Particle with incarnation %d",pdgcode); |
263 | track->Print(); |
264 | } |
265 | }//for (Int_t s = 0; s<kNSpecies; s++) |
266 | |
267 | }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks |
268 | totalNevents++; |
269 | }//for(Int_t currentEvent = 0; ;currentEvent++) -- for over ESD in single file |
270 | |
271 | }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array |
272 | |
273 | fIsRead = kTRUE; |
274 | return 0; |
275 | } |
276 | /**********************************************************/ |
277 | |
278 | TFile* AliHBTReaderESD::OpenFile(Int_t n) |
279 | { |
280 | //opens file with kine tree |
281 | |
282 | const TString& dirname = GetDirName(n); |
283 | if (dirname == "") |
284 | { |
285 | Error("OpenFiles","Can not get directory name"); |
286 | return 0x0; |
287 | } |
288 | TString filename = dirname +"/"+ fESDFileName; |
289 | TFile *ret = TFile::Open(filename.Data()); |
290 | |
291 | if ( ret == 0x0) |
292 | { |
293 | Error("OpenFiles","Can't open file %s",filename.Data()); |
294 | return 0x0; |
295 | } |
296 | if (!ret->IsOpen()) |
297 | { |
298 | Error("OpenFiles","Can't open file %s",filename.Data()); |
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 | } |