Replacing AliTOFv2FHoles by AliTOFv4T0
[u/mrichter/AliRoot.git] / TOF / AliTOFpidESD.cxx
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 #include "TError.h"
25
26 #include "AliTOFpidESD.h"
27 #include "AliESD.h"
28 #include "AliESDtrack.h"
29 #include "AliTOFdigit.h"
30
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
34 #include "AliTOFv4T0.h"
35
36 #include <stdlib.h>
37
38
39 ClassImp(AliTOFpidESD)
40
41 //_________________________________________________________________________
42 static Int_t InitGeo() {
43   //gSystem->Load("libgeant321");
44   //new TGeant3("C++ Interface to Geant3");
45
46   AliBODY *body = new AliBODY("BODY", "Alice envelop");
47   body->CreateGeometry();
48
49   AliFRAMEv2 *frame = new AliFRAMEv2("FRAME", "Space Frame");
50   frame->SetHoles(1);
51   frame->CreateGeometry();
52
53   AliTOF *tof = new AliTOFv4T0("TOF", "TOF with Holes");
54   tof->CreateGeometry();
55
56   return 0;
57
58
59 //_________________________________________________________________________
60 static Int_t DtoM(Int_t *dig, Float_t *g) {
61   const Int_t kMAX=13;
62   Int_t lnam[kMAX],lnum[kMAX];
63
64   TGeant3 *geant3=(TGeant3*)gMC;
65   if (!geant3) {::Info("DtoM","no geant3 found !"); return 1;}
66
67   strncpy((Char_t*)(lnam+0),"ALIC",4); lnum[0]=1;
68   strncpy((Char_t*)(lnam+1),"B077",4); lnum[1]=1;
69
70 //4 padx  if z<=0 then ...
71 if ((dig[1]==4)||(dig[1]==3)) dig[4]=AliTOFGeometry::NpadX()-dig[4];
72 else if (dig[1]==2) {
73    if (dig[2]>7) dig[4]=AliTOFGeometry::NpadX()-dig[4];
74    else if (dig[2]==7) {
75       if (dig[3]==1) dig[4]=AliTOFGeometry::NpadX()-dig[4];
76       else dig[4]+=1; 
77    } else dig[4]+=1; 
78 } else dig[4]+=1;
79
80 //3 padz
81 if ((dig[1]==3)||(dig[1]==4)) dig[3]=AliTOFGeometry::NpadZ()-dig[3];
82 else dig[3]+=1;
83
84 //2 strip
85 if (dig[1]==0) dig[2]=AliTOFGeometry::NStripC()-dig[2]; 
86 else if (dig[1]==1) dig[2]=AliTOFGeometry::NStripB()-dig[2];
87 else dig[2]+=1; 
88
89 //1 module
90 dig[1]=AliTOFGeometry::NPlates()-dig[1];
91
92 //0 sector
93 dig[0]+=(AliTOFGeometry::NSectors()/2); dig[0]%=AliTOFGeometry::NSectors();
94 dig[0]+=1;
95
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:
179          ::Info("DtoM","Wrong sector number : %d !",dig[0]);
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:
211         ::Info("DtoM","Wrong module number : %d !",dig[1]);
212         return 3;
213   }
214
215   //strip
216   if (dig[2]<0/*1*/ || dig[2]>19) {
217      ::Info("DtoM","Wrong strip number : %d !",dig[2]);
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) {
225      ::Info("DtoM","Wrong z-division number : %d !",dig[3]);
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) {
232      ::Info("DtoM","Wrong x-division number : %d !",dig[4]);
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
240   if (err) {
241     ::Info("DtoM","Wrong volume !"); return 3;
242   }
243   Float_t l[3]={0.,0.,0}; geant3->Gdtom(l,g,1);
244   return 0;
245 }
246
247
248 //_________________________________________________________________________
249 AliTOFpidESD::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 //_________________________________________________________________________
272 Int_t AliTOFpidESD::LoadClusters(const TFile *df) {
273
274   //--------------------------------------------------------------------
275   //This function loads the TOF clusters
276   //--------------------------------------------------------------------
277
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
312     // fTOFGeometry->GetPos(dig,g);   // uncomment this 
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
320     AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
321     InsertCluster(cl);
322   }  
323
324   delete dTree;
325   return 0;
326 }
327
328 //_________________________________________________________________________
329 Int_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
355     //fTOFGeometry->GetPos(dig,g);
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
363     AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
364     InsertCluster(cl);
365   }  
366
367   return 0;
368 }
369
370 //_________________________________________________________________________
371 void 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
379 //_________________________________________________________________________
380 Int_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
397 //_________________________________________________________________________
398 Int_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
413 static 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
423 //_________________________________________________________________________
424 Int_t AliTOFpidESD::MakePID(AliESD *event)
425 {
426   //
427   //  This function calculates the "detector response" PID probabilities
428   //                Just for a bare hint... 
429
430   static const Double_t kMasses[]={
431     0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
432   };
433
434   Int_t ntrk=event->GetNumberOfTracks();
435   AliESDtrack **tracks=new AliESDtrack*[ntrk];
436
437   Int_t i;
438   for (i=0; i<ntrk; i++) {
439     AliESDtrack *t=event->GetTrack(i);
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];
447
448     if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue;
449
450     Double_t x,par[5]; t->GetExternalParameters(x,par);
451     Double_t cov[15]; t->GetExternalCovariance(cov);
452
453 Double_t dphi=(5*TMath::Sqrt(cov[0]) + 0.5*fDy + 2.5*TMath::Abs(par[2]))/fR; 
454 Double_t dz=5*TMath::Sqrt(cov[2]) + 0.5*fDz + 2.5*TMath::Abs(par[3]);
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
469       Double_t dph=TMath::Abs(c->GetPhi()-phi);
470       if (dph>TMath::Pi()) dph-=TMath::Pi();
471       if (dph>dphi) continue;
472
473       Double_t d2=dph*dph*fR*fR + (c->GetZ()-z)*(c->GetZ()-z);
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];
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;
495
496     Double_t time[10]; t->GetIntegratedTimes(time);
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;
502
503     Double_t p[10];
504     Double_t mom=t->GetP();
505     for (Int_t j=0; j<AliESDtrack::kSPECIES; j++) {
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));
510       sigma=TMath::Sqrt(sigma*sigma + fSigma*fSigma);
511
512       time[j]+=dlt/3e-2*TMath::Sqrt(mom*mom+mass*mass)/mom;
513
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     }
520     /*
521     if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) {
522        cerr<<"Wrong matching: "<<t->GetLabel()<<endl;
523        continue;
524     }
525     */
526     t->SetTOFpid(p);
527
528   }
529
530   Info("MakePID","Number of matched ESD track: %d",nmatch);
531
532   delete[] tracks;
533
534   return 0;
535 }
536