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