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