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