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 #include "AliTOFGeometry.h"
  31
  32 #include <stdlib.h>
  33
  34
  35 ClassImp(AliTOFpidESD)
  36
  37 //_________________________________________________________________________
  38 AliTOFpidESD::AliTOFpidESD(Double_t *param) {
  39   //
  40   //  The main constructor
  41   //
  42   fR=378.;
  43   fDy=AliTOFGeometry::XPad(); fDz=AliTOFGeometry::ZPad();
  44   fN=0; fEventN=0;
  45
  46   fSigma=param[0];
  47   fRange=param[1];
  48
  49 }
  50
  51 //_________________________________________________________________________
  52 Int_t AliTOFpidESD::LoadClusters(TTree *dTree, AliTOFGeometry *geom) {
  53   //--------------------------------------------------------------------
  54   //This function loads the TOF clusters
  55   //--------------------------------------------------------------------
  56   TBranch *branch=dTree->GetBranch("TOF");
  57   if (!branch) {
  58     Error("LoadClusters"," can't get the branch with the TOF digits !\n");
  59     return 1;
  60   }
  61
  62   TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
  63   branch->SetAddress(&digits);
  64
  65   dTree->GetEvent(0);
  66   Int_t nd=digits->GetEntriesFast();
  67   Info("LoadClusters","number of digits: %d",nd);
  68
  69   for (Int_t i=0; i<nd; i++) {
  70     AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
  71     Int_t dig[5]; Float_t g[3];
  72     dig[0]=d->GetSector();
  73     dig[1]=d->GetPlate();
  74     dig[2]=d->GetStrip();
  75     dig[3]=d->GetPadz();
  76     dig[4]=d->GetPadx();
  77
  78     geom->GetPos(dig,g);
  79
  80     Double_t h[5];
  81     h[0]=TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
  82     h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2];
  83     h[3]=d->GetTdc(); h[4]=d->GetAdc();
  84
  85     AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
  86     InsertCluster(cl);
  87   }
  88
  89   return 0;
  90 }
  91
  92 //_________________________________________________________________________
  93 void AliTOFpidESD::UnloadClusters() {
  94   //--------------------------------------------------------------------
  95   //This function unloads TOF clusters
  96   //--------------------------------------------------------------------
  97   for (Int_t i=0; i<fN; i++) delete fClusters[i];
  98   fN=0;
  99 }
 100
 101 //_________________________________________________________________________
 102 Int_t AliTOFpidESD::InsertCluster(AliTOFcluster *c) {
 103   //--------------------------------------------------------------------
 104   //This function adds a cluster to the array of clusters sorted in Z
 105   //--------------------------------------------------------------------
 106   if (fN==kMaxCluster) {
 107     Error("InsertCluster","Too many clusters !\n");
 108     return 1;
 109   }
 110
 111   if (fN==0) {fClusters[fN++]=c; return 0;}
 112   Int_t i=FindClusterIndex(c->GetZ());
 113   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTOFcluster*));
 114   fClusters[i]=c; fN++;
 115
 116   return 0;
 117 }
 118
 119 //_________________________________________________________________________
 120 Int_t AliTOFpidESD::FindClusterIndex(Double_t z) const {
 121   //--------------------------------------------------------------------
 122   // This function returns the index of the nearest cluster
 123   //--------------------------------------------------------------------
 124   if (fN==0) return 0;
 125   if (z <= fClusters[0]->GetZ()) return 0;
 126   if (z > fClusters[fN-1]->GetZ()) return fN;
 127   Int_t b=0, e=fN-1, m=(b+e)/2;
 128   for (; b<e; m=(b+e)/2) {
 129     if (z > fClusters[m]->GetZ()) b=m+1;
 130     else e=m;
 131   }
 132   return m;
 133 }
 134
 135 static int cmp(const void *p1, const void *p2) {
 136   AliESDtrack *t1=*((AliESDtrack**)p1);
 137   AliESDtrack *t2=*((AliESDtrack**)p2);
 138   Double_t c1[15]; t1->GetExternalCovariance(c1);
 139   Double_t c2[15]; t2->GetExternalCovariance(c2);
 140   if (c1[0]*c1[2] <c2[0]*c2[2]) return -1;
 141   if (c1[0]*c1[2]==c2[0]*c2[2]) return 0;
 142   return 1;
 143 }
 144
 145 //_________________________________________________________________________
 146 Int_t AliTOFpidESD::MakePID(AliESD *event)
 147 {
 148   //
 149   //  This function calculates the "detector response" PID probabilities
 150   //                Just for a bare hint...
 151
 152   static const Double_t kMasses[]={
 153     0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
 154   };
 155
 156   Int_t ntrk=event->GetNumberOfTracks();
 157   AliESDtrack **tracks=new AliESDtrack*[ntrk];
 158
 159   Int_t i;
 160   for (i=0; i<ntrk; i++) {
 161     AliESDtrack *t=event->GetTrack(i);
 162     tracks[i]=t;
 163   }
 164   qsort(tracks,ntrk,sizeof(AliESDtrack*),cmp);
 165
 166   Int_t nmatch=0;
 167   for (i=0; i<ntrk; i++) {
 168     AliESDtrack *t=tracks[i];
 169
 170     if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue;
 171     if ((t->GetStatus()&AliESDtrack::kTRDStop)!=0) continue;
 172
 173     Double_t x,par[5]; t->GetExternalParameters(x,par);
 174     Double_t cov[15]; t->GetExternalCovariance(cov);
 175
 176 Double_t dphi=(5*TMath::Sqrt(cov[0]) + 0.5*fDy + 2.5*TMath::Abs(par[2]))/fR;
 177 Double_t dz=5*TMath::Sqrt(cov[2]) + 0.5*fDz + 2.5*TMath::Abs(par[3]);
 178
 179     Double_t phi=TMath::ATan2(par[0],x) + t->GetAlpha();
 180     if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
 181     if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
 182     Double_t z=par[1];
 183
 184     Double_t d2max=1000.;
 185     Int_t index=-1;
 186     for (Int_t k=FindClusterIndex(z-dz); k<fN; k++) {
 187       AliTOFcluster *c=fClusters[k];
 188
 189       if (c->GetZ() > z+dz) break;
 190       if (c->IsUsed()) continue;
 191
 192       Double_t dph=TMath::Abs(c->GetPhi()-phi);
 193       if (dph>TMath::Pi()) dph-=TMath::Pi();
 194       if (dph>dphi) continue;
 195
 196       Double_t d2=dph*dph*fR*fR + (c->GetZ()-z)*(c->GetZ()-z);
 197       if (d2 > d2max) continue;
 198
 199       d2max=d2;
 200       index=k;
 201     }
 202
 203     if (index<0) {
 204       //Info("MakePID","matching failed ! %d",TMath::Abs(t->GetLabel()));
 205        continue;
 206     }
 207
 208     nmatch++;
 209
 210     AliTOFcluster *c=fClusters[index];
 211     c->Use();
 212
 213     Double_t tof=50*c->GetTDC()+32; // in ps
 214     t->SetTOFsignal(tof);
 215     t->SetTOFcluster(c->GetIndex());
 216
 217     if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
 218
 219     Double_t time[10]; t->GetIntegratedTimes(time);
 220
 221     //track length correction
 222     Double_t rc=TMath::Sqrt(c->GetR()*c->GetR() + c->GetZ()*c->GetZ());
 223     Double_t rt=TMath::Sqrt(x*x + par[0]*par[0] + par[1]*par[1]);
 224     Double_t dlt=rc-rt;
 225
 226     Double_t p[10];
 227     Double_t mom=t->GetP();
 228     for (Int_t j=0; j<AliESDtrack::kSPECIES; j++) {
 229
 230       p[j]=1.00/AliESDtrack::kSPECIES;  // don't do anything
 231       if (mom<0.7) continue;            // with the low momenta
 232
 233       Double_t mass=kMasses[j];
 234       Double_t dpp=0.01;      //mean relative pt resolution;
 235       if (mom>0.5) dpp=0.01*mom;
 236       Double_t sigma=dpp*time[j]/(1.+ mom*mom/(mass*mass));
 237       sigma=TMath::Sqrt(sigma*sigma + fSigma*fSigma);
 238
 239       time[j]+=dlt/3e-2*TMath::Sqrt(mom*mom+mass*mass)/mom;
 240
 241       if (TMath::Abs(tof-time[j]) > fRange*sigma) {
 242         p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 243         continue;
 244       }
 245       p[j]=TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sigma*sigma))/sigma;
 246     }
 247     /*
 248     if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) {
 249        cerr<<"Wrong matching: "<<t->GetLabel()<<endl;
 250        continue;
 251     }
 252     */
 253     t->SetTOFpid(p);
 254
 255   }
 256
 257   Info("MakePID","Number of matched ESD track: %d",nmatch);
 258
 259   delete[] tracks;
 260
 261   return 0;
 262 }
 263