]>
Commit | Line | Data |
---|---|---|
c630aafd | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | //----------------------------------------------------------------- | |
17 | // Implementation of the TOF PID class | |
18 | // Very naive one... Should be made better by the detector experts... | |
19 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
20 | //----------------------------------------------------------------- | |
21 | #include "TFile.h" | |
22 | #include "TTree.h" | |
23 | #include "TClonesArray.h" | |
24 | ||
25 | #include "AliTOFpidESD.h" | |
26 | #include "AliESD.h" | |
27 | #include "AliESDtrack.h" | |
28 | #include "AliTOFdigit.h" | |
29 | ||
30 | #include <TGeant3.h> | |
31 | #include <TVirtualMC.h> | |
32 | #include "../STRUCT/AliBODY.h" | |
33 | #include "../STRUCT/AliFRAMEv2.h" | |
34 | #include "AliTOFv2FHoles.h" | |
35 | ||
36 | ||
37 | ClassImp(AliTOFpidESD) | |
38 | ||
39 | static Int_t InitGeo() { | |
40 | //gSystem->Load("libgeant321"); | |
41 | //new TGeant3("C++ Interface to Geant3"); | |
42 | ||
43 | AliBODY *BODY = new AliBODY("BODY", "Alice envelop"); | |
44 | BODY->CreateGeometry(); | |
45 | ||
46 | AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame"); | |
47 | FRAME->SetHoles(1); | |
48 | FRAME->CreateGeometry(); | |
49 | ||
50 | AliTOF *TOF = new AliTOFv2FHoles("TOF", "TOF with Holes"); | |
51 | TOF->CreateGeometry(); | |
52 | ||
53 | return 0; | |
54 | } | |
55 | ||
56 | //_________________________________________________________________________ | |
57 | AliTOFpidESD::AliTOFpidESD(Double_t *param) throw (const Char_t *) { | |
58 | // | |
59 | // The main constructor | |
60 | // | |
61 | if (InitGeo()) throw "AliTOFpidESD: can not initialize the geomtry !\n"; | |
62 | ||
63 | fR=376.; fDy=2.5; fDz=3.5; fN=0; fEventN=0; | |
64 | ||
65 | fSigma=param[0]; | |
66 | fRange=param[1]; | |
67 | ||
68 | } | |
69 | ||
70 | static Int_t DtoM(Int_t *dig, Float_t *g) { | |
71 | const Int_t MAX=13; | |
72 | Int_t lnam[MAX],lnum[MAX]; | |
73 | ||
74 | extern TVirtualMC *gMC; | |
75 | ||
76 | TGeant3 *geant3=(TGeant3*)gMC; | |
77 | if (!geant3) {cerr<<"no geant3 found !\n"; return 1;} | |
78 | ||
79 | strncpy((Char_t*)(lnam+0),"ALIC",4); lnum[0]=1; | |
80 | strncpy((Char_t*)(lnam+1),"B077",4); lnum[1]=1; | |
81 | ||
82 | //sector | |
83 | switch (dig[0]) { | |
84 | case 1: | |
85 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=9; | |
86 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
87 | break; | |
88 | case 2: | |
89 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=10; | |
90 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
91 | break; | |
92 | case 3: | |
93 | strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=4; | |
94 | strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; | |
95 | if (dig[1]>3) lnum[3]=2; | |
96 | break; | |
97 | case 4: | |
98 | strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=5; | |
99 | strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; | |
100 | if (dig[1]>3) lnum[3]=2; | |
101 | break; | |
102 | case 5: | |
103 | strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=1; | |
104 | strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; | |
105 | if (dig[1]>3) lnum[3]=2; | |
106 | break; | |
107 | case 6: | |
108 | strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=2; | |
109 | strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; | |
110 | if (dig[1]>3) lnum[3]=2; | |
111 | break; | |
112 | case 7: | |
113 | strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=3; | |
114 | strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1; | |
115 | if (dig[1]>3) lnum[3]=2; | |
116 | break; | |
117 | case 8: | |
118 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=1; | |
119 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
120 | break; | |
121 | case 9: | |
122 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=2; | |
123 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
124 | break; | |
125 | case 10: | |
126 | strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=1; | |
127 | strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1; | |
128 | if (dig[1]>4) lnum[3]=2; | |
129 | break; | |
130 | case 11: | |
131 | strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=2; | |
132 | strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1; | |
133 | if (dig[1]>4) lnum[3]=2; | |
134 | break; | |
135 | case 12: | |
136 | strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=3; | |
137 | strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1; | |
138 | if (dig[1]>4) lnum[3]=2; | |
139 | break; | |
140 | case 13: | |
141 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=3; | |
142 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
143 | break; | |
144 | case 14: | |
145 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=4; | |
146 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
147 | break; | |
148 | case 15: | |
149 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=5; | |
150 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
151 | break; | |
152 | case 16: | |
153 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=6; | |
154 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
155 | break; | |
156 | case 17: | |
157 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=7; | |
158 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
159 | break; | |
160 | case 18: | |
161 | strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=8; | |
162 | strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1; | |
163 | break; | |
164 | default: | |
165 | cerr<<"Wrong sector number : "<<dig[0]<<" !\n"; | |
166 | return 2; | |
167 | } | |
168 | //module | |
169 | switch (dig[1]) { | |
170 | case 1: | |
171 | strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1; | |
172 | if (dig[0]==1 || dig[0]==2 ) lnum[4]=2; | |
173 | if (dig[0]>=8 && dig[0]<=9) lnum[4]=2; | |
174 | if (dig[0]>=13 && dig[0]<=18) lnum[4]=2; | |
175 | strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0; | |
176 | break; | |
177 | case 2: | |
178 | strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1; | |
179 | if (dig[0]==1 || dig[0]==2 ) lnum[4]=2; | |
180 | if (dig[0]>=8 && dig[0]<=9) lnum[4]=2; | |
181 | if (dig[0]>=13 && dig[0]<=18) lnum[4]=2; | |
182 | strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0; | |
183 | break; | |
184 | case 3: | |
185 | strncpy((Char_t*)(lnam+4),"FTOA",4); lnum[4]=0; | |
186 | strncpy((Char_t*)(lnam+5),"FLTA",4); lnum[5]=0; | |
187 | break; | |
188 | case 4: | |
189 | strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1; | |
190 | strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0; | |
191 | break; | |
192 | case 5: | |
193 | strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1; | |
194 | strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0; | |
195 | break; | |
196 | default: | |
197 | cerr<<"Wrong module number : "<<dig[1]<<" !\n"; | |
198 | return 3; | |
199 | } | |
200 | ||
201 | //strip | |
202 | if (dig[2]<1 || dig[2]>19) { | |
203 | cerr<<"Wrong strip number : "<<dig[2]<<" !\n"; | |
204 | return 3; | |
205 | } | |
206 | strncpy((Char_t*)(lnam+6),"FSTR",4); lnum[6]=dig[2]; | |
207 | strncpy((Char_t*)(lnam+7),"FSEN",4); lnum[7]=0; | |
208 | ||
209 | //divz | |
210 | if (dig[3]<1 || dig[3]>2) { | |
211 | cerr<<"Wrong z-division number : "<<dig[3]<<" !\n"; | |
212 | return 3; | |
213 | } | |
214 | strncpy((Char_t*)(lnam+8),"FSEZ",4); lnum[8]=dig[3]; | |
215 | ||
216 | //divx | |
217 | if (dig[4]<1 || dig[4]>48) { | |
218 | cerr<<"Wrong x-division number : "<<dig[4]<<" !\n"; | |
219 | return 3; | |
220 | } | |
221 | strncpy((Char_t*)(lnam+9),"FSEX",4); lnum[9]=dig[4]; | |
222 | strncpy((Char_t*)(lnam+10),"FPAD",4); lnum[10]=0; | |
223 | ||
224 | Gcvolu_t *gcvolu=geant3->Gcvolu(); gcvolu->nlevel=0; | |
225 | Int_t err=geant3->Glvolu(11,lnam,lnum); //11-th level | |
226 | if (err) {cout<<err<<" Wrong volume !\n"; return 3;} | |
227 | Float_t l[3]={0.,0.,0}; geant3->Gdtom(l,g,1); | |
228 | return 0; | |
229 | } | |
230 | ||
231 | Int_t AliTOFpidESD::LoadClusters(const TFile *df) { | |
232 | //-------------------------------------------------------------------- | |
233 | //This function loads the TOF clusters | |
234 | //-------------------------------------------------------------------- | |
235 | if (!((TFile *)df)->IsOpen()) { | |
236 | Error("LoadClusters","file with the TOF digits has not been open !\n"); | |
237 | return 1; | |
238 | } | |
239 | ||
240 | Char_t name[100]; | |
241 | sprintf(name,"TreeD%d",GetEventNumber()); | |
242 | TTree *dTree=(TTree*)((TFile *)df)->Get(name); | |
243 | if (!dTree) { | |
244 | Error("LoadClusters"," can't get the tree with digits !\n"); | |
245 | return 1; | |
246 | } | |
247 | TBranch *branch=dTree->GetBranch("TOF"); | |
248 | if (!branch) { | |
249 | Error("LoadClusters"," can't get the branch with the TOF digits !\n"); | |
250 | return 1; | |
251 | } | |
252 | ||
253 | TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy; | |
254 | branch->SetAddress(&digits); | |
255 | ||
256 | dTree->GetEvent(0); | |
257 | Int_t nd=digits->GetEntriesFast(); | |
258 | Info("LoadClusters","number of digits: %d",nd); | |
259 | ||
260 | for (Int_t i=0; i<nd; i++) { | |
261 | AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i); | |
262 | Int_t dig[5]; Float_t g[3]; | |
263 | dig[0]=d->GetSector(); | |
264 | dig[1]=d->GetPlate(); | |
265 | dig[2]=d->GetStrip(); | |
266 | dig[3]=d->GetPadz(); | |
267 | dig[4]=d->GetPadx(); | |
268 | ||
269 | DtoM(dig,g); | |
270 | ||
271 | Double_t h[5]; | |
272 | h[0]=TMath::Sqrt(g[0]*g[0]+g[1]*g[1]); | |
273 | h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2]; | |
274 | h[3]=d->GetTdc(); h[4]=d->GetAdc(); | |
275 | ||
276 | AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks()); | |
277 | InsertCluster(cl); | |
278 | } | |
279 | ||
280 | delete dTree; | |
281 | return 0; | |
282 | } | |
283 | ||
284 | Int_t AliTOFpidESD::LoadClusters(TTree *dTree) { | |
285 | //-------------------------------------------------------------------- | |
286 | //This function loads the TOF clusters | |
287 | //-------------------------------------------------------------------- | |
288 | TBranch *branch=dTree->GetBranch("TOF"); | |
289 | if (!branch) { | |
290 | Error("LoadClusters"," can't get the branch with the TOF digits !\n"); | |
291 | return 1; | |
292 | } | |
293 | ||
294 | TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy; | |
295 | branch->SetAddress(&digits); | |
296 | ||
297 | dTree->GetEvent(0); | |
298 | Int_t nd=digits->GetEntriesFast(); | |
299 | Info("LoadClusters","number of digits: %d",nd); | |
300 | ||
301 | for (Int_t i=0; i<nd; i++) { | |
302 | AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i); | |
303 | Int_t dig[5]; Float_t g[3]; | |
304 | dig[0]=d->GetSector(); | |
305 | dig[1]=d->GetPlate(); | |
306 | dig[2]=d->GetStrip(); | |
307 | dig[3]=d->GetPadz(); | |
308 | dig[4]=d->GetPadx(); | |
309 | ||
310 | DtoM(dig,g); | |
311 | ||
312 | Double_t h[5]; | |
313 | h[0]=TMath::Sqrt(g[0]*g[0]+g[1]*g[1]); | |
314 | h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2]; | |
315 | h[3]=d->GetTdc(); h[4]=d->GetAdc(); | |
316 | ||
317 | AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks()); | |
318 | InsertCluster(cl); | |
319 | } | |
320 | ||
321 | return 0; | |
322 | } | |
323 | ||
324 | void AliTOFpidESD::UnloadClusters() { | |
325 | //-------------------------------------------------------------------- | |
326 | //This function unloads TOF clusters | |
327 | //-------------------------------------------------------------------- | |
328 | for (Int_t i=0; i<fN; i++) delete fClusters[i]; | |
329 | fN=0; | |
330 | } | |
331 | ||
332 | Int_t AliTOFpidESD::InsertCluster(AliTOFcluster *c) { | |
333 | //-------------------------------------------------------------------- | |
334 | //This function adds a cluster to the array of clusters sorted in Z | |
335 | //-------------------------------------------------------------------- | |
336 | if (fN==kMaxCluster) { | |
337 | Error("InsertCluster","Too many clusters !\n"); | |
338 | return 1; | |
339 | } | |
340 | ||
341 | if (fN==0) {fClusters[fN++]=c; return 0;} | |
342 | Int_t i=FindClusterIndex(c->GetZ()); | |
343 | memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTOFcluster*)); | |
344 | fClusters[i]=c; fN++; | |
345 | ||
346 | return 0; | |
347 | } | |
348 | ||
349 | Int_t AliTOFpidESD::FindClusterIndex(Double_t z) const { | |
350 | //-------------------------------------------------------------------- | |
351 | // This function returns the index of the nearest cluster | |
352 | //-------------------------------------------------------------------- | |
353 | if (fN==0) return 0; | |
354 | if (z <= fClusters[0]->GetZ()) return 0; | |
355 | if (z > fClusters[fN-1]->GetZ()) return fN; | |
356 | Int_t b=0, e=fN-1, m=(b+e)/2; | |
357 | for (; b<e; m=(b+e)/2) { | |
358 | if (z > fClusters[m]->GetZ()) b=m+1; | |
359 | else e=m; | |
360 | } | |
361 | return m; | |
362 | } | |
363 | ||
364 | //_________________________________________________________________________ | |
365 | Int_t AliTOFpidESD::MakePID(AliESD *event) | |
366 | { | |
367 | // | |
368 | // This function calculates the "detector response" PID probabilities | |
369 | // Just for a bare hint... | |
370 | ||
371 | static const Double_t masses[]={ | |
372 | 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613 | |
373 | }; | |
374 | ||
375 | Int_t ntrk=event->GetNumberOfTracks(); | |
376 | Int_t nmatch=0; | |
377 | for (Int_t i=0; i<ntrk; i++) { | |
378 | AliESDtrack *t=event->GetTrack(i); | |
379 | ||
380 | if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue; | |
381 | if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue; | |
382 | ||
383 | Double_t x,par[5]; t->GetExternalParameters(x,par); | |
384 | Double_t cov[15]; t->GetExternalCovariance(cov); | |
385 | ||
386 | Double_t dphi2=3*3*(cov[0] + fDy*fDy/12.)/fR/fR; | |
387 | Double_t dz=3*TMath::Sqrt(cov[2] + fDz*fDz/12.); | |
388 | ||
389 | Double_t phi=TMath::ATan2(par[0],x) + t->GetAlpha(); | |
390 | if (phi<-TMath::Pi()) phi+=2*TMath::Pi(); | |
391 | if (phi>=TMath::Pi()) phi-=2*TMath::Pi(); | |
392 | Double_t z=par[1]; | |
393 | ||
394 | Double_t d2max=1000.; | |
395 | Int_t index=-1; | |
396 | for (Int_t k=FindClusterIndex(z-dz); k<fN; k++) { | |
397 | AliTOFcluster *c=fClusters[k]; | |
398 | ||
399 | if (c->GetZ() > z+dz) break; | |
400 | if (c->IsUsed()) continue; | |
401 | ||
402 | Double_t dphi=TMath::Abs(c->GetPhi()-phi); | |
403 | if (dphi>TMath::Pi()) dphi-=TMath::Pi(); | |
404 | if (dphi*dphi > dphi2) continue; | |
405 | ||
406 | Double_t d2=dphi*dphi*fR*fR + (c->GetZ()-z)*(c->GetZ()-z); | |
407 | if (d2 > d2max) continue; | |
408 | ||
409 | d2max=d2; | |
410 | index=k; | |
411 | } | |
412 | ||
413 | if (index<0) { | |
414 | //Info("MakePID","matching failed ! %d",TMath::Abs(t->GetLabel())); | |
415 | continue; | |
416 | } | |
417 | ||
418 | nmatch++; | |
419 | ||
420 | AliTOFcluster *c=fClusters[index]; | |
421 | ||
422 | Double_t time[10]; t->GetIntegratedTimes(time); | |
423 | Double_t tof=50*c->GetTDC(); // in ps | |
424 | ||
425 | Double_t p[10]; | |
426 | Double_t mom=t->GetP(); | |
427 | for (Int_t j=0; j<AliESDtrack::kSPECIES; j++) { | |
428 | Double_t mass=masses[j]; | |
429 | Double_t dp_p=0.03; //mean relative pt resolution; | |
430 | Double_t sigma=dp_p*time[j]/(1.+ mom*mom/(mass*mass)); | |
431 | sigma=TMath::Sqrt(sigma*sigma + fSigma*fSigma); | |
432 | if (TMath::Abs(tof-time[j]) > fRange*sigma) { | |
433 | p[j]=TMath::Exp(-0.5*fRange*fRange); | |
434 | continue; | |
435 | } | |
436 | p[j]=TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sigma*sigma)); | |
437 | } | |
438 | ||
439 | //if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) continue; | |
440 | ||
441 | t->SetTOFsignal(tof); | |
442 | t->SetTOFcluster(index); | |
443 | t->SetTOFpid(p); | |
444 | ||
445 | c->Use(); | |
446 | } | |
447 | ||
448 | Info("MakePID","Number of matched ESD track: %d",nmatch); | |
449 | ||
450 | return 0; | |
451 | } |