Implementing ESD functionality in the NewIO (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
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 }