1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 //-----------------------------------------------------------------
23 #include "TClonesArray.h"
25 #include "AliTOFpidESD.h"
27 #include "AliESDtrack.h"
28 #include "AliTOFdigit.h"
31 #include <TVirtualMC.h>
32 #include "../STRUCT/AliBODY.h"
33 #include "../STRUCT/AliFRAMEv2.h"
34 #include "AliTOFv2FHoles.h"
37 ClassImp(AliTOFpidESD)
39 static Int_t InitGeo() {
40 //gSystem->Load("libgeant321");
41 //new TGeant3("C++ Interface to Geant3");
43 AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
44 BODY->CreateGeometry();
46 AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
48 FRAME->CreateGeometry();
50 AliTOF *TOF = new AliTOFv2FHoles("TOF", "TOF with Holes");
51 TOF->CreateGeometry();
56 //_________________________________________________________________________
57 AliTOFpidESD::AliTOFpidESD(Double_t *param) throw (const Char_t *) {
59 // The main constructor
61 if (InitGeo()) throw "AliTOFpidESD: can not initialize the geomtry !\n";
63 fR=376.; fDy=2.5; fDz=3.5; fN=0; fEventN=0;
70 static Int_t DtoM(Int_t *dig, Float_t *g) {
72 Int_t lnam[MAX],lnum[MAX];
74 extern TVirtualMC *gMC;
76 TGeant3 *geant3=(TGeant3*)gMC;
77 if (!geant3) {cerr<<"no geant3 found !\n"; return 1;}
79 strncpy((Char_t*)(lnam+0),"ALIC",4); lnum[0]=1;
80 strncpy((Char_t*)(lnam+1),"B077",4); lnum[1]=1;
85 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=9;
86 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
89 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=10;
90 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
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;
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;
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;
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;
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;
118 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=1;
119 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
122 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=2;
123 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
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;
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;
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;
141 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=3;
142 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
145 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=4;
146 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
149 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=5;
150 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
153 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=6;
154 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
157 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=7;
158 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
161 strncpy((Char_t*)(lnam+2),"B071",4); lnum[2]=8;
162 strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
165 cerr<<"Wrong sector number : "<<dig[0]<<" !\n";
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;
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;
185 strncpy((Char_t*)(lnam+4),"FTOA",4); lnum[4]=0;
186 strncpy((Char_t*)(lnam+5),"FLTA",4); lnum[5]=0;
189 strncpy((Char_t*)(lnam+4),"FTOB",4); lnum[4]=1;
190 strncpy((Char_t*)(lnam+5),"FLTB",4); lnum[5]=0;
193 strncpy((Char_t*)(lnam+4),"FTOC",4); lnum[4]=1;
194 strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0;
197 cerr<<"Wrong module number : "<<dig[1]<<" !\n";
202 if (dig[2]<1 || dig[2]>19) {
203 cerr<<"Wrong strip number : "<<dig[2]<<" !\n";
206 strncpy((Char_t*)(lnam+6),"FSTR",4); lnum[6]=dig[2];
207 strncpy((Char_t*)(lnam+7),"FSEN",4); lnum[7]=0;
210 if (dig[3]<1 || dig[3]>2) {
211 cerr<<"Wrong z-division number : "<<dig[3]<<" !\n";
214 strncpy((Char_t*)(lnam+8),"FSEZ",4); lnum[8]=dig[3];
217 if (dig[4]<1 || dig[4]>48) {
218 cerr<<"Wrong x-division number : "<<dig[4]<<" !\n";
221 strncpy((Char_t*)(lnam+9),"FSEX",4); lnum[9]=dig[4];
222 strncpy((Char_t*)(lnam+10),"FPAD",4); lnum[10]=0;
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);
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");
241 sprintf(name,"TreeD%d",GetEventNumber());
242 TTree *dTree=(TTree*)((TFile *)df)->Get(name);
244 Error("LoadClusters"," can't get the tree with digits !\n");
247 TBranch *branch=dTree->GetBranch("TOF");
249 Error("LoadClusters"," can't get the branch with the TOF digits !\n");
253 TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
254 branch->SetAddress(&digits);
257 Int_t nd=digits->GetEntriesFast();
258 Info("LoadClusters","number of digits: %d",nd);
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();
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();
276 AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks());
284 Int_t AliTOFpidESD::LoadClusters(TTree *dTree) {
285 //--------------------------------------------------------------------
286 //This function loads the TOF clusters
287 //--------------------------------------------------------------------
288 TBranch *branch=dTree->GetBranch("TOF");
290 Error("LoadClusters"," can't get the branch with the TOF digits !\n");
294 TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
295 branch->SetAddress(&digits);
298 Int_t nd=digits->GetEntriesFast();
299 Info("LoadClusters","number of digits: %d",nd);
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();
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();
317 AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks());
324 void AliTOFpidESD::UnloadClusters() {
325 //--------------------------------------------------------------------
326 //This function unloads TOF clusters
327 //--------------------------------------------------------------------
328 for (Int_t i=0; i<fN; i++) delete fClusters[i];
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");
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++;
349 Int_t AliTOFpidESD::FindClusterIndex(Double_t z) const {
350 //--------------------------------------------------------------------
351 // This function returns the index of the nearest cluster
352 //--------------------------------------------------------------------
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;
364 //_________________________________________________________________________
365 Int_t AliTOFpidESD::MakePID(AliESD *event)
368 // This function calculates the "detector response" PID probabilities
369 // Just for a bare hint...
371 static const Double_t masses[]={
372 0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
375 Int_t ntrk=event->GetNumberOfTracks();
377 for (Int_t i=0; i<ntrk; i++) {
378 AliESDtrack *t=event->GetTrack(i);
380 if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue;
381 if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
383 Double_t x,par[5]; t->GetExternalParameters(x,par);
384 Double_t cov[15]; t->GetExternalCovariance(cov);
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.);
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();
394 Double_t d2max=1000.;
396 for (Int_t k=FindClusterIndex(z-dz); k<fN; k++) {
397 AliTOFcluster *c=fClusters[k];
399 if (c->GetZ() > z+dz) break;
400 if (c->IsUsed()) continue;
402 Double_t dphi=TMath::Abs(c->GetPhi()-phi);
403 if (dphi>TMath::Pi()) dphi-=TMath::Pi();
404 if (dphi*dphi > dphi2) continue;
406 Double_t d2=dphi*dphi*fR*fR + (c->GetZ()-z)*(c->GetZ()-z);
407 if (d2 > d2max) continue;
414 //Info("MakePID","matching failed ! %d",TMath::Abs(t->GetLabel()));
420 AliTOFcluster *c=fClusters[index];
422 Double_t time[10]; t->GetIntegratedTimes(time);
423 Double_t tof=50*c->GetTDC(); // in ps
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);
436 p[j]=TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sigma*sigma));
439 //if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) continue;
441 t->SetTOFsignal(tof);
442 t->SetTOFcluster(index);
448 Info("MakePID","Number of matched ESD track: %d",nmatch);