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