TOF geometry updating (addition of AliTOFGeometry)
[u/mrichter/AliRoot.git] / TOF / AliTOFpidESD.cxx
CommitLineData
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"
3f83f224 24#include "TError.h"
c630aafd 25
26#include "AliTOFpidESD.h"
27#include "AliESD.h"
28#include "AliESDtrack.h"
29#include "AliTOFdigit.h"
30
0f4a7374 31#include <TGeant3.h> // For DtoM and InitGeo, should
32#include "../STRUCT/AliBODY.h" // be removed after we switch to
33#include "../STRUCT/AliFRAMEv2.h" // AliTOGeometry::GetPos
c630aafd 34#include "AliTOFv2FHoles.h"
35
3f83f224 36#include <stdlib.h>
37
c630aafd 38
39ClassImp(AliTOFpidESD)
40
0f4a7374 41//_________________________________________________________________________
c630aafd 42static Int_t InitGeo() {
43 //gSystem->Load("libgeant321");
44 //new TGeant3("C++ Interface to Geant3");
45
3f83f224 46 AliBODY *body = new AliBODY("BODY", "Alice envelop");
47 body->CreateGeometry();
c630aafd 48
3f83f224 49 AliFRAMEv2 *frame = new AliFRAMEv2("FRAME", "Space Frame");
50 frame->SetHoles(1);
51 frame->CreateGeometry();
c630aafd 52
3f83f224 53 AliTOF *tof = new AliTOFv2FHoles("TOF", "TOF with Holes");
54 tof->CreateGeometry();
c630aafd 55
56 return 0;
57}
58
59//_________________________________________________________________________
c630aafd 60static Int_t DtoM(Int_t *dig, Float_t *g) {
3f83f224 61 const Int_t kMAX=13;
62 Int_t lnam[kMAX],lnum[kMAX];
c630aafd 63
64 TGeant3 *geant3=(TGeant3*)gMC;
3f83f224 65 if (!geant3) {::Info("DtoM","no geant3 found !"); return 1;}
c630aafd 66
67 strncpy((Char_t*)(lnam+0),"ALIC",4); lnum[0]=1;
68 strncpy((Char_t*)(lnam+1),"B077",4); lnum[1]=1;
69
3f83f224 70//4 padx if z<=0 then ...
0f4a7374 71if ((dig[1]==4)||(dig[1]==3)) dig[4]=AliTOFGeometry::NpadX()-dig[4];
3f83f224 72else if (dig[1]==2) {
0f4a7374 73 if (dig[2]>7) dig[4]=AliTOFGeometry::NpadX()-dig[4];
3f83f224 74 else if (dig[2]==7) {
0f4a7374 75 if (dig[3]==1) dig[4]=AliTOFGeometry::NpadX()-dig[4];
3f83f224 76 else dig[4]+=1;
77 } else dig[4]+=1;
78} else dig[4]+=1;
79
80//3 padz
0f4a7374 81if ((dig[1]==3)||(dig[1]==4)) dig[3]=AliTOFGeometry::NpadZ()-dig[3];
3f83f224 82else dig[3]+=1;
83
84//2 strip
0f4a7374 85if (dig[1]==0) dig[2]=AliTOFGeometry::NStripC()-dig[2];
86else if (dig[1]==1) dig[2]=AliTOFGeometry::NStripB()-dig[2];
3f83f224 87else dig[2]+=1;
88
89//1 module
0f4a7374 90dig[1]=AliTOFGeometry::NPlates()-dig[1];
3f83f224 91
92//0 sector
0f4a7374 93dig[0]+=(AliTOFGeometry::NSectors()/2); dig[0]%=AliTOFGeometry::NSectors();
3f83f224 94dig[0]+=1;
95
c630aafd 96 //sector
97 switch (dig[0]) {
98 case 1:
99 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=9;
100 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
101 break;
102 case 2:
103 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=10;
104 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
105 break;
106 case 3:
107 strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=4;
108 strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
109 if (dig[1]>3) lnum[3]=2;
110 break;
111 case 4:
112 strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=5;
113 strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
114 if (dig[1]>3) lnum[3]=2;
115 break;
116 case 5:
117 strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=1;
118 strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
119 if (dig[1]>3) lnum[3]=2;
120 break;
121 case 6:
122 strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=2;
123 strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
124 if (dig[1]>3) lnum[3]=2;
125 break;
126 case 7:
127 strncpy((Char_t*)(lnam+2),"B074",4); lnum[2]=3;
128 strncpy((Char_t*)(lnam+3),"BTO2",4); lnum[3]=1;
129 if (dig[1]>3) lnum[3]=2;
130 break;
131 case 8:
132 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=1;
133 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
134 break;
135 case 9:
136 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=2;
137 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
138 break;
139 case 10:
140 strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=1;
141 strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1;
142 if (dig[1]>4) lnum[3]=2;
143 break;
144 case 11:
145 strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=2;
146 strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1;
147 if (dig[1]>4) lnum[3]=2;
148 break;
149 case 12:
150 strncpy((Char_t*)(lnam+2),"B075",4); lnum[2]=3;
151 strncpy((Char_t*)(lnam+3),"BTO3",4); lnum[3]=1;
152 if (dig[1]>4) lnum[3]=2;
153 break;
154 case 13:
155 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=3;
156 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
157 break;
158 case 14:
159 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=4;
160 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
161 break;
162 case 15:
163 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=5;
164 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
165 break;
166 case 16:
167 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=6;
168 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
169 break;
170 case 17:
171 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=7;
172 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
173 break;
174 case 18:
175 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=8;
176 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
177 break;
178 default:
3f83f224 179 ::Info("DtoM","Wrong sector number : %d !",dig[0]);
c630aafd 180 return 2;
181 }
182 //module
183 switch (dig[1]) {
184 case 1:
185 strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1;
186 if (dig[0]==1 || dig[0]==2 ) lnum[4]=2;
187 if (dig[0]>=8 && dig[0]<=9) lnum[4]=2;
188 if (dig[0]>=13 && dig[0]<=18) lnum[4]=2;
189 strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0;
190 break;
191 case 2:
192 strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1;
193 if (dig[0]==1 || dig[0]==2 ) lnum[4]=2;
194 if (dig[0]>=8 && dig[0]<=9) lnum[4]=2;
195 if (dig[0]>=13 && dig[0]<=18) lnum[4]=2;
196 strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0;
197 break;
198 case 3:
199 strncpy((Char_t*)(lnam+4),"FTOA",4); lnum[4]=0;
200 strncpy((Char_t*)(lnam+5),"FLTA",4); lnum[5]=0;
201 break;
202 case 4:
203 strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1;
204 strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0;
205 break;
206 case 5:
207 strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1;
208 strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0;
209 break;
210 default:
3f83f224 211 ::Info("DtoM","Wrong module number : %d !",dig[1]);
c630aafd 212 return 3;
213 }
214
215 //strip
3f83f224 216 if (dig[2]<0/*1*/ || dig[2]>19) {
217 ::Info("DtoM","Wrong strip number : %d !",dig[2]);
c630aafd 218 return 3;
219 }
220 strncpy((Char_t*)(lnam+6),"FSTR",4); lnum[6]=dig[2];
221 strncpy((Char_t*)(lnam+7),"FSEN",4); lnum[7]=0;
222
223 //divz
224 if (dig[3]<1 || dig[3]>2) {
3f83f224 225 ::Info("DtoM","Wrong z-division number : %d !",dig[3]);
c630aafd 226 return 3;
227 }
228 strncpy((Char_t*)(lnam+8),"FSEZ",4); lnum[8]=dig[3];
229
230 //divx
231 if (dig[4]<1 || dig[4]>48) {
3f83f224 232 ::Info("DtoM","Wrong x-division number : %d !",dig[4]);
c630aafd 233 return 3;
234 }
235 strncpy((Char_t*)(lnam+9),"FSEX",4); lnum[9]=dig[4];
236 strncpy((Char_t*)(lnam+10),"FPAD",4); lnum[10]=0;
237
238 Gcvolu_t *gcvolu=geant3->Gcvolu(); gcvolu->nlevel=0;
239 Int_t err=geant3->Glvolu(11,lnam,lnum); //11-th level
3f83f224 240 if (err) {
241 ::Info("DtoM","Wrong volume !"); return 3;
242 }
c630aafd 243 Float_t l[3]={0.,0.,0}; geant3->Gdtom(l,g,1);
244 return 0;
245}
246
0f4a7374 247
248//_________________________________________________________________________
249AliTOFpidESD::AliTOFpidESD(Double_t *param) throw (const Char_t *) {
250 //
251 // The main constructor
252 //
253 if (InitGeo()) throw "AliTOFpidESD: can not initialize the geometry !\n";
254
255
256 if (fTOFGeometry) delete fTOFGeometry;
257 fTOFGeometry = new AliTOFGeometry();
258
259 fR=378.;
260 fDy=fTOFGeometry->XPad(); fDz=fTOFGeometry->ZPad();
261 fN=0; fEventN=0;
262
263 fSigma=param[0];
264 fRange=param[1];
265
266}
267
268
269//_________________________________________________________________________
270
271//_________________________________________________________________________
c630aafd 272Int_t AliTOFpidESD::LoadClusters(const TFile *df) {
0f4a7374 273
c630aafd 274 //--------------------------------------------------------------------
275 //This function loads the TOF clusters
276 //--------------------------------------------------------------------
0f4a7374 277
c630aafd 278 if (!((TFile *)df)->IsOpen()) {
279 Error("LoadClusters","file with the TOF digits has not been open !\n");
280 return 1;
281 }
282
283 Char_t name[100];
284 sprintf(name,"TreeD%d",GetEventNumber());
285 TTree *dTree=(TTree*)((TFile *)df)->Get(name);
286 if (!dTree) {
287 Error("LoadClusters"," can't get the tree with digits !\n");
288 return 1;
289 }
290 TBranch *branch=dTree->GetBranch("TOF");
291 if (!branch) {
292 Error("LoadClusters"," can't get the branch with the TOF digits !\n");
293 return 1;
294 }
295
296 TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
297 branch->SetAddress(&digits);
298
299 dTree->GetEvent(0);
300 Int_t nd=digits->GetEntriesFast();
301 Info("LoadClusters","number of digits: %d",nd);
302
303 for (Int_t i=0; i<nd; i++) {
304 AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
305 Int_t dig[5]; Float_t g[3];
306 dig[0]=d->GetSector();
307 dig[1]=d->GetPlate();
308 dig[2]=d->GetStrip();
309 dig[3]=d->GetPadz();
310 dig[4]=d->GetPadx();
311
0f4a7374 312 // fTOFGeometry->GetPos(dig,g); // uncomment this
c630aafd 313 DtoM(dig,g);
314
315 Double_t h[5];
316 h[0]=TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
317 h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2];
318 h[3]=d->GetTdc(); h[4]=d->GetAdc();
319
3f83f224 320 AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
c630aafd 321 InsertCluster(cl);
322 }
323
324 delete dTree;
325 return 0;
326}
327
0f4a7374 328//_________________________________________________________________________
c630aafd 329Int_t AliTOFpidESD::LoadClusters(TTree *dTree) {
330 //--------------------------------------------------------------------
331 //This function loads the TOF clusters
332 //--------------------------------------------------------------------
333 TBranch *branch=dTree->GetBranch("TOF");
334 if (!branch) {
335 Error("LoadClusters"," can't get the branch with the TOF digits !\n");
336 return 1;
337 }
338
339 TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
340 branch->SetAddress(&digits);
341
342 dTree->GetEvent(0);
343 Int_t nd=digits->GetEntriesFast();
344 Info("LoadClusters","number of digits: %d",nd);
345
346 for (Int_t i=0; i<nd; i++) {
347 AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
348 Int_t dig[5]; Float_t g[3];
349 dig[0]=d->GetSector();
350 dig[1]=d->GetPlate();
351 dig[2]=d->GetStrip();
352 dig[3]=d->GetPadz();
353 dig[4]=d->GetPadx();
354
0f4a7374 355 //fTOFGeometry->GetPos(dig,g);
c630aafd 356 DtoM(dig,g);
357
358 Double_t h[5];
359 h[0]=TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
360 h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2];
361 h[3]=d->GetTdc(); h[4]=d->GetAdc();
362
3f83f224 363 AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
c630aafd 364 InsertCluster(cl);
365 }
366
367 return 0;
368}
369
0f4a7374 370//_________________________________________________________________________
c630aafd 371void AliTOFpidESD::UnloadClusters() {
372 //--------------------------------------------------------------------
373 //This function unloads TOF clusters
374 //--------------------------------------------------------------------
375 for (Int_t i=0; i<fN; i++) delete fClusters[i];
376 fN=0;
377}
378
0f4a7374 379//_________________________________________________________________________
c630aafd 380Int_t AliTOFpidESD::InsertCluster(AliTOFcluster *c) {
381 //--------------------------------------------------------------------
382 //This function adds a cluster to the array of clusters sorted in Z
383 //--------------------------------------------------------------------
384 if (fN==kMaxCluster) {
385 Error("InsertCluster","Too many clusters !\n");
386 return 1;
387 }
388
389 if (fN==0) {fClusters[fN++]=c; return 0;}
390 Int_t i=FindClusterIndex(c->GetZ());
391 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTOFcluster*));
392 fClusters[i]=c; fN++;
393
394 return 0;
395}
396
0f4a7374 397//_________________________________________________________________________
c630aafd 398Int_t AliTOFpidESD::FindClusterIndex(Double_t z) const {
399 //--------------------------------------------------------------------
400 // This function returns the index of the nearest cluster
401 //--------------------------------------------------------------------
402 if (fN==0) return 0;
403 if (z <= fClusters[0]->GetZ()) return 0;
404 if (z > fClusters[fN-1]->GetZ()) return fN;
405 Int_t b=0, e=fN-1, m=(b+e)/2;
406 for (; b<e; m=(b+e)/2) {
407 if (z > fClusters[m]->GetZ()) b=m+1;
408 else e=m;
409 }
410 return m;
411}
412
3f83f224 413static int cmp(const void *p1, const void *p2) {
414 AliESDtrack *t1=*((AliESDtrack**)p1);
415 AliESDtrack *t2=*((AliESDtrack**)p2);
416 Double_t c1[15]; t1->GetExternalCovariance(c1);
417 Double_t c2[15]; t2->GetExternalCovariance(c2);
418 if (c1[0]*c1[2] <c2[0]*c2[2]) return -1;
419 if (c1[0]*c1[2]==c2[0]*c2[2]) return 0;
420 return 1;
421}
422
c630aafd 423//_________________________________________________________________________
424Int_t AliTOFpidESD::MakePID(AliESD *event)
425{
426 //
427 // This function calculates the "detector response" PID probabilities
428 // Just for a bare hint...
429
3f83f224 430 static const Double_t kMasses[]={
c630aafd 431 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
432 };
433
434 Int_t ntrk=event->GetNumberOfTracks();
3f83f224 435 AliESDtrack **tracks=new AliESDtrack*[ntrk];
436
437 Int_t i;
438 for (i=0; i<ntrk; i++) {
c630aafd 439 AliESDtrack *t=event->GetTrack(i);
3f83f224 440 tracks[i]=t;
441 }
442 qsort(tracks,ntrk,sizeof(AliESDtrack*),cmp);
443
444 Int_t nmatch=0;
445 for (i=0; i<ntrk; i++) {
446 AliESDtrack *t=tracks[i];
c630aafd 447
448 if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue;
c630aafd 449
450 Double_t x,par[5]; t->GetExternalParameters(x,par);
451 Double_t cov[15]; t->GetExternalCovariance(cov);
452
3f83f224 453Double_t dphi=(5*TMath::Sqrt(cov[0]) + 0.5*fDy + 2.5*TMath::Abs(par[2]))/fR;
454Double_t dz=5*TMath::Sqrt(cov[2]) + 0.5*fDz + 2.5*TMath::Abs(par[3]);
c630aafd 455
456 Double_t phi=TMath::ATan2(par[0],x) + t->GetAlpha();
457 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
458 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
459 Double_t z=par[1];
460
461 Double_t d2max=1000.;
462 Int_t index=-1;
463 for (Int_t k=FindClusterIndex(z-dz); k<fN; k++) {
464 AliTOFcluster *c=fClusters[k];
465
466 if (c->GetZ() > z+dz) break;
467 if (c->IsUsed()) continue;
468
3f83f224 469 Double_t dph=TMath::Abs(c->GetPhi()-phi);
470 if (dph>TMath::Pi()) dph-=TMath::Pi();
471 if (dph>dphi) continue;
c630aafd 472
3f83f224 473 Double_t d2=dph*dph*fR*fR + (c->GetZ()-z)*(c->GetZ()-z);
c630aafd 474 if (d2 > d2max) continue;
475
476 d2max=d2;
477 index=k;
478 }
479
480 if (index<0) {
481 //Info("MakePID","matching failed ! %d",TMath::Abs(t->GetLabel()));
482 continue;
483 }
484
485 nmatch++;
486
487 AliTOFcluster *c=fClusters[index];
3f83f224 488 c->Use();
489
490 Double_t tof=50*c->GetTDC()+32; // in ps
491 t->SetTOFsignal(tof);
492 t->SetTOFcluster(c->GetIndex());
493
494 if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
c630aafd 495
496 Double_t time[10]; t->GetIntegratedTimes(time);
3f83f224 497
498 //track length correction
499 Double_t rc=TMath::Sqrt(c->GetR()*c->GetR() + c->GetZ()*c->GetZ());
500 Double_t rt=TMath::Sqrt(x*x + par[0]*par[0] + par[1]*par[1]);
501 Double_t dlt=rc-rt;
c630aafd 502
503 Double_t p[10];
504 Double_t mom=t->GetP();
505 for (Int_t j=0; j<AliESDtrack::kSPECIES; j++) {
3f83f224 506 Double_t mass=kMasses[j];
507 Double_t dpp=0.01; //mean relative pt resolution;
508 if (mom>0.5) dpp=0.01*mom;
509 Double_t sigma=dpp*time[j]/(1.+ mom*mom/(mass*mass));
c630aafd 510 sigma=TMath::Sqrt(sigma*sigma + fSigma*fSigma);
3f83f224 511
512 time[j]+=dlt/3e-2*TMath::Sqrt(mom*mom+mass*mass)/mom;
513
c630aafd 514 if (TMath::Abs(tof-time[j]) > fRange*sigma) {
515 p[j]=TMath::Exp(-0.5*fRange*fRange);
516 continue;
517 }
518 p[j]=TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sigma*sigma));
519 }
3f83f224 520 /*
521 if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) {
522 cerr<<"Wrong matching: "<<t->GetLabel()<<endl;
523 continue;
524 }
525 */
c630aafd 526 t->SetTOFpid(p);
527
c630aafd 528 }
529
530 Info("MakePID","Number of matched ESD track: %d",nmatch);
531
3f83f224 532 delete[] tracks;
533
c630aafd 534 return 0;
535}
0f4a7374 536