TOF/AliTOFtrackerV2.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 // AliTOFtrackerV2 Class                                              //
  19 // Task: Perform association of the ESD tracks to TOF Clusters        //
  20 // and Update ESD track with associated TOF Cluster parameters        //
  21 //                                                                    //
  22 // -- Authors : S. Arcelli, C. Zampolli (Bologna University and INFN) //
  23 // -- Contacts: Annalisa.De.Caro@cern.ch                              //
  24 // --         : Chiara.Zampolli@bo.infn.it                            //
  25 // --         : Silvia.Arcelli@bo.infn.it                             //
  26 //                                                                    //
  27 //--------------------------------------------------------------------//
  28
  29 #include <Rtypes.h>
  30 #include <TROOT.h>
  31
  32 #include <TClonesArray.h>
  33 #include <TObjArray.h>
  34 #include <TGeoManager.h>
  35 #include <TTree.h>
  36
  37 #include "AliGeomManager.h"
  38 #include "AliESDtrack.h"
  39 #include "AliESDEvent.h"
  40 #include "AliESDpid.h"
  41 #include "AliESDTOFcluster.h"
  42 #include "AliLog.h"
  43 #include "AliTrackPointArray.h"
  44 #include "AliCDBManager.h"
  45
  46 #include "AliTOFRecoParam.h"
  47 #include "AliTOFReconstructor.h"
  48 #include "AliTOFcluster.h"
  49 #include "AliTOFGeometry.h"
  50 #include "AliTOFtrackerV2.h"
  51 #include "AliTOFtrack.h"
  52
  53 extern TGeoManager *gGeoManager;
  54
  55 ClassImp(AliTOFtrackerV2)
  56
  57 //_____________________________________________________________________________
  58 AliTOFtrackerV2::AliTOFtrackerV2():
  59   fkRecoParam(0x0),
  60   fGeom(0x0),
  61   fN(0),
  62   fNseeds(0),
  63   fNseedsTOF(0),
  64   fnunmatch(0),
  65   fnmatch(0),
  66   fTracks(new TClonesArray("AliTOFtrack")),
  67   fSeeds(new TObjArray(100)),
  68   fClusters(0x0)
  69 {
  70   //AliTOFtrackerV2 main Ctor
  71
  72   // Getting the geometry
  73   fGeom = new AliTOFGeometry();
  74
  75 }
  76 //_____________________________________________________________________________
  77 AliTOFtrackerV2::~AliTOFtrackerV2() {
  78   //
  79   // Dtor
  80   //
  81
  82   if(!(AliCDBManager::Instance()->GetCacheFlag())){
  83     delete fkRecoParam;
  84   }
  85   delete fGeom;
  86   if (fTracks){
  87     fTracks->Delete();
  88     delete fTracks;
  89     fTracks=0x0;
  90   }
  91   if (fSeeds){
  92     fSeeds->Delete();
  93     delete fSeeds;
  94     fSeeds=0x0;
  95   }
  96
  97   if(fClusters){
  98     delete[] fClusters;
  99     fClusters = NULL;
 100   }
 101
 102 }
 103 //_____________________________________________________________________________
 104 void AliTOFtrackerV2::GetPidSettings(AliESDpid *esdPID) {
 105   //
 106   // Sets TOF resolution from RecoParams
 107   //
 108   if (fkRecoParam)
 109     esdPID->GetTOFResponse().SetTimeResolution(fkRecoParam->GetTimeResolution());
 110   else
 111     AliWarning("fkRecoParam not yet set; cannot set PID settings");
 112 }
 113 //_____________________________________________________________________________
 114 Int_t AliTOFtrackerV2::PropagateBack(AliESDEvent * const event) {
 115   //
 116   // Gets seeds from ESD event and Match with TOF Clusters
 117   //
 118
 119   //Update the matched ESD tracks
 120   // needed in case of call of TOF info before of the selection of matching and in case of no clusters available at all
 121   if(fN==0)
 122     event->SetTOFcluster(1,fClusters);
 123   else
 124     event->SetTOFcluster(fN,fClusters);
 125
 126   if (fN==0) {
 127     AliInfo("No TOF recPoints to be matched with reconstructed tracks");
 128     return 0;
 129   }
 130
 131   // initialize RecoParam for current event
 132   AliDebug(1,"Initializing params for TOF");
 133
 134   fkRecoParam = AliTOFReconstructor::GetRecoParam();  // instantiate reco param from STEER...
 135
 136   if (fkRecoParam == 0x0) {
 137     AliFatal("No Reco Param found for TOF!!!");
 138   }
 139
 140   //Initialise some counters
 141
 142   fNseeds=0;
 143   fNseedsTOF=0;
 144   fnunmatch=0;
 145   fnmatch=0;
 146
 147   Int_t ntrk=event->GetNumberOfTracks();
 148   fNseeds = ntrk;
 149
 150   //Load ESD tracks into a local Array of ESD Seeds
 151   for (Int_t i=0; i<fNseeds; i++){
 152     fSeeds->AddLast(event->GetTrack(i));
 153     event->GetTrack(i)->SetESDEvent(event);
 154   }
 155   //Prepare ESD tracks candidates for TOF Matching
 156   CollectESD();
 157
 158   if (fNseeds==0 || fNseedsTOF==0) {
 159     AliInfo("No seeds to try TOF match");
 160     return 0;
 161   }
 162
 163   // clusterize before of matching
 164   Clusterize();
 165
 166   //Second Step with Looser Matching Criterion
 167   MatchTracks();
 168
 169   AliInfo(Form("Number of matched tracks = %d",fnmatch));
 170
 171   for (Int_t i=0; i<ntrk; i++) {
 172     AliESDtrack *t=event->GetTrack(i);
 173     //t->SetESDEvent(event);
 174     AliESDtrack *seed =(AliESDtrack*)fSeeds->At(i);
 175
 176     if ( (seed->GetStatus()&AliESDtrack::kTOFin)!=0 ) {
 177       t->SetStatus(AliESDtrack::kTOFin);
 178       if ( (seed->GetStatus()&AliESDtrack::kTOFout)!=0 ) {
 179         t->SetStatus(AliESDtrack::kTOFout);
 180         //t->SetTOFclusterArray(seed->GetNTOFclusters(),seed->GetTOFclusterArray());
 181         t->SortTOFcluster();
 182
 183         // Make attention, please:
 184         //      AliESDtrack::fTOFInfo array does not be stored in the AliESDs.root file
 185         //      it is there only for a check during the reconstruction step.
 186         Float_t info[10]; seed->GetTOFInfo(info);
 187         t->SetTOFInfo(info);
 188         AliDebug(3,Form(" distance=%f; residual in the pad reference frame: dX=%f, dZ=%f", info[0],info[1],info[2]));
 189
 190         /*
 191         Double_t alphaA = (Double_t)t->GetAlpha();
 192         Double_t xA = (Double_t)t->GetX();
 193         Double_t yA = (Double_t)t->GetY();
 194         Double_t zA = (Double_t)t->GetZ();
 195         Double_t p1A = (Double_t)t->GetSnp();
 196         Double_t p2A = (Double_t)t->GetTgl();
 197         Double_t p3A = (Double_t)t->GetSigned1Pt();
 198         const Double_t *covA = (Double_t*)t->GetCovariance();
 199
 200         // Check done:
 201         //       by calling the AliESDtrack::UpdateTrackParams,
 202         //       the current track parameters are changed
 203         //       and it could cause refit problems.
 204         //       We need to update only the following track parameters:
 205         //            the track length and expected times.
 206         //       Removed AliESDtrack::UpdateTrackParams call
 207         //       Called AliESDtrack::SetIntegratedTimes(...) and
 208         //       AliESDtrack::SetIntegratedLength() routines.
 209
 210         AliTOFtrack *track = new AliTOFtrack(*seed);
 211         t->UpdateTrackParams(track,AliESDtrack::kTOFout); // to be checked - AdC
 212         delete track;
 213         Double_t time[AliPID::kSPECIESC]; t->GetIntegratedTimes(time);
 214         Double_t alphaB = (Double_t)t->GetAlpha();
 215         Double_t xB = (Double_t)t->GetX();
 216         Double_t yB = (Double_t)t->GetY();
 217         Double_t zB = (Double_t)t->GetZ();
 218         Double_t p1B = (Double_t)t->GetSnp();
 219         Double_t p2B = (Double_t)t->GetTgl();
 220         Double_t p3B = (Double_t)t->GetSigned1Pt();
 221         const Double_t *covB = (Double_t*)t->GetCovariance();
 222         AliDebug(2,"Track params -now(before)-:");
 223         AliDebug(2,Form("    X: %f(%f), Y: %f(%f), Z: %f(%f) --- alpha: %f(%f)",
 224                         xB,xA,
 225                         yB,yA,
 226                         zB,zA,
 227                         alphaB,alphaA));
 228         AliDebug(2,Form("    p1: %f(%f), p2: %f(%f), p3: %f(%f)",
 229                         p1B,p1A,
 230                         p2B,p2A,
 231                         p3B,p3A));
 232         AliDebug(2,Form("    cov1: %f(%f), cov2: %f(%f), cov3: %f(%f)"
 233                         " cov4: %f(%f), cov5: %f(%f), cov6: %f(%f)"
 234                         " cov7: %f(%f), cov8: %f(%f), cov9: %f(%f)"
 235                         " cov10: %f(%f), cov11: %f(%f), cov12: %f(%f)"
 236                         " cov13: %f(%f), cov14: %f(%f), cov15: %f(%f)",
 237                         covB[0],covA[0],
 238                         covB[1],covA[1],
 239                         covB[2],covA[2],
 240                         covB[3],covA[3],
 241                         covB[4],covA[4],
 242                         covB[5],covA[5],
 243                         covB[6],covA[6],
 244                         covB[7],covA[7],
 245                         covB[8],covA[8],
 246                         covB[9],covA[9],
 247                         covB[10],covA[10],
 248                         covB[11],covA[11],
 249                         covB[12],covA[12],
 250                         covB[13],covA[13],
 251                         covB[14],covA[14]
 252                         ));
 253         */
 254         Double_t time[AliPID::kSPECIESC]; t->GetIntegratedTimes(time,AliPID::kSPECIESC);
 255         AliDebug(2,Form(" TOF params: %6d  %f %f %f %f %f %6d %3d %f",
 256                         i,
 257                         t->GetTOFsignalRaw(),t->GetTOFsignal(),t->GetTOFsignalToT(),
 258                         t->GetTOFsignalDz(),t->GetTOFsignalDx(),t->GetTOFCalChannel(),
 259                         t->GetTOFcluster(),t->GetIntegratedLength()));
 260         AliDebug(2,Form("  %f %f %f %f %f %f %f %f %f",
 261                         time[0], time[1], time[2], time[3], time[4], time[5], time[6], time[7], time[8]));
 262       }
 263     }
 264   }
 265
 266   fSeeds->Clear();
 267   fTracks->Clear();
 268
 269   AliInfo(Form("Number of cluster to be checked = %d",fN));
 270   if(fN){
 271     Int_t *matchmap = new Int_t[fN];
 272     event->SetTOFcluster(fN,fClusters,matchmap);
 273     for (Int_t i=0; i<ntrk; i++) { // remapping after TOF matching selection
 274       AliESDtrack *t=event->GetTrack(i);
 275       t->ReMapTOFcluster(fN,matchmap);
 276     }
 277
 278     delete[] matchmap;
 279   }
 280
 281
 282
 283   return 0;
 284
 285 }
 286 //_________________________________________________________________________
 287 void AliTOFtrackerV2::CollectESD() {
 288    //prepare the set of ESD tracks to be matched to clusters in TOF
 289
 290   Int_t seedsTOF1=0;
 291   Int_t seedsTOF3=0;
 292   Int_t seedsTOF2=0;
 293
 294   TClonesArray &aTOFTrack = *fTracks;
 295   for (Int_t i=0; i<fNseeds; i++) {
 296
 297     AliESDtrack *t =(AliESDtrack*)fSeeds->At(i);
 298     if ((t->GetStatus()&AliESDtrack::kTPCout)==0)continue;
 299
 300     AliTOFtrack *track = new AliTOFtrack(*t); // New
 301     Float_t x = (Float_t)track->GetX(); //New
 302
 303     // TRD 'good' tracks
 304     if ( ( (t->GetStatus()&AliESDtrack::kTRDout)!=0 ) ) {
 305
 306       AliDebug(1,Form(" Before propagation till inner TOF radius, ESDtrackLength=%f, TOFtrackLength=%f",t->GetIntegratedLength(),track->GetIntegratedLength()));
 307
 308       // TRD 'good' tracks, already propagated at 371 cm
 309       if ( x >= AliTOFGeometry::Rmin() ) {
 310
 311         if ( track->PropagateToInnerTOF() ) {
 312
 313           AliDebug(1,Form(" TRD propagated track till rho = %fcm."
 314                           " And then the track has been propagated till rho = %fcm.",
 315                           x, (Float_t)track->GetX()));
 316
 317           track->SetSeedIndex(i);
 318           t->UpdateTrackParams(track,AliESDtrack::kTOFin);
 319           new(aTOFTrack[fNseedsTOF]) AliTOFtrack(*track);
 320           fNseedsTOF++;
 321           seedsTOF1++;
 322
 323           AliDebug(1,Form(" After propagation till inner TOF radius, ESDtrackLength=%f, TOFtrackLength=%f",t->GetIntegratedLength(),track->GetIntegratedLength()));
 324         }
 325         delete track;
 326
 327       }
 328       else { // TRD 'good' tracks, propagated rho<371cm
 329
 330         if  ( track->PropagateToInnerTOF() ) {
 331
 332           AliDebug(1,Form(" TRD propagated track till rho = %fcm."
 333                           " And then the track has been propagated till rho = %fcm.",
 334                           x, (Float_t)track->GetX()));
 335
 336           track->SetSeedIndex(i);
 337           t->UpdateTrackParams(track,AliESDtrack::kTOFin);
 338           new(aTOFTrack[fNseedsTOF]) AliTOFtrack(*track);
 339           fNseedsTOF++;
 340           seedsTOF3++;
 341
 342           AliDebug(1,Form(" After propagation till inner TOF radius, ESDtrackLength=%f, TOFtrackLength=%f",t->GetIntegratedLength(),track->GetIntegratedLength()));
 343         }
 344         delete track;
 345
 346       }
 347       //delete track;
 348     }
 349
 350     else { // Propagate the rest of TPCbp
 351
 352       AliDebug(1,Form(" Before propagation till inner TOF radius, ESDtrackLength=%f, TOFtrackLength=%f",t->GetIntegratedLength(),track->GetIntegratedLength()));
 353
 354       if ( track->PropagateToInnerTOF() ) {
 355
 356         AliDebug(1,Form(" TPC propagated track till rho = %fcm."
 357                         " And then the track has been propagated till rho = %fcm.",
 358                         x, (Float_t)track->GetX()));
 359
 360         track->SetSeedIndex(i);
 361         t->UpdateTrackParams(track,AliESDtrack::kTOFin);
 362         new(aTOFTrack[fNseedsTOF]) AliTOFtrack(*track);
 363         fNseedsTOF++;
 364         seedsTOF2++;
 365
 366         AliDebug(1,Form(" After propagation till inner TOF radius, ESDtrackLength=%f, TOFtrackLength=%f",t->GetIntegratedLength(),track->GetIntegratedLength()));
 367       }
 368       delete track;
 369     }
 370   }
 371
 372   AliInfo(Form("Number of TOF seeds = %d (kTRDout371 = %d, kTRDoutLess371 = %d, !kTRDout = %d)",fNseedsTOF,seedsTOF1,seedsTOF3,seedsTOF2));
 373
 374   // Sort according uncertainties on track position
 375   fTracks->Sort();
 376
 377 }
 378
 379 //_________________________________________________________________________
 380 void AliTOFtrackerV2::MatchTracks() {
 381   //
 382   //Match ESD tracks to clusters in TOF
 383   //
 384
 385   // Parameters used/regulating the reconstruction
 386   static Float_t detDepth=18.;
 387   static Float_t padDepth=0.5;
 388
 389   const Float_t kSpeedOfLight= 2.99792458e-2; // speed of light [cm/ps]
 390
 391   Float_t dY=AliTOFGeometry::XPad();
 392   Float_t dZ=AliTOFGeometry::ZPad();
 393
 394   Float_t sensRadius = fkRecoParam->GetSensRadius();
 395   Float_t stepSize   = fkRecoParam->GetStepSize();
 396   Float_t scaleFact  = fkRecoParam->GetWindowScaleFact();
 397   Float_t dyMax=fkRecoParam->GetWindowSizeMaxY();
 398   Float_t dzMax=fkRecoParam->GetWindowSizeMaxZ();
 399   Float_t dCut=10.;//fkRecoParam->GetDistanceCut(); // This is to be loaded by OCDB. It should be 10cm always.
 400   Double_t maxChi2=fkRecoParam->GetMaxChi2TRD();
 401   Bool_t timeWalkCorr = fkRecoParam->GetTimeWalkCorr();
 402   AliDebug(1,"++++++++++++++TOF Reconstruction Parameters:++++++++++++++");
 403   AliDebug(1,Form("TOF sens radius: %f",sensRadius));
 404   AliDebug(1,Form("TOF step size: %f",stepSize));
 405   AliDebug(1,Form("TOF Window scale factor: %f",scaleFact));
 406   AliDebug(1,Form("TOF Window max dy: %f",dyMax));
 407   AliDebug(1,Form("TOF Window max dz: %f",dzMax));
 408   AliDebug(1,Form("TOF distance Cut: %f",dCut));
 409   AliDebug(1,Form("TOF Max Chi2: %f",maxChi2));
 410   AliDebug(1,Form("Time Walk Correction? : %d",timeWalkCorr));
 411
 412   //Match ESD tracks to clusters in TOF
 413
 414   // Get the number of propagation steps
 415   Int_t nSteps=(Int_t)(detDepth/stepSize);
 416   AliDebug(1,Form(" Number of steps to be done %d",nSteps));
 417
 418   AliDebug(1,"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
 419
 420   //PH Arrays (moved outside of the loop)
 421   Float_t * trackPos[4];
 422   for (Int_t ii=0; ii<4; ii++) trackPos[ii] = new Float_t[nSteps];
 423   Int_t * clind = new Int_t[fN];
 424
 425   // Some init
 426   const Int_t kNclusterMax = 1000; // related to fN value
 427   TGeoHMatrix global[kNclusterMax];
 428
 429   //The matching loop
 430   for (Int_t iseed=0; iseed<fNseedsTOF; iseed++) {
 431
 432     for (Int_t ii=0; ii<fN; ii++) clind[ii]=-1;
 433     for (Int_t ii=0; ii<kNclusterMax; ii++) global[ii] = 0x0;
 434     for (Int_t ii=0; ii<4; ii++)
 435       for (Int_t jj=0; jj<nSteps; jj++) trackPos[ii][jj]=0.;
 436
 437     AliTOFtrack *track =(AliTOFtrack*)fTracks->UncheckedAt(iseed);
 438     AliESDtrack *t =(AliESDtrack*)fSeeds->At(track->GetSeedIndex());
 439     AliTOFtrack *trackTOFin = new AliTOFtrack(*track);
 440
 441     Double_t timesOr[AliPID::kSPECIESC]; t->GetIntegratedTimes(timesOr,AliPID::kSPECIESC); // in ps
 442
 443     // Determine a window around the track
 444     Double_t x,par[5];
 445     trackTOFin->GetExternalParameters(x,par);
 446     Double_t cov[15];
 447     trackTOFin->GetExternalCovariance(cov);
 448
 449     if (cov[0]<0. || cov[2]<0.) {
 450       AliWarning(Form("Very strange track (%d)! At least one of its covariance matrix diagonal elements is negative!",iseed));
 451       delete trackTOFin;
 452       continue;
 453     }
 454
 455     Double_t dphi=
 456       scaleFact*
 457       ((5*TMath::Sqrt(TMath::Abs(cov[0])) + 0.5*dY + 2.5*TMath::Abs(par[2]))/sensRadius);
 458     Double_t dz=
 459        scaleFact*
 460        (5*TMath::Sqrt(TMath::Abs(cov[2])) + 0.5*dZ + 2.5*TMath::Abs(par[3]));
 461
 462     Double_t phi=TMath::ATan2(par[0],x) + trackTOFin->GetAlpha();
 463     if (phi<-TMath::Pi())phi+=2*TMath::Pi();
 464     if (phi>=TMath::Pi())phi-=2*TMath::Pi();
 465     Double_t z=par[1];
 466
 467     //upper limit on window's size.
 468     if (dz> dzMax) dz=dzMax;
 469     if (dphi*sensRadius> dyMax) dphi=dyMax/sensRadius;
 470
 471
 472     // find the clusters in the window of the track
 473     Int_t nc=0;
 474     for (Int_t k=FindClusterIndex(z-dz); k<fN; k++) {
 475
 476       if (nc>=kNclusterMax) {
 477         AliWarning("No more matchable clusters can be stored! Please, increase the corresponding vectors size.");
 478         break;
 479       }
 480
 481       AliESDTOFcluster *c=&(fClusters[k]);
 482       if (c->GetZ() > z+dz) break;
 483       if (!c->GetStatus()) {
 484         AliDebug(1,"Cluster in channel declared bad!");
 485         continue; // skip bad channels as declared in OCDB
 486       }
 487
 488       Double_t dph=TMath::Abs(c->GetPhi()-phi);
 489       if (dph>TMath::Pi()) dph-=2.*TMath::Pi();
 490       if (TMath::Abs(dph)>dphi) continue;
 491
 492       Double_t yc=(c->GetPhi() - trackTOFin->GetAlpha())*c->GetR();
 493       Double_t p[2]={yc, c->GetZ()};
 494       Double_t cov2[3]= {dY*dY/12., 0., dZ*dZ/12.};
 495       if (trackTOFin->AliExternalTrackParam::GetPredictedChi2(p,cov2) > maxChi2)continue;
 496
 497       clind[nc] = k;
 498       Char_t path[200];
 499       Int_t ind[5]; fGeom->GetVolumeIndices(c->GetTOFchannel(),ind);
 500       fGeom->GetVolumePath(ind,path);
 501       gGeoManager->cd(path);
 502       global[nc] = *gGeoManager->GetCurrentMatrix();
 503       nc++;
 504     }
 505
 506     if (nc == 0 ) {
 507       AliDebug(1,Form("No available clusters for the track number %d",iseed));
 508       fnunmatch++;
 509       delete trackTOFin;
 510       continue;
 511     }
 512
 513     AliDebug(1,Form(" Number of available TOF clusters for the track number %d: %d",iseed,nc));
 514
 515     //start fine propagation
 516
 517     Double_t *times[AliPID::kSPECIESC];
 518     for(Int_t isp=0;isp < AliPID::kSPECIESC;isp++){
 519       times[isp] = new Double_t[nSteps];
 520     }
 521
 522     Int_t nStepsDone = 0;
 523     for( Int_t istep=0; istep<nSteps; istep++){
 524
 525       // First of all, propagate the track...
 526       Float_t xs = AliTOFGeometry::RinTOF()+istep*stepSize;
 527       if (!(trackTOFin->PropagateTo(xs))) break;
 528
 529       //  ...and then, if necessary, rotate the track
 530       Double_t ymax = xs*TMath::Tan(0.5*AliTOFGeometry::GetAlpha());
 531       Double_t ysect = trackTOFin->GetY();
 532       if (ysect > ymax) {
 533         if (!(trackTOFin->Rotate(AliTOFGeometry::GetAlpha()))) break;
 534       } else if (ysect <-ymax) {
 535         if (!(trackTOFin->Rotate(-AliTOFGeometry::GetAlpha()))) break;
 536       }
 537
 538       Double_t mom = trackTOFin->P();
 539
 540       if(istep == 0){
 541         for(Int_t isp=0;isp<AliPID::kSPECIESC;isp++){
 542           Double_t mass=AliPID::ParticleMass(isp);
 543           Double_t momz = mom*AliPID::ParticleCharge(isp);
 544           times[isp][nStepsDone] = stepSize/kSpeedOfLight*TMath::Sqrt(momz*momz+mass*mass)/momz;
 545         }
 546       }
 547       else{
 548         for(Int_t isp=0;isp<AliPID::kSPECIESC;isp++){
 549           Double_t mass=AliPID::ParticleMass(isp);
 550           Double_t momz = mom*AliPID::ParticleCharge(isp);
 551           times[isp][nStepsDone] = times[isp][nStepsDone-1] + (trackTOFin->GetIntegratedLength()-trackPos[3][nStepsDone-1])/kSpeedOfLight*TMath::Sqrt(momz*momz+mass*mass)/momz;
 552         }
 553       }
 554
 555       // store the running point (Globalrf) - fine propagation
 556
 557       Double_t r[3]; trackTOFin->GetXYZ(r);
 558       trackPos[0][nStepsDone]= (Float_t) r[0];
 559       trackPos[1][nStepsDone]= (Float_t) r[1];
 560       trackPos[2][nStepsDone]= (Float_t) r[2];
 561       trackPos[3][nStepsDone]= trackTOFin->GetIntegratedLength();
 562
 563       nStepsDone++;
 564       AliDebug(3,Form(" current step %d (%d) - nStepsDone=%d",istep,nSteps,nStepsDone));
 565     }
 566
 567     if ( nStepsDone == 0 ) {
 568       AliDebug(1,Form(" No track points for track number %d",iseed));
 569       fnunmatch++;
 570       delete trackTOFin;
 571       continue;
 572     }
 573
 574     AliDebug(3,Form(" Number of steps done for the track number %d: %d",iseed,nStepsDone));
 575
 576     Int_t *isClusterMatchable = NULL;
 577     if(nc){
 578       isClusterMatchable = new Int_t[nc];
 579       for (Int_t i=0; i<nc; i++) isClusterMatchable[i] = kFALSE;
 580     }
 581
 582     Int_t nfound = 0;
 583     Bool_t accept = kFALSE;
 584     Bool_t isInside = kFALSE;
 585     for (Int_t istep=0; istep<nStepsDone; istep++) {
 586
 587       Bool_t gotInsideCluster = kFALSE;
 588       Int_t trackInsideCluster = -1;
 589
 590       Float_t ctrackPos[3];
 591       ctrackPos[0] = trackPos[0][istep];
 592       ctrackPos[1] = trackPos[1][istep];
 593       ctrackPos[2] = trackPos[2][istep];
 594
 595       //now see whether the track matches any of the TOF clusters
 596
 597       Float_t dist3d[3]={0.,0.,0.};
 598       accept = kFALSE;
 599
 600       for (Int_t i=0; i<nc; i++) {
 601
 602         // ***** NEW *****
 603         /* check whether track was inside another cluster
 604          * and in case inhibit this cluster.
 605          * this will allow to only go on and add track points for
 606          * that cluster where the track got inside first */
 607         if (gotInsideCluster && trackInsideCluster != i) {
 608           AliDebug(3,Form(" A - istep=%d ~ %d %d ~ nfound=%d",istep,trackInsideCluster,i,nfound));
 609           continue;
 610         }
 611         AliDebug(3,Form(" B - istep=%d ~ %d %d ~ nfound=%d",istep,trackInsideCluster,i,nfound));
 612
 613         /* check whether track is inside this cluster */
 614         for (Int_t hh=0; hh<3; hh++) dist3d[hh]=0.;
 615         isInside = fGeom->IsInsideThePad((TGeoHMatrix*)(&global[i]),ctrackPos,dist3d);
 616
 617         // ***** NEW *****
 618         /* if track is inside this cluster set flags which will then
 619          * inhibit to add track points for the other clusters */
 620         if (isInside) {
 621           gotInsideCluster = kTRUE;
 622           trackInsideCluster = i;
 623         }
 624
 625         Float_t yLoc = dist3d[1];
 626         Float_t rLoc = TMath::Sqrt(dist3d[0]*dist3d[0]+dist3d[2]*dist3d[2]);
 627         accept = (TMath::Abs(yLoc)<padDepth*0.5 && rLoc<dCut);
 628
 629         //***** NEW *****
 630         /* add point everytime that:
 631          * - the track is inside the cluster
 632          * - the track got inside the cluster, even when it eventually exited the cluster
 633          * - the tracks is within dCut from the cluster
 634          */
 635         if (accept || isInside || gotInsideCluster) {
 636
 637           Double_t timesCurrent[AliPID::kSPECIESC];
 638           AliDebug(3,Form(" Momentum for track %d -> %f", iseed,t->P()));
 639           for (Int_t j=0;j<AliPID::kSPECIESC;j++) {
 640             timesCurrent[j] = timesOr[j] + times[j][istep];
 641           }
 642
 643           if (TMath::Abs(dist3d[1])<stepSize && !isClusterMatchable[i]) {
 644             isClusterMatchable[i] = kTRUE;
 645             fClusters[clind[i]].Update(t->GetID(),dist3d[1],dist3d[0],dist3d[2],trackPos[3][istep],timesCurrent);//x,y,z -> tracking RF
 646             t->AddTOFcluster(clind[i]);
 647             t->SetStatus(AliESDtrack::kTOFout);
 648           }
 649           AliDebug(2,Form(" dist3dLoc[0] = %f, dist3dLoc[1] = %f, dist3dLoc[2] = %f ",dist3d[0],dist3d[1],dist3d[2]));
 650
 651           nfound++;
 652
 653           AliDebug(3,Form(" C - istep=%d ~ %d %d ~ nfound=%d",istep,trackInsideCluster,i,nfound));
 654
 655           // ***** NEW *****
 656         }//end if accept
 657
 658       } //end for on the clusters
 659     } //end for on the steps
 660     if(nc) delete[] isClusterMatchable;
 661
 662     for(Int_t isp=0;isp < AliPID::kSPECIESC;isp++){
 663       delete[] times[isp];
 664     }
 665
 666
 667     if (nfound == 0 ) {
 668       AliDebug(1,Form(" No matchable track points for the track number %d",iseed));
 669       fnunmatch++;
 670       delete trackTOFin;
 671       continue;
 672     }
 673
 674     AliDebug(1,Form(" Number of track points for the track number %d: %d",iseed,nfound));
 675
 676     Int_t nMatchedClusters = t->GetNTOFclusters();
 677
 678     if (nMatchedClusters==0) {
 679       AliDebug(1,Form("Reconstructed track %d doesn't match any TOF cluster", iseed));
 680       fnunmatch++;
 681       delete trackTOFin;
 682       continue;
 683     }
 684
 685     AliDebug(1,Form(" %d - matched (%d)",track->GetSeedIndex()/*iseed*/,nMatchedClusters));
 686
 687     fnmatch++;
 688
 689     /*
 690     AliTOFcluster cTOF = AliTOFcluster(volIdClus,
 691     (Float_t)posClus[0],(Float_t)posClus[1],(Float_t)posClus[2],
 692     (Float_t)covClus[0],(Float_t)covClus[1],(Float_t)covClus[2],
 693     (Float_t)covClus[3],(Float_t)covClus[4],(Float_t)covClus[5],
 694     tofLabels,volIndices,parClu,kTRUE,index[i]);
 695
 696     // Fill the track residual histograms.
 697     FillResiduals(trackTOFin,c,kFALSE);
 698     */
 699
 700     delete trackTOFin;
 701
 702   } // loop on fSeeds
 703
 704   for (Int_t ii=0; ii<4; ii++) delete [] trackPos[ii];
 705   delete [] clind;
 706
 707 }
 708 //_________________________________________________________________________
 709 Int_t AliTOFtrackerV2::LoadClusters(TTree *cTree) {
 710   //--------------------------------------------------------------------
 711   //This function loads the TOF clusters
 712   //--------------------------------------------------------------------
 713
 714   TBranch *branch=cTree->GetBranch("TOF");
 715   if (!branch) {
 716     AliError("can't get the branch with the TOF clusters !");
 717     return 1;
 718   }
 719
 720   static TClonesArray dummy("AliTOFcluster",10000);
 721   dummy.Clear();
 722   TClonesArray *clusters=&dummy;
 723   branch->SetAddress(&clusters);
 724
 725   cTree->GetEvent(0);
 726   fN=clusters->GetEntriesFast();
 727   AliInfo(Form("Number of clusters: %d",fN));
 728
 729   if(fClusters){
 730     delete[] fClusters;
 731     fClusters = NULL;
 732   }
 733
 734   if(fN)
 735     fClusters = new AliESDTOFcluster[fN];
 736   else{
 737     fClusters = new AliESDTOFcluster[1];
 738     fN = 1;
 739     return 0;
 740   }
 741
 742   for (Int_t i=0; i<fN; i++) {
 743     AliTOFcluster *c=(AliTOFcluster*)clusters->UncheckedAt(i);
 744     Int_t ind[5];
 745     ind[0]=c->GetDetInd(0);
 746     ind[1]=c->GetDetInd(1);
 747     ind[2]=c->GetDetInd(2);
 748     ind[3]=c->GetDetInd(3);
 749     ind[4]=c->GetDetInd(4);
 750     Int_t calindex = AliTOFGeometry::GetIndex(ind);
 751     Int_t tofLabels[3]={c->GetLabel(0),c->GetLabel(1),c->GetLabel(2)};
 752     AliESDTOFcluster esdTOFclus(i,calindex,
 753                                 AliTOFGeometry::TdcBinWidth()*c->GetTDC()/*ps*/,
 754                                 AliTOFGeometry::TdcBinWidth()*c->GetTDCRAW()/*ps*/,
 755                                 AliTOFGeometry::ToTBinWidth()*c->GetToT()*1E-3/*ns*/,
 756                                 tofLabels,
 757                                 c->GetDeltaBC(),c->GetL0L1Latency(),
 758                                 c->GetStatus(),c->GetZ(),c->GetPhi(),c->GetR());
 759
 760     fClusters[i] = esdTOFclus;
 761
 762   }
 763
 764   return 0;
 765 }
 766 //_________________________________________________________________________
 767 void AliTOFtrackerV2::UnloadClusters() {
 768   //--------------------------------------------------------------------
 769   //This function unloads TOF clusters
 770   //--------------------------------------------------------------------
 771
 772   // don't delete TOF clusters here because they should be written
 773 }
 774
 775 //_________________________________________________________________________
 776 Int_t AliTOFtrackerV2::FindClusterIndex(Double_t z) const {
 777   //--------------------------------------------------------------------
 778   // This function returns the index of the nearest cluster
 779   //--------------------------------------------------------------------
 780   if (fN==0) return 0;
 781   if (z <= fClusters[0].GetZ()) return 0;
 782   if (z > fClusters[fN-1].GetZ()) return fN;
 783   Int_t b=0, e=fN-1, m=(b+e)/2;
 784   for (; b<e; m=(b+e)/2) {
 785     if (z > fClusters[m].GetZ()) b=m+1;
 786
 787     else e=m;
 788   }
 789   return m;
 790 }
 791
 792 //_________________________________________________________________________
 793 Bool_t AliTOFtrackerV2::GetTrackPoint(Int_t index, AliTrackPoint& p) const
 794 {
 795   // Get track space point with index i
 796   // Coordinates are in the global system
 797   AliESDTOFcluster *cl = &(fClusters[index]);
 798   Float_t xyz[3];
 799   xyz[0] = cl->GetR()*TMath::Cos(cl->GetPhi());
 800   xyz[1] = cl->GetR()*TMath::Sin(cl->GetPhi());
 801   xyz[2] = cl->GetZ();
 802   Float_t phiangle = (Int_t(cl->GetPhi()*TMath::RadToDeg()/20.)+0.5)*20.*TMath::DegToRad();
 803   Float_t sinphi = TMath::Sin(phiangle), cosphi = TMath::Cos(phiangle);
 804   Int_t tofChannel=cl->GetTOFchannel();
 805   Int_t ind[5]; fGeom->GetVolumeIndices(tofChannel,ind);
 806   Float_t tiltangle = AliTOFGeometry::GetAngles(ind[1],ind[2])*TMath::DegToRad();
 807   Float_t sinth = TMath::Sin(tiltangle), costh = TMath::Cos(tiltangle);
 808   Float_t sigmay2 = AliTOFGeometry::XPad()*AliTOFGeometry::XPad()/12.;
 809   Float_t sigmaz2 = AliTOFGeometry::ZPad()*AliTOFGeometry::ZPad()/12.;
 810   Float_t cov[6];
 811   cov[0] = sinphi*sinphi*sigmay2 + cosphi*cosphi*sinth*sinth*sigmaz2;
 812   cov[1] = -sinphi*cosphi*sigmay2 + sinphi*cosphi*sinth*sinth*sigmaz2;
 813   cov[2] = -cosphi*sinth*costh*sigmaz2;
 814   cov[3] = cosphi*cosphi*sigmay2 + sinphi*sinphi*sinth*sinth*sigmaz2;
 815   cov[4] = -sinphi*sinth*costh*sigmaz2;
 816   cov[5] = costh*costh*sigmaz2;
 817   p.SetXYZ(xyz[0],xyz[1],xyz[2],cov);
 818
 819   // Detector numbering scheme
 820   Int_t nSector = AliTOFGeometry::NSectors();
 821   Int_t nPlate  = AliTOFGeometry::NPlates();
 822   Int_t nStripA = AliTOFGeometry::NStripA();
 823   Int_t nStripB = AliTOFGeometry::NStripB();
 824   Int_t nStripC = AliTOFGeometry::NStripC();
 825
 826   Int_t isector = ind[0];//cl->GetDetInd(0);
 827   if (isector >= nSector)
 828     AliError(Form("Wrong sector number in TOF (%d) !",isector));
 829   Int_t iplate = ind[1];//cl->GetDetInd(1);
 830   if (iplate >= nPlate)
 831     AliError(Form("Wrong plate number in TOF (%d) !",iplate));
 832   Int_t istrip = ind[2];//cl->GetDetInd(2);
 833
 834   Int_t stripOffset = 0;
 835   switch (iplate) {
 836   case 0:
 837     stripOffset = 0;
 838     break;
 839   case 1:
 840     stripOffset = nStripC;
 841     break;
 842   case 2:
 843     stripOffset = nStripC+nStripB;
 844     break;
 845   case 3:
 846     stripOffset = nStripC+nStripB+nStripA;
 847     break;
 848   case 4:
 849     stripOffset = nStripC+nStripB+nStripA+nStripB;
 850     break;
 851   default:
 852     AliError(Form("Wrong plate number in TOF (%d) !",iplate));
 853     break;
 854   };
 855
 856   Int_t idet = (2*(nStripC+nStripB)+nStripA)*isector +
 857                stripOffset +
 858                istrip;
 859   UShort_t volid = AliGeomManager::LayerToVolUID(AliGeomManager::kTOF,idet);
 860   p.SetVolumeID((UShort_t)volid);
 861   return kTRUE;
 862 }
 863 //_________________________________________________________________________
 864
 865 void AliTOFtrackerV2::FillClusterArray(TObjArray* arr) const
 866 {
 867   //
 868   // Returns the TOF cluster array
 869   //
 870
 871   if (fN==0)
 872     arr = 0x0;
 873   else
 874     for (Int_t i=0; i<fN; ++i) arr->Add(&(fClusters[i]));
 875
 876 }
 877 //_________________________________________________________________________
 878 Float_t AliTOFtrackerV2::CorrectTimeWalk( Float_t dist, Float_t tof) const {
 879
 880   //dummy, for the moment
 881   Float_t tofcorr=0.;
 882   if(dist<AliTOFGeometry::ZPad()*0.5){
 883     tofcorr=tof;
 884     //place here the actual correction
 885   }else{
 886     tofcorr=tof;
 887   }
 888   return tofcorr;
 889 }
 890 //_________________________________________________________________________
 891 void AliTOFtrackerV2::Clusterize(){
 892   Int_t detId[5];
 893   for(Int_t i=0; i < fN-1;i++){
 894     AliESDTOFcluster *c1=&(fClusters[i]);
 895     if(!c1->GetStatus()) continue;
 896
 897     Int_t chan1 = c1->GetTOFchannel();
 898     AliTOFGeometry::GetVolumeIndices(chan1, detId); // Get volume index from channel index
 899
 900     Int_t ieta = detId[2]/*strip*/*2 + detId[3]/*pad Z*/;
 901     if(detId[1]/*module*/ == 0) ieta += 0;
 902     else if(detId[1] == 1) ieta += 38;
 903     else if(detId[1] == 2) ieta += 76;
 904     else if(detId[1] == 3) ieta += 106;
 905     else if(detId[1] == 4) ieta += 144;
 906     Int_t iphi = detId[0]/*phi sector*/*48 + detId[4]/*pad x*/;
 907
 908
 909     for(Int_t j=i+1; j < fN;j++){
 910       AliESDTOFcluster *c2=&(fClusters[j]);
 911       if(!c2->GetStatus()) continue;
 912
 913       Int_t chan2 = c2->GetTOFchannel();
 914
 915       // check if the two TOF hits are in the same strip
 916       if(chan1/96 != chan2/96) continue;
 917
 918       AliTOFGeometry::GetVolumeIndices(chan2, detId); // Get volume index from channel index
 919       Int_t ieta2 = detId[2]/*strip*/*2 + detId[3]/*pad Z*/;
 920       if(detId[1]/*module*/ == 0) ieta2 += 0;
 921       else if(detId[1] == 1) ieta2 += 38;
 922       else if(detId[1] == 2) ieta2 += 76;
 923       else if(detId[1] == 3) ieta2 += 106;
 924       else if(detId[1] == 4) ieta2 += 144;
 925       Int_t iphi2 = detId[0]/*phi sector*/*48 + detId[4]/*pad x*/;
 926
 927       // check if the fired pad are close in space
 928       if(TMath::Abs(iphi-iphi2)>1 || TMath::Abs(ieta-ieta2)>1) continue;
 929
 930       // check if the TOF time are close enough to be merged
 931       if(TMath::Abs(c1->GetTime() - c2->GetTime()) > 500/*in ps*/) continue;
 932
 933       // merge them
 934       Int_t label[3] = {c2->GetLabel(0),c2->GetLabel(1),c2->GetLabel(2)};
 935       fClusters[i].AddTOFhit(c2->GetClusterIndex(),chan2,c2->GetTime(),c2->GetTimeRaw(),c2->GetTOT(),label,
 936                                        c2->GetDeltaBC(),c2->GetL0L1Latency(),1,c2->GetZ(),c2->GetPhi(),c2->GetR());
 937
 938       c2->SetStatus(0); // only the merged one should be used
 939       j = fN; // cluster "i" merged go to the next one ("i+1")
 940     }
 941   }
 942
 943   // second step of clusterization
 944   for(Int_t i=0; i < fN-1;i++){
 945     AliESDTOFcluster *c1=&(fClusters[i]);
 946     if(!c1->GetStatus()) continue;
 947
 948     Int_t chan1 = c1->GetTOFchannel(0);
 949     AliTOFGeometry::GetVolumeIndices(chan1, detId); // Get volume index from channel index
 950
 951     Int_t ieta = detId[2]/*strip*/*2 + detId[3]/*pad Z*/;
 952     if(detId[1]/*module*/ == 0) ieta += 0;
 953     else if(detId[1] == 1) ieta += 38;
 954     else if(detId[1] == 2) ieta += 76;
 955     else if(detId[1] == 3) ieta += 106;
 956     else if(detId[1] == 4) ieta += 144;
 957     Int_t iphi = detId[0]/*phi sector*/*48 + detId[4]/*pad x*/;
 958
 959     Int_t ieta2,iphi2,chan2=chan1;
 960     if(c1->GetNTOFhits() > 1){
 961       chan2 = c1->GetTOFchannel(1);
 962       AliTOFGeometry::GetVolumeIndices(chan2, detId); // Get volume index from channel index
 963
 964       ieta2 = detId[2]/*strip*/*2 + detId[3]/*pad Z*/;
 965       if(detId[1]/*module*/ == 0) ieta2 += 0;
 966       else if(detId[1] == 1) ieta2 += 38;
 967       else if(detId[1] == 2) ieta2 += 76;
 968       else if(detId[1] == 3) ieta2 += 106;
 969       else if(detId[1] == 4) ieta2 += 144;
 970       iphi2 = detId[0]/*phi sector*/*48 + detId[4]/*pad x*/;
 971     }
 972     else{
 973       iphi2=iphi;
 974       ieta2=ieta;
 975     }
 976
 977     for(Int_t j=i+1; j < i;j++){
 978       AliESDTOFcluster *c2=&(fClusters[j]);
 979       if(!c2->GetStatus()) continue;
 980
 981       Int_t chan3 = c2->GetTOFchannel();
 982
 983       // check if the two TOF hits are in the same strip
 984       if(chan1/96 != chan3/96) continue;
 985
 986       AliTOFGeometry::GetVolumeIndices(chan3, detId); // Get volume index from channel index
 987       Int_t ieta3 = detId[2]/*strip*/*2 + detId[3]/*pad Z*/;
 988       if(detId[1]/*module*/ == 0) ieta3 += 0;
 989       else if(detId[1] == 1) ieta3 += 38;
 990       else if(detId[1] == 2) ieta3 += 76;
 991       else if(detId[1] == 3) ieta3 += 106;
 992       else if(detId[1] == 4) ieta3 += 144;
 993       Int_t iphi3 = detId[0]/*phi sector*/*48 + detId[4]/*pad x*/;
 994
 995       // check if the fired pad are close in space
 996       if((TMath::Abs(iphi-iphi3)>1 && TMath::Abs(iphi2-iphi3)>1) || (TMath::Abs(ieta-ieta3)>1 && TMath::Abs(ieta2-ieta3)>1))
 997 continue;
 998
 999       // check if the TOF time are close enough to be merged
1000       if(TMath::Abs(c1->GetTime() - c2->GetTime()) > 500/*in ps*/) continue;
1001
1002       // merge them
1003       Int_t label[3] = {c2->GetLabel(0),c2->GetLabel(1),c2->GetLabel(2)};
1004       fClusters[i].AddTOFhit(c2->GetClusterIndex(),chan2,c2->GetTime(),c2->GetTimeRaw(),c2->GetTOT(),label,
1005                                        c2->GetDeltaBC(),c2->GetL0L1Latency(),1,c2->GetZ(),c2->GetPhi(),c2->GetR());
1006
1007       if(c2->GetNTOFhits() > 1){ // in case also the second cluster has two hits
1008         Int_t label2[3] = {c2->GetLabel(0,1),c2->GetLabel(1,1),c2->GetLabel(2,1)};
1009         fClusters[i].AddTOFhit(c2->GetClusterIndex(1),c2->GetTOFchannel(1),c2->GetTime(1),c2->GetTimeRaw(1),
1010                                          c2->GetTOT(1),label2,c2->GetDeltaBC(2),c2->GetL0L1Latency(2),1,c2->GetZ(),
1011                                          c2->GetPhi(),c2->GetR());
1012       }
1013
1014       c1->SetStatus(0); // only the merged one should be used
1015       j = fN; // cluster "i" merged go to the next one ("i+1")
1016     }
1017   }
1018 }
1019