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