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