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