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 }