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