TOF/AliTOFtrackerMI.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 //                                                                 //
  18 //   AliTOFtracker Class                                           //
  19 //   Task: Perform association of the ESD tracks to TOF Clusters   //
  20 //   and Update ESD track with associated TOF Cluster parameters   //
  21 //                                                                 //
  22 //-----------------------------------------------------------------//
  23
  24 #include <Rtypes.h>
  25
  26 #include "TClonesArray.h"
  27 #include "TGeoManager.h"
  28 #include "TTree.h"
  29 #include "TTreeStream.h"
  30
  31 #include "AliRun.h"
  32 #include "AliESD.h"
  33 #include "AliESDtrack.h"
  34
  35 #include "AliTOFcalib.h"
  36 #include "AliTOFcluster.h"
  37 #include "AliTOFGeometry.h"
  38 #include "AliTOFtrackerMI.h"
  39 #include "AliTOFtrack.h"
  40
  41 extern TGeoManager *gGeoManager;
  42
  43 ClassImp(AliTOFtrackerMI)
  44
  45 //_____________________________________________________________________________
  46 AliTOFtrackerMI::AliTOFtrackerMI(AliTOFGeometry * geom, Double_t parPID[2]):
  47   fGeom(geom),
  48   fTOFpid(new AliTOFpidESD(parPID)),
  49   fHoles(kFALSE),
  50   fN(0),
  51   fNseeds(0),
  52   fNseedsTOF(0),
  53   fngoodmatch(0),
  54   fnbadmatch(0),
  55   fnunmatch(0),
  56   fnmatch(0),
  57   fR(378.),
  58   fTOFHeigth(15.3),
  59   fdCut(3.),
  60   fDx(1.5),
  61   fDy(0),
  62   fDz(0),
  63   fTracks(0x0),
  64   fSeeds(0x0),
  65   fDebugStreamer(0x0)
  66  {
  67   //AliTOFtrackerMI main Ctor
  68
  69   //fHoles=true;
  70   fDy=AliTOFGeometry::XPad();
  71   fDz=AliTOFGeometry::ZPad();
  72   fDebugStreamer = new TTreeSRedirector("TOFdebug.root");
  73   //Init(); // temporary solution to know about Holes/no Holes
  74   fHoles=geom->GetHoles();
  75 }
  76 //_____________________________________________________________________________
  77 AliTOFtrackerMI::AliTOFtrackerMI(const AliTOFtrackerMI &t):
  78   AliTracker(),
  79   fGeom(0x0),
  80   fTOFpid(0x0),
  81   fHoles(kFALSE),
  82   fN(0),
  83   fNseeds(0),
  84   fNseedsTOF(0),
  85   fngoodmatch(0),
  86   fnbadmatch(0),
  87   fnunmatch(0),
  88   fnmatch(0),
  89   fR(378.),
  90   fTOFHeigth(15.3),
  91   fdCut(3.),
  92   fDx(1.5),
  93   fDy(0),
  94   fDz(0),
  95   fTracks(0x0),
  96   fSeeds(0x0),
  97   fDebugStreamer(0x0)
  98  {
  99   //AliTOFtrackerMI copy Ctor
 100
 101   fHoles=t.fHoles;
 102   fNseeds=t.fNseeds;
 103   fNseedsTOF=t.fNseedsTOF;
 104   fngoodmatch=t.fngoodmatch;
 105   fnbadmatch=t.fnbadmatch;
 106   fnunmatch=t.fnunmatch;
 107   fnmatch=t.fnmatch;
 108   fGeom = t.fGeom;
 109   fTOFpid = t.fTOFpid;
 110   fR=t.fR;
 111   fTOFHeigth=t.fTOFHeigth;
 112   fdCut=t.fdCut;
 113   fDy=t.fDy;
 114   fDz=t.fDz;
 115   fDx=1.5;
 116   fSeeds=t.fSeeds;
 117   fTracks=t.fTracks;
 118   fN=t.fN;
 119 }
 120
 121 //_________________________________________________________________________________
 122 AliTOFtrackerMI& AliTOFtrackerMI::operator=(const AliTOFtrackerMI &t)
 123 {
 124   //AliTOFtrackerMI assignment operator
 125
 126   this->fHoles=t.fHoles;
 127   this->fNseeds=t.fNseeds;
 128   this->fNseedsTOF=t.fNseedsTOF;
 129   this->fngoodmatch=t.fngoodmatch;
 130   this->fnbadmatch=t.fnbadmatch;
 131   this->fnunmatch=t.fnunmatch;
 132   this->fnmatch=t.fnmatch;
 133   this->fGeom = t.fGeom;
 134   this->fTOFpid = t.fTOFpid;
 135   this->fR=t.fR;
 136   this->fTOFHeigth=t.fTOFHeigth;
 137   this->fdCut=t.fdCut;
 138   this->fDy=t.fDy;
 139   this->fDz=t.fDz;
 140   this->fDx=t.fDx;
 141   this->fSeeds=t.fSeeds;
 142   this->fTracks=t.fTracks;
 143   this->fN=t.fN;
 144   return *this;
 145
 146 }
 147
 148 //_____________________________________________________________________________
 149 AliTOFtrackerMI::~AliTOFtrackerMI(){
 150   //
 151   //
 152   //
 153   if (fDebugStreamer) {
 154     //fDebugStreamer->Close();
 155     delete fDebugStreamer;
 156   }
 157   delete fTOFpid;
 158 }
 159
 160 //_____________________________________________________________________________
 161 /*
 162 void AliTOFtrackerMI::Init() {
 163
 164 // temporary solution to know about Holes/no Holes, will be implemented as
 165 // an AliTOFGeometry getter
 166
 167   AliModule* frame=gAlice->GetModule("FRAME");
 168
 169   if(!frame) {
 170     AliError("Could Not load FRAME! Assume Frame with Holes ");
 171     fHoles=true;
 172   } else{
 173     if(frame->IsVersion()==1) {fHoles=false;}
 174     else {fHoles=true;}
 175   }
 176 }
 177 */
 178 //_____________________________________________________________________________
 179 Int_t AliTOFtrackerMI::PropagateBack(AliESD* event) {
 180   //
 181   // Gets seeds from ESD event and Match with TOF Clusters
 182   //
 183
 184
 185   //Initialise some counters
 186
 187   fNseeds=0;
 188   fNseedsTOF=0;
 189   fngoodmatch=0;
 190   fnbadmatch=0;
 191   fnunmatch=0;
 192   fnmatch=0;
 193
 194   Int_t ntrk=event->GetNumberOfTracks();
 195   fNseeds = ntrk;
 196   fSeeds= new TClonesArray("AliESDtrack");
 197   TClonesArray &aESDTrack = *fSeeds;
 198
 199
 200   //Load ESD tracks into a local Array of ESD Seeds
 201
 202   for (Int_t i=0; i<fNseeds; i++) {
 203     AliESDtrack *t=event->GetTrack(i);
 204     new(aESDTrack[i]) AliESDtrack(*t);
 205   }
 206
 207   //Prepare ESD tracks candidates for TOF Matching
 208   CollectESD();
 209
 210   //First Step with Strict Matching Criterion
 211   //MatchTracks(kFALSE);
 212
 213   //Second Step with Looser Matching Criterion
 214   //MatchTracks(kTRUE);
 215   MatchTracksMI(kFALSE);  // assign track to clusters
 216   MatchTracksMI(kTRUE);   // assign clusters to esd
 217
 218   Info("PropagateBack","Number of matched tracks: %d",fnmatch);
 219   Info("PropagateBack","Number of good matched tracks: %d",fngoodmatch);
 220   Info("PropagateBack","Number of bad  matched tracks: %d",fnbadmatch);
 221
 222   //Update the matched ESD tracks
 223
 224   for (Int_t i=0; i<ntrk; i++) {
 225     AliESDtrack *t=event->GetTrack(i);
 226     AliESDtrack *seed =(AliESDtrack*)fSeeds->UncheckedAt(i);
 227     if(seed->GetTOFsignal()>0){
 228       t->SetTOFsignal(seed->GetTOFsignal());
 229       t->SetTOFcluster(seed->GetTOFcluster());
 230       Int_t tlab[3];
 231       seed->GetTOFLabel(tlab);
 232       t->SetTOFLabel(tlab);
 233       AliTOFtrack *track = new AliTOFtrack(*seed);
 234       Float_t info[10];
 235       Double_t times[10];
 236       seed->GetTOFInfo(info);
 237       seed->GetIntegratedTimes(times);
 238       t->UpdateTrackParams(track,AliESDtrack::kTOFout);
 239       t->SetIntegratedLength(seed->GetIntegratedLength());
 240       t->SetIntegratedTimes(times);
 241       t->SetTOFsignalToT(seed->GetTOFsignalToT());
 242       t->SetTOFCalChannel(seed->GetTOFCalChannel());
 243       //
 244       t->SetTOFInfo(info);
 245       delete track;
 246     }
 247   }
 248
 249
 250   //Make TOF PID
 251   fTOFpid->MakePID(event);
 252
 253   if (fSeeds) {
 254     fSeeds->Delete();
 255     delete fSeeds;
 256     fSeeds = 0x0;
 257   }
 258   if (fTracks) {
 259     fTracks->Delete();
 260     delete fTracks;
 261     fTracks = 0x0;
 262   }
 263   return 0;
 264
 265 }
 266 //_________________________________________________________________________
 267 void AliTOFtrackerMI::CollectESD() {
 268    //prepare the set of ESD tracks to be matched to clusters in TOF
 269
 270   fTracks= new TClonesArray("AliTOFtrack");
 271   TClonesArray &aTOFTrack = *fTracks;
 272   Int_t c0=0;
 273   Int_t c1=0;
 274   for (Int_t i=0; i<fNseeds; i++) {
 275
 276     AliESDtrack *t =(AliESDtrack*)fSeeds->UncheckedAt(i);
 277     if ((t->GetStatus()&AliESDtrack::kTPCout)==0)continue;
 278
 279     // TRD good tracks, already propagated at 371 cm
 280
 281     AliTOFtrack *track = new AliTOFtrack(*t); // New
 282     Double_t x = track->GetX(); //New
 283
 284     if (((t->GetStatus()&AliESDtrack::kTRDout)!=0 ) &&
 285          ( x >= fGeom->RinTOF()) ){
 286       track->SetSeedIndex(i);
 287       t->UpdateTrackParams(track,AliESDtrack::kTOFout);
 288       new(aTOFTrack[fNseedsTOF]) AliTOFtrack(*track);
 289       fNseedsTOF++;
 290       c0++;
 291       delete track;
 292     }
 293
 294     // Propagate the rest of TPCbp
 295
 296     else {
 297       if(track->PropagateToInnerTOF(fHoles)){ // temporary solution
 298         //      if(track->PropagateToInnerTOF(fGeom->GetHoles())){
 299         track->SetSeedIndex(i);
 300         t->UpdateTrackParams(track,AliESDtrack::kTOFout);
 301         new(aTOFTrack[fNseedsTOF]) AliTOFtrack(*track);
 302         fNseedsTOF++;
 303         c1++;
 304       }
 305       delete track;
 306     }
 307   }
 308   //
 309   //
 310   printf("TRD\tOn\t%d\tOff\t%d\n",c0,c1);
 311
 312   // Sort according uncertainties on track position
 313   fTracks->Sort();
 314
 315 }
 316
 317
 318
 319
 320
 321
 322
 323 //
 324 //
 325 //_________________________________________________________________________
 326 void AliTOFtrackerMI::MatchTracks( Bool_t /*mLastStep*/){
 327   return;
 328 }
 329 //
 330 //
 331 //_________________________________________________________________________
 332 void AliTOFtrackerMI::MatchTracksMI(Bool_t mLastStep){
 333
 334   //Match ESD tracks to clusters in TOF
 335   const Float_t kTofOffset = 26;  // time offset
 336   const Float_t kMinQuality        = -6.; // minimal quality
 337   const Float_t kMaxQualityD       = 1.;  // max delta quality if cluster used
 338   const Float_t kForbiddenR        = 0.1;  // minimal PID according TPC
 339
 340   static const Double_t kMasses[]={
 341     0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
 342   };
 343
 344   Int_t nSteps=(Int_t)(fTOFHeigth/0.1);
 345
 346   AliTOFcalib *calib = new AliTOFcalib(fGeom);
 347
 348   //PH Arrays (moved outside of the loop)
 349   Float_t * trackPos[4];
 350   for (Int_t ii=0; ii<4; ii++) trackPos[ii] = new Float_t[nSteps];
 351   Int_t * clind[6];
 352   for (Int_t ii=0;ii<6;ii++) clind[ii] = new Int_t[fN];
 353
 354   // Some init
 355
 356   Int_t         index[1000];
 357   Float_t       dist3D[1000][6];
 358   Double_t      times[1000][6];
 359   Float_t       mintimedist[1000];
 360   Float_t       likelihood[1000];
 361   Float_t       length[1000];
 362   AliTOFcluster *clusters[1000];
 363   Double_t       tpcpid[5];
 364   dist3D[0][0]=1;
 365
 366   for (Int_t i=0; i<fNseedsTOF; i++) {
 367
 368     AliTOFtrack *track =(AliTOFtrack*)fTracks->UncheckedAt(i);
 369     AliESDtrack *t =(AliESDtrack*)fSeeds->UncheckedAt(track->GetSeedIndex());
 370     Bool_t hasTime = ( (t->GetStatus()& AliESDtrack::kTIME)>0) ? kTRUE:kFALSE;   // did we integrate time
 371     Float_t trdquality = t->GetTRDQuality();
 372     //
 373     // Normalize tpc pid
 374     //
 375     t->GetTPCpid(tpcpid);
 376     Double_t sumpid=0;
 377     for (Int_t ipid=0;ipid<5;ipid++){
 378       sumpid+=tpcpid[ipid];
 379     }
 380     for (Int_t ipid=0;ipid<5;ipid++){
 381       if (sumpid>0) tpcpid[ipid]/=sumpid;
 382       else{
 383         tpcpid[ipid]=0.2;
 384       }
 385     }
 386
 387     if (trdquality<0) continue; // no chance
 388     //
 389     AliTOFtrack *trackTOFin =new AliTOFtrack(*track);
 390     //
 391     //propagat track to the middle of TOF
 392     //
 393     Float_t xs = 378.2;  // should be defined in the TOF geometry
 394     Double_t ymax=xs*TMath::Tan(0.5*AliTOFGeometry::GetAlpha());
 395     Bool_t skip=kFALSE;
 396     Double_t ysect=trackTOFin->GetYat(xs,skip);
 397     if (skip){
 398       xs = 372.;
 399       ymax=xs*TMath::Tan(0.5*AliTOFGeometry::GetAlpha());
 400       ysect=trackTOFin->GetYat(xs,skip);
 401     }
 402     if (ysect > ymax) {
 403       if (!trackTOFin->Rotate(AliTOFGeometry::GetAlpha())) {
 404         continue;
 405       }
 406     } else if (ysect <-ymax) {
 407       if (!trackTOFin->Rotate(-AliTOFGeometry::GetAlpha())) {
 408         continue;
 409       }
 410     }
 411     if(!trackTOFin->PropagateTo(xs)) {
 412       continue;
 413     }
 414     //
 415     // Determine a window around the track
 416     //
 417     Double_t x,par[5];
 418     trackTOFin->GetExternalParameters(x,par);
 419     Double_t cov[15];
 420     trackTOFin->GetExternalCovariance(cov);
 421     Float_t scalefact=3.;
 422     Double_t dphi=
 423       scalefact*
 424       ((5*TMath::Sqrt(cov[0]) + 3.*fDy +10.*TMath::Abs(par[2]))/fR);
 425     Double_t dz=
 426       scalefact*
 427       (5*TMath::Sqrt(cov[2]) + 3.*fDz  +10.*TMath::Abs(par[3]));
 428
 429     Double_t phi=TMath::ATan2(par[0],x) + trackTOFin->GetAlpha();
 430     if (phi<-TMath::Pi())phi+=2*TMath::Pi();
 431     if (phi>=TMath::Pi())phi-=2*TMath::Pi();
 432     Double_t z=par[1];
 433
 434     Int_t nc     =0;
 435     Int_t nfound =0;
 436     // find the clusters in the window of the track
 437
 438     for (Int_t k=FindClusterIndex(z-dz); k<fN; k++) {
 439       AliTOFcluster *c=fClusters[k];
 440       if (c->GetZ() > z+dz) break;
 441       //      if (c->IsUsed()) continue;
 442
 443       Double_t dph=TMath::Abs(c->GetPhi()-phi);
 444       if (dph>TMath::Pi()) dph-=2.*TMath::Pi();
 445       if (TMath::Abs(dph)>dphi) continue;
 446
 447       clind[0][nc] = c->GetDetInd(0);
 448       clind[1][nc] = c->GetDetInd(1);
 449       clind[2][nc] = c->GetDetInd(2);
 450       clind[3][nc] = c->GetDetInd(3);
 451       clind[4][nc] = c->GetDetInd(4);
 452       clind[5][nc] = k;
 453       nc++;
 454     }
 455
 456     //
 457     // select close clusters
 458     //
 459     Double_t mom=t->GetP();
 460     //    Bool_t dump = kTRUE;
 461     for (Int_t icl=0; icl<nc; icl++){
 462       Float_t distances[5];
 463       if (nfound>=1000) break;
 464       index[nfound]=clind[5][icl];
 465       AliTOFcluster *cluster = fClusters[clind[5][icl]];
 466       GetLinearDistances(trackTOFin,cluster, distances);
 467       dist3D[nfound][0] = distances[4];
 468       dist3D[nfound][1] = distances[1];
 469       dist3D[nfound][2] = distances[2];
 470       // cut on distance
 471       if (TMath::Abs(dist3D[nfound][1])>20 || TMath::Abs(dist3D[nfound][2])>20) continue;
 472       //
 473       GetLikelihood(distances[1],distances[2],cov,trackTOFin, dist3D[nfound][3], dist3D[nfound][4]);
 474       //
 475       // cut on likelihood
 476       if (dist3D[nfound][3]*dist3D[nfound][4]<0.00000000000001) continue;  // log likelihood
 477       if (TMath::Log(dist3D[nfound][3]*dist3D[nfound][4])<-9) continue;  // log likelihood
 478       //
 479       clusters[nfound] = cluster;
 480       //
 481       //length  and TOF updates
 482       trackTOFin->GetIntegratedTimes(times[nfound]);
 483       length[nfound] = trackTOFin->GetIntegratedLength();
 484       length[nfound]+=distances[4];
 485       mintimedist[nfound]=1000;
 486       Double_t tof2=AliTOFGeometry::TdcBinWidth()*cluster->GetTDC()+kTofOffset; // in ps
 487       // Float_t tgamma = TMath::Sqrt(cluster->GetR()*cluster->GetR()+cluster->GetZ()*cluster->GetZ())/0.03;  //time for "primary" gamma
 488       //if (trackTOFin->GetPt()<0.7 && TMath::Abs(tgamma-tof2)<350) continue;  // gamma conversion candidate - TEMPORARY
 489       for(Int_t j=0;j<=5;j++){
 490         Double_t mass=kMasses[j];
 491         times[nfound][j]+=distances[4]/3e-2*TMath::Sqrt(mom*mom+mass*mass)/mom;   // add time distance
 492         if ( TMath::Abs(times[nfound][j]-tof2)<mintimedist[nfound] && tpcpid[j]>kForbiddenR){
 493           mintimedist[nfound]=TMath::Abs(times[nfound][j]-tof2);
 494         }
 495       }
 496       //
 497       Float_t   liketime =  TMath::Exp(-mintimedist[nfound]/90.)+0.25*TMath::Exp(-mintimedist[nfound]/180.);
 498       if (!hasTime)  liketime=0.2;
 499       likelihood[nfound] = TMath::Log(dist3D[nfound][3]*dist3D[nfound][4]*liketime);
 500
 501       if (TMath::Log(dist3D[nfound][3]*dist3D[nfound][4])<-1){
 502         if (likelihood[nfound]<-9.) continue;
 503       }
 504       //
 505       nfound++;
 506     }
 507     if (nfound == 0 ) {
 508       fnunmatch++;
 509       delete trackTOFin;
 510       continue;
 511     }
 512     //
 513     //choose the best cluster
 514     //
 515     Float_t quality[1000];
 516     Int_t   index[1000];
 517     //
 518     AliTOFcluster * cgold=0;
 519     Int_t igold =-1;
 520     for (Int_t icl=0;icl<nfound;icl++){
 521       quality[icl] = dist3D[icl][3]*dist3D[icl][4];
 522     }
 523     TMath::Sort(nfound,quality,index,kTRUE);
 524     igold =  index[0];
 525     cgold =  clusters[igold];
 526     if (nfound>1 &&likelihood[index[1]]>likelihood[index[0]]){
 527       if ( quality[index[0]]<quality[index[1]]+0.5){
 528         igold = index[1];
 529         cgold =  clusters[igold];
 530       }
 531     }
 532     //
 533     //
 534     Float_t qualityGold = TMath::Log(0.0000001+(quality[igold]*(0.1+TMath::Exp(-mintimedist[igold]/80.))*(0.2+trdquality)));
 535     if (!mLastStep){
 536       if (cgold->GetQuality()<qualityGold) cgold->SetQuality(qualityGold);
 537       continue;
 538     }
 539     //
 540     if (mLastStep){
 541       //signed better cluster
 542       if (cgold->GetQuality()>qualityGold+kMaxQualityD) continue;
 543       if (2.*qualityGold-cgold->GetQuality()<kMinQuality) continue;
 544     }
 545
 546     Int_t inonfake=-1;
 547
 548     for (Int_t icl=0;icl<nfound;icl++){
 549       if (
 550           (clusters[index[icl]]->GetLabel(0)==TMath::Abs(trackTOFin->GetLabel()))
 551           ||
 552           (clusters[index[icl]]->GetLabel(1)==TMath::Abs(trackTOFin->GetLabel()))
 553           ||
 554           (clusters[index[icl]]->GetLabel(2)==TMath::Abs(trackTOFin->GetLabel()))
 555           ) {
 556         inonfake = icl;
 557         break;
 558       }
 559     }
 560     fnmatch++;;
 561     if (inonfake==0) fngoodmatch++;
 562     else{
 563       fnbadmatch++;
 564     }
 565
 566     Int_t tlab[3];
 567     tlab[0]=cgold->GetLabel(0);
 568     tlab[1]=cgold->GetLabel(1);
 569     tlab[2]=cgold->GetLabel(2);
 570     //    Double_t tof2=25.*cgold->GetTDC()-350; // in ps
 571     Double_t tof2=AliTOFGeometry::TdcBinWidth()*cgold->GetTDC()+kTofOffset; // in ps
 572     Float_t tgamma = TMath::Sqrt(cgold->GetR()*cgold->GetR()+cgold->GetZ()*cgold->GetZ())/0.03;
 573     Float_t info[11]={dist3D[igold][0],dist3D[igold][1],dist3D[igold][2],dist3D[igold][3],dist3D[igold][4],mintimedist[igold],
 574                       -1,tgamma, qualityGold,cgold->GetQuality(),0};
 575     //    GetLinearDistances(trackTOFin,cgold,&info[6]);
 576     if (inonfake>=0){
 577       info[6] = inonfake;
 578       //      info[7] = mintimedist[index[inonfake]];
 579     }
 580     //
 581     //  Store quantities to be used for TOF Calibration
 582     Float_t tToT=cgold->GetToT(); // in ps
 583     t->SetTOFsignalToT(tToT);
 584     Int_t ind[5];
 585     ind[0]=cgold->GetDetInd(0);
 586     ind[1]=cgold->GetDetInd(1);
 587     ind[2]=cgold->GetDetInd(2);
 588     ind[3]=cgold->GetDetInd(3);
 589     ind[4]=cgold->GetDetInd(4);
 590     Int_t calindex = calib->GetIndex(ind);
 591     t->SetTOFCalChannel(calindex);
 592
 593     t->SetTOFInfo(info);
 594     t->SetTOFsignal(tof2);
 595     t->SetTOFcluster(cgold->GetIndex());
 596     AliTOFtrack *trackTOFout = new AliTOFtrack(*t);
 597     trackTOFout->PropagateTo(378.);
 598     t->UpdateTrackParams(trackTOFout,AliESDtrack::kTOFout);
 599     t->SetIntegratedLength(length[igold]);
 600     t->SetIntegratedTimes(times[igold]);
 601     t->SetTOFLabel(tlab);
 602     //
 603     delete trackTOFin;
 604     delete trackTOFout;
 605     //
 606   }
 607   //
 608   //
 609   //
 610   for (Int_t ii=0; ii<4; ii++) delete [] trackPos[ii];
 611   for (Int_t ii=0;ii<6;ii++) delete [] clind[ii];
 612   delete calib;
 613 }
 614
 615
 616
 617
 618 // //_________________________________________________________________________
 619 // Int_t AliTOFtrackerMI::LoadClusters(TTree *dTree) {
 620 //   //--------------------------------------------------------------------
 621 //   //This function loads the TOF clusters
 622 //   //--------------------------------------------------------------------
 623
 624 //   TBranch *branch=dTree->GetBranch("TOF");
 625 //   if (!branch) {
 626 //     AliError(" can't get the branch with the TOF digits !");
 627 //     return 1;
 628 //   }
 629
 630 //   TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
 631 //   branch->SetAddress(&digits);
 632
 633 //   dTree->GetEvent(0);
 634 //   Int_t nd=digits->GetEntriesFast();
 635 //   Info("LoadClusters","number of digits: %d",nd);
 636
 637 //   for (Int_t i=0; i<nd; i++) {
 638 //     AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
 639 //     Int_t dig[5]; Float_t g[3];
 640 //     dig[0]=d->GetSector();
 641 //     dig[1]=d->GetPlate();
 642 //     dig[2]=d->GetStrip();
 643 //     dig[3]=d->GetPadz();
 644 //     dig[4]=d->GetPadx();
 645
 646 //     fGeom->GetPos(dig,g);
 647
 648 //     Double_t h[5];
 649 //     h[0]=TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
 650 //     h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2];
 651 //     h[3]=d->GetTdc(); h[4]=d->GetAdc();
 652
 653 //     AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),dig,i);
 654 //     InsertCluster(cl);
 655 //   }
 656
 657 //   return 0;
 658 // }
 659 // //_________________________________________________________________________
 660 // void AliTOFtrackerMI::UnloadClusters() {
 661 //   //--------------------------------------------------------------------
 662 //   //This function unloads TOF clusters
 663 //   //--------------------------------------------------------------------
 664 //   for (Int_t i=0; i<fN; i++) delete fClusters[i];
 665 //   fN=0;
 666 // }
 667
 668
 669
 670 //_________________________________________________________________________
 671 Int_t AliTOFtrackerMI::LoadClusters(TTree *cTree) {
 672   //--------------------------------------------------------------------
 673   //This function loads the TOF clusters
 674   //--------------------------------------------------------------------
 675
 676   TBranch *branch=cTree->GetBranch("TOF");
 677   if (!branch) {
 678     AliError("can't get the branch with the TOF clusters !");
 679     return 1;
 680   }
 681
 682   TClonesArray dummy("AliTOFcluster",10000), *clusters=&dummy;
 683   branch->SetAddress(&clusters);
 684
 685   cTree->GetEvent(0);
 686   Int_t nc=clusters->GetEntriesFast();
 687   AliInfo(Form("Number of clusters: %d",nc));
 688
 689   for (Int_t i=0; i<nc; i++) {
 690     AliTOFcluster *c=(AliTOFcluster*)clusters->UncheckedAt(i);
 691     /*
 692     for (Int_t jj=0; jj<5; jj++) dig[jj]=c->GetDetInd(jj);
 693
 694     Double_t h[5];
 695     h[0]=c->GetR();
 696     h[1]=c->GetPhi();
 697     h[2]=c->GetZ();
 698     h[3]=c->GetTDC();
 699     h[4]=c->GetADC();
 700
 701     Int_t indexDig[3];
 702     for (Int_t jj=0; jj<3; jj++) indexDig[jj] = c->GetLabel(jj);
 703
 704     AliTOFcluster *cl=new AliTOFcluster(h,c->GetTracks(),dig,i);
 705     */
 706
 707     //    fClusters[i]=c; fN++;
 708     fClusters[i]=new AliTOFcluster(*c); fN++;
 709
 710     //AliInfo(Form("%4i %4i  %f %f %f  %f %f   %2i %1i %2i %1i %2i",i, fClusters[i]->GetIndex(),fClusters[i]->GetZ(),fClusters[i]->GetR(),fClusters[i]->GetPhi(), fClusters[i]->GetTDC(),fClusters[i]->GetADC(),fClusters[i]->GetDetInd(0),fClusters[i]->GetDetInd(1),fClusters[i]->GetDetInd(2),fClusters[i]->GetDetInd(3),fClusters[i]->GetDetInd(4)));
 711     //AliInfo(Form("%i %f",i, fClusters[i]->GetZ()));
 712   }
 713
 714   //AliInfo(Form("Number of clusters: %d",fN));
 715
 716   return 0;
 717 }
 718 //_________________________________________________________________________
 719 void AliTOFtrackerMI::UnloadClusters() {
 720   //--------------------------------------------------------------------
 721   //This function unloads TOF clusters
 722   //--------------------------------------------------------------------
 723   for (Int_t i=0; i<fN; i++) {
 724     delete fClusters[i];
 725     fClusters[i] = 0x0;
 726   }
 727   fN=0;
 728 }
 729
 730
 731
 732
 733 //_________________________________________________________________________
 734 Int_t AliTOFtrackerMI::InsertCluster(AliTOFcluster *c) {
 735   //--------------------------------------------------------------------
 736   //This function adds a cluster to the array of clusters sorted in Z
 737   //--------------------------------------------------------------------
 738   if (fN==kMaxCluster) {
 739     AliError("Too many clusters !");
 740     return 1;
 741   }
 742
 743   if (fN==0) {fClusters[fN++]=c; return 0;}
 744   Int_t i=FindClusterIndex(c->GetZ());
 745   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTOFcluster*));
 746   fClusters[i]=c; fN++;
 747
 748   return 0;
 749 }
 750
 751 //_________________________________________________________________________
 752 Int_t AliTOFtrackerMI::FindClusterIndex(Double_t z) const {
 753   //--------------------------------------------------------------------
 754   // This function returns the index of the nearest cluster
 755   //--------------------------------------------------------------------
 756   if (fN==0) return 0;
 757   if (z <= fClusters[0]->GetZ()) return 0;
 758   if (z > fClusters[fN-1]->GetZ()) return fN;
 759   Int_t b=0, e=fN-1, m=(b+e)/2;
 760   for (; b<e; m=(b+e)/2) {
 761     if (z > fClusters[m]->GetZ()) b=m+1;
 762     else e=m;
 763   }
 764   return m;
 765 }
 766
 767
 768
 769 Float_t AliTOFtrackerMI::GetLinearDistances(AliTOFtrack * track, AliTOFcluster *cluster, Float_t distances[5])
 770 {
 771   //
 772   // calclates distance between cluster and track
 773   // use linear aproximation
 774   //
 775   const Float_t kRaddeg = 180/3.14159265358979312;
 776   //
 777   //  Float_t tiltangle  = fGeom->GetAngles(cluster->fdetIndex[1],cluster->fdetIndex[2])/kRaddeg;  //tiltangle
 778   Int_t cind[5];
 779   cind[0]= cluster->GetDetInd(0);
 780   cind[1]= cluster->GetDetInd(1);
 781   cind[2]= cluster->GetDetInd(2);
 782   cind[3]= cluster->GetDetInd(3);
 783   cind[4]= cluster->GetDetInd(4);
 784   Float_t tiltangle  = fGeom->GetAngles(cluster->GetDetInd(1),cluster->GetDetInd(2))/kRaddeg;  //tiltangle
 785
 786   Float_t cpos[3];  //cluster position
 787   Float_t cpos0[3];  //cluster position
 788   //  fGeom->GetPos(cluster->fdetIndex,cpos);
 789   //fGeom->GetPos(cluster->fdetIndex,cpos0);
 790   //
 791   fGeom->GetPos(cind,cpos);
 792   fGeom->GetPos(cind,cpos0);
 793
 794   Float_t phi = TMath::ATan2(cpos[1],cpos[0]);
 795   if(phi<0) phi=2.*TMath::Pi()+phi;
 796   //  Get the local angle in the sector philoc
 797   Float_t phiangle   = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
 798   //
 799   Double_t v0[3];
 800   Double_t dir[3];
 801   track->GetXYZ(v0);
 802   track->GetPxPyPz(dir);
 803   dir[0]/=track->GetP();
 804   dir[1]/=track->GetP();
 805   dir[2]/=track->GetP();
 806   //
 807   //
 808   //rotate 0
 809   Float_t sinphi = TMath::Sin(phiangle);
 810   Float_t cosphi = TMath::Cos(phiangle);
 811   Float_t sinth  = TMath::Sin(tiltangle);
 812   Float_t costh  = TMath::Cos(tiltangle);
 813   //
 814   Float_t temp;
 815   temp    =  cpos[0]*cosphi+cpos[1]*sinphi;
 816   cpos[1] = -cpos[0]*sinphi+cpos[1]*cosphi;
 817   cpos[0] = temp;
 818   temp    =  v0[0]*cosphi+v0[1]*sinphi;
 819   v0[1] = -v0[0]*sinphi+v0[1]*cosphi;
 820   v0[0] = temp;
 821   //
 822   temp    =  cpos[0]*costh+cpos[2]*sinth;
 823   cpos[2] = -cpos[0]*sinth+cpos[2]*costh;
 824   cpos[0] = temp;
 825   temp    =  v0[0]*costh+v0[2]*sinth;
 826   v0[2]   = -v0[0]*sinth+v0[2]*costh;
 827   v0[0]   = temp;
 828   //
 829   //
 830   //rotate direction vector
 831   //
 832   temp    =  dir[0]*cosphi+dir[1]*sinphi;
 833   dir[1] = -dir[0]*sinphi+dir[1]*cosphi;
 834   dir[0] = temp;
 835   //
 836   temp    =  dir[0]*costh+dir[2]*sinth;
 837   dir[2]  = -dir[0]*sinth+dir[2]*costh;
 838   dir[0]  = temp;
 839   //
 840   Float_t v3[3];
 841   Float_t k = (cpos[0]-v0[0])/dir[0];
 842   v3[0] = v0[0]+k*dir[0];
 843   v3[1] = v0[1]+k*dir[1];
 844   v3[2] = v0[2]+k*dir[2];
 845   //
 846   distances[0] = v3[0]-cpos[0];
 847   distances[1] = v3[1]-cpos[1];
 848   distances[2] = v3[2]-cpos[2];
 849   distances[3] = TMath::Sqrt( distances[0]*distances[0]+distances[1]*distances[1]+distances[2]*distances[2]); //distance
 850   distances[4] = k;  //length
 851
 852   //
 853   // Debuging part of the matching
 854   //
 855   if (track->GetLabel()==cluster->GetLabel(0) ||track->GetLabel()==cluster->GetLabel(1) ){
 856     TTreeSRedirector& cstream = *fDebugStreamer;
 857     Float_t tdc = cluster->GetTDC();
 858     cstream<<"Tracks"<<
 859       "TOF.="<<track<<
 860       "Cx="<<cpos0[0]<<
 861       "Cy="<<cpos0[1]<<
 862       "Cz="<<cpos0[2]<<
 863       "Dist="<<k<<
 864       "Dist0="<<distances[0]<<
 865       "Dist1="<<distances[1]<<
 866       "Dist2="<<distances[2]<<
 867       "TDC="<<tdc<<
 868       "\n";
 869   }
 870   return distances[3];
 871 }
 872
 873
 874
 875 void    AliTOFtrackerMI::GetLikelihood(Float_t dy, Float_t dz, const Double_t *cov, AliTOFtrack * /*track*/, Float_t & py, Float_t &pz)
 876 {
 877   //
 878   //  get likelihood - track covariance taken
 879   //  75 % of gauss with expected sigma
 880   //  25 % of gauss with extended sigma
 881
 882   Double_t kMaxSigmaY  = 0.6;  // ~ 90% of TRD tracks
 883   Double_t kMaxSigmaZ  = 1.2;  // ~ 90% of TRD tracks
 884   Double_t kMeanSigmaY = 0.25; // mean TRD sigma
 885   Double_t kMeanSigmaZ = 0.5;  // mean TRD sigma
 886
 887
 888   Float_t normwidth, normd, p0,p1;
 889   Float_t sigmay = TMath::Max(TMath::Sqrt(cov[0]+kMeanSigmaY*kMeanSigmaY),kMaxSigmaY);
 890   Float_t sigmaz = TMath::Max(TMath::Sqrt(cov[2]+kMeanSigmaZ*kMeanSigmaZ),kMaxSigmaZ);
 891
 892   py=0;
 893   pz=0;
 894   //
 895   // py calculation   - 75% admixture of original sigma - 25% tails
 896   //
 897   normwidth = fDy/sigmay;
 898   normd     = dy/sigmay;
 899   p0 = 0.5*(1+TMath::Erf(normd-normwidth*0.5));
 900   p1 = 0.5*(1+TMath::Erf(normd+normwidth*0.5));
 901   py+= 0.75*(p1-p0);
 902   //
 903   normwidth = fDy/(3.*sigmay);
 904   normd     = dy/(3.*sigmay);
 905   p0 = 0.5*(1+TMath::Erf(normd-normwidth*0.5));
 906   p1 = 0.5*(1+TMath::Erf(normd+normwidth*0.5));
 907   py+= 0.25*(p1-p0);
 908   //
 909   // pz calculation   - 75% admixture of original sigma - 25% tails
 910   //
 911   normwidth = fDz/sigmaz;
 912   normd     = dz/sigmaz;
 913   p0 = 0.5*(1+TMath::Erf(normd-normwidth*0.5));
 914   p1 = 0.5*(1+TMath::Erf(normd+normwidth*0.5));
 915   pz+= 0.75*(p1-p0);
 916   //
 917   normwidth = fDz/(3.*sigmaz);
 918   normd     = dz/(3.*sigmaz);
 919   p0 = 0.5*(1+TMath::Erf(normd-normwidth*0.5));
 920   p1 = 0.5*(1+TMath::Erf(normd+normwidth*0.5));
 921   pz+= 0.25*(p1-p0);
 922 }