adding map include file to solve compilation issue (bug https://savannah.cern.ch...
[u/mrichter/AliRoot.git] / HLT / ITS / tracking / AliITStrackerHLT.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, 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 /* $Id: AliITStrackerHLT.cxx 32466 2009-05-20 07:51:56Z hristov $ */
16
17 //-------------------------------------------------------------------------
18 //               Implementation of the ITS tracker class
19 //    It reads AliITSRecPoint clusters and creates AliHLTITSTrack tracks
20 //                   and fills with them the ESD
21 //          Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch 
22 //          Current support and development: 
23 //                     Andrea Dainese, andrea.dainese@lnl.infn.it
24 //     dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
25 //     Params moved to AliITSRecoParam by: Andrea Dainese, INFN
26 //     Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
27 //-------------------------------------------------------------------------
28
29 //#include <TMatrixD.h>
30 #include <TTree.h>
31 #include <TDatabasePDG.h>
32 #include <TString.h>
33 #include <TRandom.h>
34 #include <TTreeStream.h>
35
36
37 #include "AliLog.h"
38 #include "AliITSCalibrationSPD.h"
39 #include "AliITSCalibrationSDD.h"
40 #include "AliITSCalibrationSSD.h"
41 #include "AliCDBEntry.h"
42 #include "AliCDBManager.h"
43 #include "AliAlignObj.h"
44 #include "AliTrackPointArray.h"
45 #include "AliESDVertex.h"
46 #include "AliESDEvent.h"
47 #include "AliESDtrack.h"
48 #include "AliV0.h"
49 #include "AliHelix.h"
50 #include "AliITSChannelStatus.h"
51 #include "AliITSRecPoint.h"
52 #include "AliITSgeomTGeo.h"
53 #include "AliITSReconstructor.h"
54 #include "AliITSClusterParam.h"
55 #include "AliITSsegmentation.h"
56 #include "AliITSCalibration.h"
57 #include "AliITSV0Finder.h"
58 #include "AliITStrackerHLT.h"
59 #include "TStopwatch.h"
60 //#include "AliHLTTPCCATrackParam.h"
61 //#include "AliHLTVertexer.h"
62 #include <vector>
63
64
65 ClassImp(AliITStrackerHLT)
66
67 Bool_t AliITStrackerHLT::TransportToX( AliExternalTrackParam *t, double x ) const
68 {
69   return t->PropagateTo( x, t->GetBz() );
70 }
71
72 Bool_t AliITStrackerHLT::TransportToPhiX( AliExternalTrackParam *t, double phi, double x ) const
73 {
74   return t->Propagate( phi, x, t->GetBz() );
75 }
76
77
78
79 AliITStrackerHLT::AliITStrackerHLT()
80   :AliTracker(),
81    fRecoParam(0),
82    fLayers(new AliHLTITSLayer[AliITSgeomTGeo::kNLayers]),
83    fUseTGeo(2),
84    fxOverX0Pipe(-1.),
85    fxTimesRhoPipe(-1.), 
86    fTracks(0),
87    fITSOutTracks(0),
88    fNTracks(0),
89    fNITSOutTracks(0),
90    fLoadTime(0),
91    fRecoTime(0),
92    fNEvents(0),
93    fClusters(0),
94    fNClusters(0)
95 {
96   //Default constructor
97   Int_t i;
98   for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
99   for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
100   for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
101 }
102 //------------------------------------------------------------------------
103 AliITStrackerHLT::AliITStrackerHLT(const Char_t *geom) 
104 : AliTracker(),
105   fRecoParam(0),
106   fLayers(new AliHLTITSLayer[AliITSgeomTGeo::kNLayers]),  
107   fUseTGeo(2),
108   fxOverX0Pipe(-1.),
109   fxTimesRhoPipe(-1.),
110   fTracks(0),
111   fITSOutTracks(0),
112   fNTracks(0),
113   fNITSOutTracks(0),
114   fLoadTime(0),
115    fRecoTime(0),
116   fNEvents(0),
117   fClusters(0),
118   fNClusters(0)
119 {
120   //--------------------------------------------------------------------
121   //This is the AliITStrackerHLT constructor
122   //--------------------------------------------------------------------
123   if (geom) {
124     AliWarning("\"geom\" is actually a dummy argument !");
125   }
126
127   if(AliGeomManager::GetGeometry()==NULL){
128     AliGeomManager::LoadGeometry();
129   }
130
131   fRecoParam = AliITSReconstructor::GetRecoParam();
132   if( !fRecoParam ){
133     AliITSReconstructor *tmp = new AliITSReconstructor();
134     tmp->Init();
135     fRecoParam = AliITSRecoParam::GetLowFluxParam();
136     tmp->AliReconstructor::SetRecoParam(fRecoParam);
137   }
138   for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
139     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
140     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
141
142     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
143     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz); 
144     Double_t poff=TMath::ATan2(y,x);
145     Double_t zoff=z;
146     Double_t r=TMath::Sqrt(x*x + y*y);
147
148     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
149     r += TMath::Sqrt(x*x + y*y);
150     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
151     r += TMath::Sqrt(x*x + y*y);
152     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
153     r += TMath::Sqrt(x*x + y*y);
154     r*=0.25;
155
156     new (fLayers+i-1) AliHLTITSLayer(r,poff,zoff,nlad,ndet);
157
158     for (Int_t j=1; j<nlad+1; j++) {
159       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
160         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
161         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
162         m.Multiply(tm);
163         Double_t txyz[3]={0.};
164         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
165         m.LocalToMaster(txyz,xyz);
166         r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
167         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
168
169         if (phi<0) phi+=TMath::TwoPi();
170         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
171
172         AliHLTITSDetector &det=fLayers[i-1].GetDetector((j-1)*ndet + k-1); 
173         new(&det) AliHLTITSDetector(r,phi); 
174         // compute the real radius (with misalignment)
175         TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
176         mmisal.Multiply(tm);
177         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
178         mmisal.LocalToMaster(txyz,xyz);
179         Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
180         det.SetRmisal(rmisal);
181         
182       } // end loop on detectors
183     } // end loop on ladders
184   } // end loop on layers
185
186   
187   Double_t xyzVtx[]={ fRecoParam->GetXVdef(),
188                       fRecoParam->GetYVdef(),
189                       fRecoParam->GetZVdef()}; 
190   Double_t ersVtx[]={ fRecoParam->GetSigmaXVdef(),
191                       fRecoParam->GetSigmaYVdef(),
192                       fRecoParam->GetSigmaZVdef()}; 
193
194   SetVertex(xyzVtx,ersVtx);
195
196   // store positions of centre of SPD modules (in z)
197   Double_t tr[3];
198   AliITSgeomTGeo::GetTranslation(1,1,1,tr);
199   fSPDdetzcentre[0] = tr[2];
200   AliITSgeomTGeo::GetTranslation(1,1,2,tr);
201   fSPDdetzcentre[1] = tr[2];
202   AliITSgeomTGeo::GetTranslation(1,1,3,tr);
203   fSPDdetzcentre[2] = tr[2];
204   AliITSgeomTGeo::GetTranslation(1,1,4,tr);
205   fSPDdetzcentre[3] = tr[2];
206
207   //fUseTGeo = fRecoParam->GetUseTGeoInTracker();
208   //if(fRecoParam->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
209   //AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
210   //fUseTGeo = 3;
211   //}
212
213   for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
214   for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
215   
216   Init();
217 }
218 //------------------------------------------------------------------------
219 AliITStrackerHLT::AliITStrackerHLT(const AliITStrackerHLT &tracker)
220 :AliTracker(tracker),
221  fRecoParam( tracker.fRecoParam),
222  fLayers(new AliHLTITSLayer[AliITSgeomTGeo::kNLayers]), 
223  fUseTGeo(tracker.fUseTGeo),
224  fxOverX0Pipe(tracker.fxOverX0Pipe),
225  fxTimesRhoPipe(tracker.fxTimesRhoPipe), 
226  fTracks(0),
227  fITSOutTracks(0),
228  fNTracks(0),
229  fNITSOutTracks(0),
230   fLoadTime(0),
231    fRecoTime(0),
232  fNEvents(0),
233  fClusters(0),
234  fNClusters(0)
235 {
236   //Copy constructor
237   Int_t i;
238   for(i=0;i<4;i++) {
239     fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
240   }
241   for(i=0;i<6;i++) {
242     fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
243     fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
244   }
245   for(i=0;i<2;i++) {
246     fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
247     fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
248   }
249   Init();
250 }
251 //------------------------------------------------------------------------
252 AliITStrackerHLT & AliITStrackerHLT::operator=(const AliITStrackerHLT &tracker){
253   //Assignment operator
254   this->~AliITStrackerHLT();
255   new(this) AliITStrackerHLT(tracker);
256   return *this;
257 }
258 //------------------------------------------------------------------------
259 AliITStrackerHLT::~AliITStrackerHLT()
260 {
261   //
262   //destructor
263   //
264   delete[] fLayers;
265   delete[] fTracks;
266   delete[] fITSOutTracks;
267   delete[] fClusters;
268 }
269
270 void AliITStrackerHLT::Init()
271 {
272   BuildMaterialLUT("Layers");  
273   BuildMaterialLUT("Pipe");  
274   BuildMaterialLUT("Shields");  
275 }
276
277
278 void AliITStrackerHLT::StartLoadClusters( Int_t NOfClusters )
279 {
280   // !
281   delete[] fClusters;
282   fClusters = new AliITSRecPoint[NOfClusters];
283   fNClusters = 0;
284 }
285
286 void AliITStrackerHLT::LoadCluster( const AliITSRecPoint &cluster) 
287 {
288   fClusters[fNClusters++] = cluster ;
289 }
290
291
292
293 //------------------------------------------------------------------------
294 Int_t AliITStrackerHLT::LoadClusters(TTree *cTree) {
295   //--------------------------------------------------------------------
296   //This function loads ITS clusters
297   //--------------------------------------------------------------------
298
299
300   TBranch *branch=cTree->GetBranch("ITSRecPoints");
301   if (!branch) { 
302     Error("LoadClusters"," can't get the branch !\n");
303     return 1;
304   }
305
306   static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
307   branch->SetAddress(&clusters);
308
309   int nClustersTotal = 0;
310   {
311     Int_t j=0;
312     for (int i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
313       int ndet=fLayers[i].GetNdetectors();
314       Int_t jmax = j + fLayers[i].GetNladders()*ndet;
315       for (; j<jmax; j++) {           
316         if (!cTree->GetEvent(j)) continue;
317         nClustersTotal+=clusters->GetEntriesFast();      
318         clusters->Delete();
319       }
320     }
321   }
322   StartLoadClusters(nClustersTotal);
323   {
324     Int_t j=0;
325     for (int i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
326       int ndet=fLayers[i].GetNdetectors();
327       Int_t jmax = j + fLayers[i].GetNladders()*ndet;
328       for (; j<jmax; j++) {           
329         if (!cTree->GetEvent(j)) continue;
330         Int_t ncl=clusters->GetEntriesFast(); 
331         while (ncl--) {
332           LoadCluster( *( (AliITSRecPoint*)clusters->UncheckedAt(ncl)));
333         }
334         clusters->Delete();
335       }
336     }
337   }
338
339   dummy.Clear();
340
341   return 0;
342 }
343
344 //------------------------------------------------------------------------
345 void AliITStrackerHLT::UnloadClusters() {
346   //--------------------------------------------------------------------
347   //This function unloads ITS clusters
348   //--------------------------------------------------------------------
349   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayers[i].ResetClusters();
350   delete[] fClusters;
351   fClusters = 0;
352   fNClusters=0;
353 }
354
355
356
357
358 void AliITStrackerHLT::Reconstruct( AliExternalTrackParam *tracksTPC, int *tracksTPCLab, int nTPCTracks )
359 {
360
361   //--------------------------------------------------------------------
362   // This functions reconstructs ITS tracks
363   //--------------------------------------------------------------------
364
365   fNEvents++;
366
367   // Init clusters
368  
369   TStopwatch timerInit;
370
371   for( int i=0; i<AliITSgeomTGeo::GetNLayers(); i++ ){ 
372     fLayers[i].ResetClusters();
373   }
374
375   for( int icl=0; icl<fNClusters; icl++ ){   
376     AliITSRecPoint &cl = fClusters[icl];
377     if (!cl.Misalign()) AliWarning("Can't misalign this cluster !"); 
378     fLayers[cl.GetLayer()].InsertCluster(&cl); 
379   }
380
381   for( int i=0; i<AliITSgeomTGeo::GetNLayers(); i++ ){ 
382     fLayers[i].ResetRoad(); //road defined by the cluster density
383     fLayers[i].SortClusters();
384   }  
385   timerInit.Stop();
386   fLoadTime+=timerInit.RealTime();
387
388
389   TStopwatch timer;
390
391   Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
392   delete[] fTracks;
393   delete[] fITSOutTracks;
394   fTracks = new AliHLTITSTrack[nTPCTracks];
395   fITSOutTracks = new AliHLTITSTrack[nTPCTracks];
396   fNTracks = 0;
397   fNITSOutTracks = 0;
398   for( int itr=0; itr<nTPCTracks; itr++ ){    
399
400     AliHLTITSTrack tMI( tracksTPC[itr] );
401     AliHLTITSTrack *t = &tMI;
402     t->SetTPCtrackId( itr );
403     t->SetMass(pimass); 
404     t->SetExpQ(0);
405     t->SetLabel(tracksTPCLab[itr]);
406
407     //if (!CorrectForTPCtoITSDeadZoneMaterial(t))  continue;
408       
409     Int_t tpcLabel=t->GetLabel(); //save the TPC track label       
410     
411     FollowProlongationTree(t); 
412     int nclu=0;
413     for(Int_t i=0; i<6; i++) {
414       if( t->GetClusterIndex(i)>=0 ) nclu++; 
415     }
416     //cout<<"N assigned ITS clusters = "<<nclu<<std::endl;
417     t->SetLabel(-1);
418     if( nclu>0 ){
419       t->SetLabel(tpcLabel);
420       t->SetFakeRatio(1.);
421       CookLabel(t,.99); //For comparison only
422       //cout<<"SG: label = "<<t->GetLabel()<<" / "<<tpcLabel<<endl;
423     }
424
425     CorrectForPipeMaterial(t);
426    
427     TransportToX(t, 0 );
428     fTracks[fNTracks++] = *t;  
429     //cout<<"SG: ITS: Bz = "<<t->GetBz()<<endl;
430
431     if(  nclu>0 ){ // construct ITSOut track
432       AliHLTITSTrack tOut(*t);
433       if( FitOutward( &tOut ) ){
434         fITSOutTracks[fNITSOutTracks++] = *t;  
435       }
436     }
437   }
438
439   timer.Stop();
440   fRecoTime+=timer.RealTime();
441 }
442
443
444
445 //------------------------------------------------------------------------
446 Int_t AliITStrackerHLT::Clusters2Tracks(AliESDEvent *event) {
447   //--------------------------------------------------------------------
448   // This functions reconstructs ITS tracks
449   // The clusters must be already loaded !
450   //--------------------------------------------------------------------
451   
452   
453   std::vector<AliExternalTrackParam> tracksTPC;
454   std::vector<int> tracksTPCLab;
455   tracksTPC.reserve(event->GetNumberOfTracks());
456
457   for( int itr=0; itr<event->GetNumberOfTracks(); itr++ ){
458
459     AliESDtrack *esdTrack = event->GetTrack(itr);
460     //esdTrack->myITS = esdTrack->myTPC;
461     if ((esdTrack->GetStatus()&AliESDtrack::kTPCin)==0) continue;
462     //if (esdTrack->GetStatus()&AliESDtrack::kTPCout) continue;
463     if (esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
464     if (esdTrack->GetKinkIndex(0)>0) continue;   //kink daughter
465     
466     AliHLTITSTrack t(*esdTrack);
467     t.SetTPCtrackId( itr );
468     tracksTPC.push_back( t );
469     tracksTPCLab.push_back(esdTrack->GetLabel());
470   }
471   //for( int iter=0; iter<100; iter++){
472   Reconstruct( &(tracksTPC[0]), &(tracksTPCLab[0]), tracksTPC.size() );
473   //}
474
475   for( int itr=0; itr<fNTracks; itr++ ){
476     AliHLTITSTrack &t = fTracks[itr];    
477     UpdateESDtrack(event->GetTrack(t.TPCtrackId()), &t, AliESDtrack::kITSin);          
478     //event->GetTrack(t.TPCtrackId())->myITS = t;
479   }
480  
481
482   //int hz = ( int ) ( (fRecoTime+fLoadTime) > 1.e-4 ? fNEvents / (fRecoTime+fLoadTime) : 0 );
483   //int hz1 = ( int ) ( fRecoTime > 1.e-4 ? fNEvents / fRecoTime : 0 );
484   //int hz2 = ( int ) ( fLoadTime > 1.e-4 ? fNEvents / fLoadTime : 0 );
485
486   //std::cout<<"\n\nSG: ITS tracker time = "<<hz2<<" Hz load / "<<hz1<<" Hz reco ="
487   //<<hz<<
488   //" Hz ("
489   //<<fLoadTime/fNEvents*1000<<"+"<<fRecoTime/fNEvents*1000.
490   //<<" = "<<(fLoadTime + fRecoTime)/fNEvents*1000. 
491   //<<" ms/ev), "<<fNEvents<<" events processed\n\n "<<std::endl;
492   return 0;
493 }
494
495
496 AliCluster *AliITStrackerHLT::GetCluster(Int_t index) const 
497 {
498   //       Return pointer to a given cluster
499   Int_t l=(index & 0xf0000000) >> 28;
500   Int_t c=(index & 0x0fffffff) >> 00;
501   return fLayers[l].GetCluster(c);
502 }
503
504
505
506
507 //------------------------------------------------------------------------
508 void AliITStrackerHLT::FollowProlongationTree(AliHLTITSTrack * track ) 
509 {
510   // FollowProlongationTree
511   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
512    
513     AliHLTITSLayer &layer=fLayers[ilayer];
514   
515     // material between SSD and SDD, SDD and SPD
516     //if (ilayer==3 && !CorrectForShieldMaterial(track,1)) continue;
517     //if (ilayer==1 && !CorrectForShieldMaterial(track,0)) continue;
518     
519     int idet;
520
521     {            
522       Double_t xloc, phi,z;
523       if( !track->GetLocalXPhiZat( layer.GetR(), xloc, phi, z ) ) return;
524       idet = layer.FindDetectorIndex(phi,z);
525     }
526
527     // track outside layer acceptance in z
528    
529     if( idet<0 ) continue;
530     
531     // propagate to the intersection with the detector plane     
532     {
533       const AliHLTITSDetector &det=layer.GetDetector( idet );
534       if (!TransportToPhiX( track, det.GetPhi(), det.GetR() ) ) return;
535       CorrectForLayerMaterial(track,ilayer);
536     }
537
538     // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
539     
540     // road in global (rphi,z) [i.e. in tracking ref. system]
541     
542     Double_t zmin,zmax,ymin,ymax;
543
544     if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
545   
546     layer.SelectClusters(zmin,zmax,ymin,ymax);     
547     
548     // Define criteria for track-cluster association
549     
550     Double_t msz = track->GetSigmaZ2() + 
551       fRecoParam->GetNSigmaZLayerForRoadZ()*
552       fRecoParam->GetNSigmaZLayerForRoadZ()*
553       fRecoParam->GetSigmaZ2(ilayer);
554     
555     Double_t msy = track->GetSigmaY2() + 
556       fRecoParam->GetNSigmaYLayerForRoadY()*
557       fRecoParam->GetNSigmaYLayerForRoadY()*
558       fRecoParam->GetSigmaY2(ilayer);
559     
560     msz *= fRecoParam->GetNSigma2RoadZNonC();
561     msy *= fRecoParam->GetNSigma2RoadYNonC(); 
562     
563     msz = 1./msz; // 1/RoadZ^2
564     msy = 1./msy; // 1/RoadY^2    
565     
566     const AliITSRecPoint *cl=0; 
567     Int_t clidx=-1;     
568     
569     // loop over clusters in the road     
570     
571     const AliITSRecPoint *bestCluster=0; 
572     double bestChi2 = 1.e10;
573     AliHLTITSTrack bestTrack( *track );
574     int bestIdx = -1;
575      
576     while( (cl=layer.GetNextCluster(clidx)) ){
577       Int_t idetc=cl->GetDetectorIndex();
578       if ( idet !=idetc ) { // new cluster's detector
579         const AliHLTITSDetector &detc=layer.GetDetector(idetc);
580         if (!TransportToPhiX( track, detc.GetPhi(),detc.GetR()) ) continue;
581         idet = idetc;
582       }  
583       //double y,z;
584       //if (! track->GetLocalYZat( layer.GetDetector(idetc).GetR() + cl->GetX(),y,z ) ) continue;
585       double dz = track->GetZ() - cl->GetZ();
586       double dy = track->GetY() - cl->GetY();
587       double chi2 = dz*dz*msz + dy*dy*msy ;       
588       if ( chi2 < bestChi2 ){
589         bestChi2 = chi2;
590         bestCluster = cl;
591         bestTrack = *track;
592         bestIdx = clidx;
593         continue;
594       }
595     }
596     
597     if( !bestCluster || bestChi2 >2*10. ) continue;
598     
599     if (!TransportToX( &bestTrack, layer.GetDetector(bestCluster->GetDetectorIndex()).GetR() + bestCluster->GetX() ) ) continue;
600     
601     Double_t par[2]={ bestCluster->GetY(), bestCluster->GetZ()};
602     Double_t cov[3]={ bestCluster->GetSigmaY2(), 0., bestCluster->GetSigmaZ2()};
603     if( !bestTrack.AliExternalTrackParam::Update(par,cov) ) continue;
604     
605     *track = bestTrack;
606     track->SetClusterIndex(track->GetNumberOfClusters(), (ilayer<<28)+bestIdx);
607     track->SetNumberOfClusters(track->GetNumberOfClusters()+1);  
608   }
609 }
610
611
612
613 Int_t AliITStrackerHLT::FitOutward(AliHLTITSTrack * track ) 
614 {
615   // FitOutward
616   track->ResetCovariance(100);
617
618   for (Int_t iTrCl=track->GetNumberOfClusters()-1; iTrCl>=0; iTrCl--) {
619     
620     Int_t index = track->GetClusterIndex(iTrCl);
621     Int_t ilayer=(index & 0xf0000000) >> 28;
622     Int_t ic=(index & 0x0fffffff) >> 00;
623     const AliHLTITSLayer &layer=fLayers[ilayer];
624     AliITSRecPoint *cl = layer.GetCluster(ic); 
625     int idet = cl->GetDetectorIndex();
626     const AliHLTITSDetector &det=layer.GetDetector( idet );
627  
628     // material between SSD and SDD, SDD and SPD
629     //if (ilayer==4 && !CorrectForShieldMaterial(track,1)) continue;
630     //if (ilayer==2 && !CorrectForShieldMaterial(track,0)) continue;
631     
632
633     // propagate to the intersection with the detector plane     
634     {
635       if (!TransportToPhiX( track, det.GetPhi(), det.GetR()+ cl->GetX() ) ) return 0;
636       CorrectForLayerMaterial(track,ilayer);
637     }
638
639     Double_t par[2]={ cl->GetY(), cl->GetZ()};
640     Double_t cov[3]={ cl->GetSigmaY2(), 0., cl->GetSigmaZ2()};
641     if( !track->AliExternalTrackParam::Update(par,cov) ) return 0;    
642   }
643   return 1;
644 }
645
646
647 //------------------------------------------------------------------------
648 AliHLTITSLayer & AliITStrackerHLT::GetLayer(Int_t layer) const
649 {
650   //--------------------------------------------------------------------
651   //
652   //
653   return fLayers[layer];
654 }
655
656
657
658 //------------------------------------------------------------------------
659 void AliITStrackerHLT::CookLabel(AliHLTITSTrack *track,Float_t wrong) const 
660 {
661   // get MC label for the track
662
663   Int_t mcLabel = -1;
664   
665   vector<int> labels;
666   Int_t nClusters = track->GetNumberOfClusters();
667   Int_t nClustersEff = 0;
668   for (Int_t i=0; i<nClusters; i++){
669     Int_t cindex = track->GetClusterIndex(i);
670     //Int_t l=(cindex & 0xf0000000) >> 28;
671     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);    
672     if ( cl->GetLabel(0) >= 0 ){ labels.push_back(cl->GetLabel(0)) ; nClustersEff++; }
673     if ( cl->GetLabel(1) >= 0 ) labels.push_back(cl->GetLabel(1)) ;
674     if ( cl->GetLabel(2) >= 0 ) labels.push_back(cl->GetLabel(2)) ;
675   }
676   std::sort( labels.begin(), labels.end() );
677
678   labels.push_back( -1 ); // put -1 to the end    
679   int labelMax = -1, labelCur = -1, nLabelsMax = 0, nLabelsCurr = 0;
680
681   for ( unsigned int iLab = 0; iLab < labels.size(); iLab++ ) {
682     if ( labels[iLab] != labelCur ) {         
683       if ( labelCur >= 0 && nLabelsMax< nLabelsCurr ) {
684         nLabelsMax = nLabelsCurr;
685         labelMax = labelCur;
686       }
687       labelCur = labels[iLab];
688       nLabelsCurr = 0;
689     }
690     nLabelsCurr++;
691   }
692
693   if( labelMax>=0 && nLabelsMax < wrong * nClustersEff ) labelMax = -labelMax;
694   
695   mcLabel = labelMax;
696                 
697   track->SetLabel( mcLabel );
698 }
699
700
701
702
703
704
705
706
707
708 //------------------------------------------------------------------------
709 void AliITStrackerHLT::UpdateESDtrack(AliESDtrack *tESD, AliHLTITSTrack* track, ULong_t flags) const
710 {
711   //
712   // Update ESD track
713   //
714   tESD->UpdateTrackParams(track,flags);
715   AliHLTITSTrack * oldtrack = (AliHLTITSTrack*)(tESD->GetITStrack());
716   if (oldtrack) delete oldtrack; 
717   tESD->SetITStrack(new AliHLTITSTrack(*track));
718 }
719
720
721
722
723 //------------------------------------------------------------------------
724 void AliITStrackerHLT::BuildMaterialLUT(TString material) {
725   //--------------------------------------------------------------------
726   // Fill a look-up table with mean material
727   //--------------------------------------------------------------------
728
729   Int_t n=1000;
730   Double_t mparam[7];
731   Double_t point1[3],point2[3];
732   Double_t phi,cosphi,sinphi,z;
733   // 0-5 layers, 6 pipe, 7-8 shields 
734   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
735   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
736
737   Int_t ifirst=0,ilast=0;  
738   if(material.Contains("Pipe")) {
739     ifirst=6; ilast=6;
740   } else if(material.Contains("Shields")) {
741     ifirst=7; ilast=8;
742   } else if(material.Contains("Layers")) {
743     ifirst=0; ilast=5;
744   } else {
745     Error("BuildMaterialLUT","Wrong layer name\n");
746   }
747
748   for(Int_t imat=ifirst; imat<=ilast; imat++) {
749     Double_t param[5]={0.,0.,0.,0.,0.};
750     for (Int_t i=0; i<n; i++) {
751       phi = 2.*TMath::Pi()*gRandom->Rndm();
752       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi); 
753       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
754       point1[0] = rmin[imat]*cosphi;
755       point1[1] = rmin[imat]*sinphi;
756       point1[2] = z;
757       point2[0] = rmax[imat]*cosphi;
758       point2[1] = rmax[imat]*sinphi;
759       point2[2] = z;
760       AliTracker::MeanMaterialBudget(point1,point2,mparam);
761       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
762     }
763     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
764     if(imat<=5) {
765       fxOverX0Layer[imat] = param[1];
766       fxTimesRhoLayer[imat] = param[0]*param[4];
767     } else if(imat==6) {
768       fxOverX0Pipe = param[1];
769       fxTimesRhoPipe = param[0]*param[4];
770     } else if(imat==7) {
771       fxOverX0Shield[0] = param[1];
772       fxTimesRhoShield[0] = param[0]*param[4];
773     } else if(imat==8) {
774       fxOverX0Shield[1] = param[1];
775       fxTimesRhoShield[1] = param[0]*param[4];
776     }
777   }
778   /*
779   printf("%s\n",material.Data());
780   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
781   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
782   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
783   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
784   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
785   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
786   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
787   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
788   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
789   */
790   return;
791 }
792
793
794
795
796 //------------------------------------------------------------------------
797 Int_t AliITStrackerHLT::CorrectForTPCtoITSDeadZoneMaterial(AliHLTITSTrack *t) {
798   //--------------------------------------------------------------------
799   // Correction for the material between the TPC and the ITS
800   //--------------------------------------------------------------------
801   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction 
802       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
803       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
804       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
805   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
806       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
807       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
808       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
809   } else {
810     printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
811     return 0;
812   }
813   
814   return 1;
815 }
816
817
818 //------------------------------------------------------------------------
819 Int_t AliITStrackerHLT::CorrectForPipeMaterial(AliHLTITSTrack *t,
820                                                bool InwardDirection) {
821   //-------------------------------------------------------------------
822   // Propagate beyond beam pipe and correct for material
823   // (material budget in different ways according to fUseTGeo value)
824   // Add time if going outward (PropagateTo or PropagateToTGeo)
825   //-------------------------------------------------------------------
826
827   // Define budget mode:
828   // 0: material from AliITSRecoParam (hard coded)
829   // 1: material from TGeo in one step (on the fly)
830   // 2: material from lut
831   // 3: material from TGeo in one step (same for all hypotheses)
832   Int_t mode;
833   switch(fUseTGeo) {
834   case 0:
835     mode=0; 
836     break;    
837   case 1:
838     mode=1;
839     break;    
840   case 2:
841     mode=2;
842     break;
843   case 3:
844     mode=3; 
845     break;
846   case 4:
847     mode=3;
848     break;
849   default:
850     mode=0;
851     break;
852   }
853   
854   Float_t  dir = (InwardDirection ? 1. : -1.);
855   Double_t rToGo= ( InwardDirection ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
856   Double_t xToGo, phi,z;
857
858   if (!t->GetLocalXPhiZat(rToGo,xToGo,phi,z)) return 0;
859
860   Double_t xOverX0,x0,lengthTimesMeanDensity;
861
862   switch(mode) {
863   case 0:
864     xOverX0 = AliITSRecoParam::GetdPipe();
865     x0 = AliITSRecoParam::GetX0Be();
866     lengthTimesMeanDensity = xOverX0*x0;
867     lengthTimesMeanDensity *= dir;
868     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
869     break;
870   case 1:
871     if (!t->PropagateToTGeo(xToGo,1)) return 0;
872     break;
873   case 2:
874     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");  
875     xOverX0 = fxOverX0Pipe;
876     lengthTimesMeanDensity = fxTimesRhoPipe;
877     lengthTimesMeanDensity *= dir;
878     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
879     break;
880   case 3:
881     double xOverX0PipeTrks, xTimesRhoPipeTrks;
882     if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
883     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
884                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
885     xOverX0PipeTrks = TMath::Abs(xOverX0)/angle;
886     xTimesRhoPipeTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
887     xOverX0 = xOverX0PipeTrks;
888     lengthTimesMeanDensity = xTimesRhoPipeTrks;
889     lengthTimesMeanDensity *= dir;
890     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
891     break;
892   }
893   
894   return 1;
895 }
896 //------------------------------------------------------------------------
897 Int_t AliITStrackerHLT::CorrectForShieldMaterial(AliHLTITSTrack *t,
898                                                  Int_t    shieldindex,
899                                                 bool InwardDirection) {
900   //-------------------------------------------------------------------
901   // Propagate beyond SPD or SDD shield and correct for material
902   // (material budget in different ways according to fUseTGeo value)
903   // Add time if going outward (PropagateTo or PropagateToTGeo)
904   //-------------------------------------------------------------------
905
906   // Define budget mode:
907   // 0: material from AliITSRecoParam (hard coded)
908   // 1: material from TGeo in steps of X cm (on the fly)
909   //    X = AliITSRecoParam::GetStepSizeTGeo()
910   // 2: material from lut
911   // 3: material from TGeo in one step (same for all hypotheses)
912   Int_t mode;
913   switch(fUseTGeo) {
914   case 0:
915     mode=0; 
916     break;    
917   case 1:
918     mode=1;
919     break;    
920   case 2:
921     mode=2;
922     break;
923   case 3:
924     mode=3;
925     break;
926   case 4:
927     mode=3;
928     break;
929   default:
930     mode=0;
931     break;
932   }
933
934
935   Float_t  dir = (InwardDirection ? 1. : -1.);
936   Double_t rToGo;
937
938   if (shieldindex==1 ) { // SDDouter
939     rToGo=(InwardDirection ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
940   } else if (shieldindex==0 ) {        // SPDouter
941     rToGo=(InwardDirection ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0)); 
942   } else {
943     Error("CorrectForShieldMaterial"," Wrong shield name\n");
944     return 0;
945   }
946   Double_t xToGo, phi,z;
947
948   if (!t->GetLocalXPhiZat(rToGo,xToGo,phi,z)) return 0;
949
950   Double_t xOverX0,x0,lengthTimesMeanDensity;
951   Int_t nsteps=1;
952
953   switch(mode) {
954   case 0:
955     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
956     x0 = AliITSRecoParam::GetX0shield(shieldindex);
957     lengthTimesMeanDensity = xOverX0*x0;
958     lengthTimesMeanDensity *= dir;
959     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
960     break;
961   case 1:
962     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/fRecoParam->GetStepSizeTGeo())+1;
963     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
964     break;
965   case 2:
966     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");  
967     xOverX0 = fxOverX0Shield[shieldindex];
968     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
969     lengthTimesMeanDensity *= dir;
970     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
971     break;
972   case 3:    
973     if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
974     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
975                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
976     double xOverX0ShieldTrks = TMath::Abs(xOverX0)/angle;
977     double xTimesRhoShieldTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
978     xOverX0 = xOverX0ShieldTrks;
979     lengthTimesMeanDensity = xTimesRhoShieldTrks;
980     lengthTimesMeanDensity *= dir;
981     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
982     break;
983   }
984
985   return 1;
986 }
987 //------------------------------------------------------------------------
988 Int_t AliITStrackerHLT::CorrectForLayerMaterial(AliHLTITSTrack *t,
989                                                Int_t layerindex,
990                                                bool InwardDirection ){
991   //-------------------------------------------------------------------
992   // Propagate beyond layer and correct for material
993   // (material budget in different ways according to fUseTGeo value)
994   // Add time if going outward (PropagateTo or PropagateToTGeo)
995   //-------------------------------------------------------------------
996
997   /*
998     Double_t r=fLayers[layerindex].GetR();
999     Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
1000     Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
1001     rToGo+= InwardDirection ?-deltar :deltar;
1002     Double_t xToGo;
1003     if (!t->GetLocalXat(rToGo,xToGo)) return 0;  
1004   */
1005   
1006   if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
1007
1008   Double_t lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
1009   if( !InwardDirection ) lengthTimesMeanDensity = -lengthTimesMeanDensity;
1010
1011   return t->CorrectForMeanMaterial(fxOverX0Layer[layerindex],lengthTimesMeanDensity,kTRUE);
1012 }
1013
1014
1015 //------------------------------------------------------------------------
1016 Bool_t AliITStrackerHLT::LocalModuleCoord(Int_t ilayer,Int_t idet,
1017                                        const AliHLTITSTrack *track,
1018                                        Float_t &xloc,Float_t &zloc) const {
1019   //-----------------------------------------------------------------
1020   // Gives position of track in local module ref. frame
1021   //-----------------------------------------------------------------
1022
1023   xloc=0.; 
1024   zloc=0.;
1025
1026   if(idet<0) return kFALSE;
1027
1028   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6 
1029
1030   Int_t lad = Int_t(idet/ndet) + 1;
1031
1032   Int_t det = idet - (lad-1)*ndet + 1;
1033
1034   Double_t xyzGlob[3],xyzLoc[3];
1035
1036   AliHLTITSDetector &detector = fLayers[ilayer].GetDetector(idet);
1037   // take into account the misalignment: xyz at real detector plane
1038   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
1039
1040   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
1041
1042   xloc = (Float_t)xyzLoc[0];
1043   zloc = (Float_t)xyzLoc[2];
1044
1045   return kTRUE;
1046 }
1047
1048 //------------------------------------------------------------------------
1049 Bool_t AliITStrackerHLT::ComputeRoad(AliHLTITSTrack* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1050   //--------------------------------------------------------------------
1051   // This function computes the rectangular road for this track
1052   //--------------------------------------------------------------------
1053
1054
1055   AliHLTITSDetector &det = fLayers[ilayer].GetDetector(idet);
1056   // take into account the misalignment: propagate track to misaligned detector plane
1057
1058   double y,z,snp,cov[3];
1059   if( !track->GetYZAtPhiX( det.GetPhi(),det.GetRmisal(), y, z, snp, cov))return 0;
1060
1061
1062   Double_t dz=fRecoParam->GetNSigmaRoadZ()*
1063                     TMath::Sqrt(cov[2] + 
1064                                 fRecoParam->GetNSigmaZLayerForRoadZ()*
1065                                 fRecoParam->GetNSigmaZLayerForRoadZ()*
1066                                 fRecoParam->GetSigmaZ2(ilayer));
1067   Double_t dy=fRecoParam->GetNSigmaRoadY()*
1068                     TMath::Sqrt(cov[0] + 
1069                                 fRecoParam->GetNSigmaYLayerForRoadY()*
1070                                 fRecoParam->GetNSigmaYLayerForRoadY()*
1071                                 fRecoParam->GetSigmaY2(ilayer));
1072       
1073   // track at boundary between detectors, enlarge road
1074   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1075   if ( (y-dy < det.GetYmin()+boundaryWidth) || 
1076        (y+dy > det.GetYmax()-boundaryWidth) || 
1077        (z-dz < det.GetZmin()+boundaryWidth) ||
1078        (z+dz > det.GetZmax()-boundaryWidth) ) {
1079     Float_t tgl = TMath::Abs(track->GetTgl());
1080     if (tgl > 1.) tgl=1.;
1081     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1082     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1083     if (snp > fRecoParam->GetMaxSnp()) return kFALSE;
1084     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1085   } // boundary
1086   
1087   // add to the road a term (up to 2-3 mm) to deal with misalignments
1088   dy = TMath::Sqrt(dy*dy + fRecoParam->GetRoadMisal()*fRecoParam->GetRoadMisal());
1089   dz = TMath::Sqrt(dz*dz + fRecoParam->GetRoadMisal()*fRecoParam->GetRoadMisal());
1090
1091   Double_t r = fLayers[ilayer].GetR();
1092   zmin = z - dz; 
1093   zmax = z + dz;
1094   ymin = y + r*det.GetPhi() - dy;
1095   ymax = y + r*det.GetPhi() + dy;
1096
1097   return kTRUE;
1098 }
1099
1100
1101
1102 Int_t AliITStrackerHLT::PropagateBack(AliESDEvent * /*event*/) 
1103 {
1104   // dummy
1105   return 0;
1106 }
1107
1108 Int_t AliITStrackerHLT::RefitInward(AliESDEvent * /*event*/ ) 
1109 {
1110   // dummy
1111   return 0;
1112 }
1113
1114
1115 Bool_t AliITStrackerHLT::GetTrackPoint(Int_t /*index*/, AliTrackPoint& /*p*/) const 
1116 {
1117   // dummy
1118   return 0;
1119 }
1120
1121 Bool_t AliITStrackerHLT::GetTrackPointTrackingError(Int_t /*index*/, 
1122                                                     AliTrackPoint& /*p*/, const AliESDtrack */*t*/) 
1123 {
1124   // dummy
1125   return 0;
1126 }