compilation error fixed
[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 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
406     //if (!CorrectForTPCtoITSDeadZoneMaterial(t))  continue;
407       
408     //Int_t tpcLabel=t->GetLabel(); //save the TPC track label       
409       
410     FollowProlongationTree(t); 
411     int nclu=0;
412     for(Int_t i=0; i<6; i++) {
413       if( t->GetClusterIndex(i)>=0 ) nclu++; 
414     }
415     //cout<<"N assigned ITS clusters = "<<nclu<<std::endl;
416     t->SetLabel(-1);
417     if( nclu>0 ){
418       t->SetLabel(-1);//tpcLabel);
419       t->SetFakeRatio(1.);
420       CookLabel(t,.99); //For comparison only
421       //cout<<"label = "<<t->GetLabel()<<" / "<<tpcLabel<<endl;
422     }
423
424     //CorrectForPipeMaterial(t);
425    
426     TransportToX(t, 0 );
427     fTracks[fNTracks++] = *t;  
428     if(  nclu>0 ){ // construct ITSOut track
429       AliHLTITSTrack tOut(*t);
430       if( FitOutward( &tOut ) ){
431         fITSOutTracks[fNITSOutTracks++] = *t;  
432       }
433     }
434   }
435
436   timer.Stop();
437   fRecoTime+=timer.RealTime();
438 }
439
440
441
442 //------------------------------------------------------------------------
443 Int_t AliITStrackerHLT::Clusters2Tracks(AliESDEvent *event) {
444   //--------------------------------------------------------------------
445   // This functions reconstructs ITS tracks
446   // The clusters must be already loaded !
447   //--------------------------------------------------------------------
448   
449   
450   std::vector<AliExternalTrackParam> tracksTPC;
451   tracksTPC.reserve(event->GetNumberOfTracks());
452
453   for( int itr=0; itr<event->GetNumberOfTracks(); itr++ ){
454
455     AliESDtrack *esdTrack = event->GetTrack(itr);
456     //esdTrack->myITS = esdTrack->myTPC;
457     if ((esdTrack->GetStatus()&AliESDtrack::kTPCin)==0) continue;
458     //if (esdTrack->GetStatus()&AliESDtrack::kTPCout) continue;
459     if (esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
460     if (esdTrack->GetKinkIndex(0)>0) continue;   //kink daughter
461     
462     AliHLTITSTrack t(*esdTrack);
463     t.SetTPCtrackId( itr );
464     tracksTPC.push_back( t );
465   }
466   //for( int iter=0; iter<100; iter++){
467   Reconstruct( &(tracksTPC[0]), tracksTPC.size() );
468   //}
469
470   for( int itr=0; itr<fNTracks; itr++ ){
471     AliHLTITSTrack &t = fTracks[itr];    
472     UpdateESDtrack(event->GetTrack(t.TPCtrackId()), &t, AliESDtrack::kITSin);          
473     //event->GetTrack(t.TPCtrackId())->myITS = t;
474   }
475  
476
477   int hz = ( int ) ( (fRecoTime+fLoadTime) > 1.e-4 ? fNEvents / (fRecoTime+fLoadTime) : 0 );
478   int hz1 = ( int ) ( fRecoTime > 1.e-4 ? fNEvents / fRecoTime : 0 );
479   int hz2 = ( int ) ( fLoadTime > 1.e-4 ? fNEvents / fLoadTime : 0 );
480
481   std::cout<<"\n\n ITS tracker time = "<<hz2<<" Hz load / "<<hz1<<" Hz reco ="
482            <<hz<<
483     " Hz ("
484            <<fLoadTime/fNEvents*1000<<"+"<<fRecoTime/fNEvents*1000.
485            <<" = "<<(fLoadTime + fRecoTime)/fNEvents*1000.
486            <<" ms/ev), "<<fNEvents<<" events processed\n\n "<<std::endl;
487   return 0;
488 }
489
490
491 AliCluster *AliITStrackerHLT::GetCluster(Int_t index) const 
492 {
493   //       Return pointer to a given cluster
494   Int_t l=(index & 0xf0000000) >> 28;
495   Int_t c=(index & 0x0fffffff) >> 00;
496   return fLayers[l].GetCluster(c);
497 }
498
499
500
501
502 //------------------------------------------------------------------------
503 void AliITStrackerHLT::FollowProlongationTree(AliHLTITSTrack * track ) 
504 {
505   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
506    
507     AliHLTITSLayer &layer=fLayers[ilayer];
508   
509     // material between SSD and SDD, SDD and SPD
510     //if (ilayer==3 && !CorrectForShieldMaterial(track,1)) continue;
511     //if (ilayer==1 && !CorrectForShieldMaterial(track,0)) continue;
512     
513     int idet;
514
515     {            
516       Double_t xloc, phi,z;
517       if( !track->GetLocalXPhiZat( layer.GetR(), xloc, phi, z ) ) return;
518       idet = layer.FindDetectorIndex(phi,z);
519     }
520
521     // track outside layer acceptance in z
522    
523     if( idet<0 ) continue;
524     
525     // propagate to the intersection with the detector plane     
526     {
527       const AliHLTITSDetector &det=layer.GetDetector( idet );
528       if (!TransportToPhiX( track, det.GetPhi(), det.GetR() ) ) return;
529       CorrectForLayerMaterial(track,ilayer);
530     }
531
532     // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
533     
534     // road in global (rphi,z) [i.e. in tracking ref. system]
535     
536     Double_t zmin,zmax,ymin,ymax;
537
538     if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
539   
540     layer.SelectClusters(zmin,zmax,ymin,ymax);     
541     
542     // Define criteria for track-cluster association
543     
544     Double_t msz = track->GetSigmaZ2() + 
545       fRecoParam->GetNSigmaZLayerForRoadZ()*
546       fRecoParam->GetNSigmaZLayerForRoadZ()*
547       fRecoParam->GetSigmaZ2(ilayer);
548     
549     Double_t msy = track->GetSigmaY2() + 
550       fRecoParam->GetNSigmaYLayerForRoadY()*
551       fRecoParam->GetNSigmaYLayerForRoadY()*
552       fRecoParam->GetSigmaY2(ilayer);
553     
554     msz *= fRecoParam->GetNSigma2RoadZNonC();
555     msy *= fRecoParam->GetNSigma2RoadYNonC(); 
556     
557     msz = 1./msz; // 1/RoadZ^2
558     msy = 1./msy; // 1/RoadY^2    
559     
560     const AliITSRecPoint *cl=0; 
561     Int_t clidx=-1;     
562     
563     // loop over clusters in the road     
564     
565     const AliITSRecPoint *bestCluster=0; 
566     double bestChi2 = 1.e10;
567     AliHLTITSTrack bestTrack( *track );
568     int bestIdx = -1;
569      
570     while( (cl=layer.GetNextCluster(clidx)) ){
571       Int_t idetc=cl->GetDetectorIndex();
572       if ( idet !=idetc ) { // new cluster's detector
573         const AliHLTITSDetector &detc=layer.GetDetector(idetc);
574         if (!TransportToPhiX( track, detc.GetPhi(),detc.GetR()) ) continue;
575         idet = idetc;
576       }  
577       //double y,z;
578       //if (! track->GetLocalYZat( layer.GetDetector(idetc).GetR() + cl->GetX(),y,z ) ) continue;
579       double dz = track->GetZ() - cl->GetZ();
580       double dy = track->GetY() - cl->GetY();
581       double chi2 = dz*dz*msz + dy*dy*msy ;       
582       if ( chi2 < bestChi2 ){
583         bestChi2 = chi2;
584         bestCluster = cl;
585         bestTrack = *track;
586         bestIdx = clidx;
587         continue;
588       }
589     }
590     
591     if( !bestCluster || bestChi2 >2*10. ) continue;
592     
593     if (!TransportToX( &bestTrack, layer.GetDetector(bestCluster->GetDetectorIndex()).GetR() + bestCluster->GetX() ) ) continue;
594     
595     Double_t par[2]={ bestCluster->GetY(), bestCluster->GetZ()};
596     Double_t cov[3]={ bestCluster->GetSigmaY2(), 0., bestCluster->GetSigmaZ2()};
597     if( !bestTrack.AliExternalTrackParam::Update(par,cov) ) continue;
598     
599     *track = bestTrack;
600     track->SetClusterIndex(track->GetNumberOfClusters(), (ilayer<<28)+bestIdx);
601     track->SetNumberOfClusters(track->GetNumberOfClusters()+1);  
602   }
603 }
604
605
606
607 Int_t AliITStrackerHLT::FitOutward(AliHLTITSTrack * track ) 
608 {
609
610   track->ResetCovariance(100);
611
612   for (Int_t iTrCl=track->GetNumberOfClusters()-1; iTrCl>=0; iTrCl--) {
613     
614     Int_t index = track->GetClusterIndex(iTrCl);
615     Int_t ilayer=(index & 0xf0000000) >> 28;
616     Int_t ic=(index & 0x0fffffff) >> 00;
617     const AliHLTITSLayer &layer=fLayers[ilayer];
618     AliITSRecPoint *cl = layer.GetCluster(ic); 
619     int idet = cl->GetDetectorIndex();
620     const AliHLTITSDetector &det=layer.GetDetector( idet );
621  
622     // material between SSD and SDD, SDD and SPD
623     //if (ilayer==4 && !CorrectForShieldMaterial(track,1)) continue;
624     //if (ilayer==2 && !CorrectForShieldMaterial(track,0)) continue;
625     
626
627     // propagate to the intersection with the detector plane     
628     {
629       if (!TransportToPhiX( track, det.GetPhi(), det.GetR()+ cl->GetX() ) ) return 0;
630       CorrectForLayerMaterial(track,ilayer);
631     }
632
633     Double_t par[2]={ cl->GetY(), cl->GetZ()};
634     Double_t cov[3]={ cl->GetSigmaY2(), 0., cl->GetSigmaZ2()};
635     if( !track->AliExternalTrackParam::Update(par,cov) ) return 0;    
636   }
637   return 1;
638 }
639
640
641 //------------------------------------------------------------------------
642 AliHLTITSLayer & AliITStrackerHLT::GetLayer(Int_t layer) const
643 {
644   //--------------------------------------------------------------------
645   //
646   //
647   return fLayers[layer];
648 }
649
650
651
652 //------------------------------------------------------------------------
653 void AliITStrackerHLT::CookLabel(AliHLTITSTrack *track,Float_t wrong) const 
654 {
655   // get MC label for the track
656
657   Int_t mcLabel = -1;
658   
659   vector<int> labels;
660   Int_t nClusters = track->GetNumberOfClusters();
661
662   for (Int_t i=0; i<nClusters; i++){
663     Int_t cindex = track->GetClusterIndex(i);
664     //Int_t l=(cindex & 0xf0000000) >> 28;
665     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);    
666     if ( cl->GetLabel(0) >= 0 ) labels.push_back(cl->GetLabel(0)) ;
667     if ( cl->GetLabel(1) >= 0 ) labels.push_back(cl->GetLabel(1)) ;
668     if ( cl->GetLabel(2) >= 0 ) labels.push_back(cl->GetLabel(2)) ;
669   }
670   std::sort( labels.begin(), labels.end() );      
671   labels.push_back( -1 ); // put -1 to the end    
672   int labelMax = -1, labelCur = -1, nLabelsMax = 0, nLabelsCurr = 0;
673   for ( unsigned int iLab = 0; iLab < labels.size(); iLab++ ) {
674     if ( labels[iLab] != labelCur ) {         
675       if ( labelCur >= 0 && nLabelsMax< nLabelsCurr ) {
676         nLabelsMax = nLabelsCurr;
677         labelMax = labelCur;
678       }
679       labelCur = labels[iLab];
680       nLabelsCurr = 0;
681     }
682     nLabelsCurr++;
683   }
684           
685   if( labelMax>=0 && nLabelsMax < wrong * nClusters ) labelMax = -labelMax;
686   
687   mcLabel = labelMax;
688                 
689   track->SetLabel( mcLabel );
690 }
691
692
693
694
695
696
697
698
699
700 //------------------------------------------------------------------------
701 void AliITStrackerHLT::UpdateESDtrack(AliESDtrack *tESD, AliHLTITSTrack* track, ULong_t flags) const
702 {
703   //
704   // Update ESD track
705   //
706   tESD->UpdateTrackParams(track,flags);
707   AliHLTITSTrack * oldtrack = (AliHLTITSTrack*)(tESD->GetITStrack());
708   if (oldtrack) delete oldtrack; 
709   tESD->SetITStrack(new AliHLTITSTrack(*track));
710 }
711
712
713
714
715 //------------------------------------------------------------------------
716 void AliITStrackerHLT::BuildMaterialLUT(TString material) {
717   //--------------------------------------------------------------------
718   // Fill a look-up table with mean material
719   //--------------------------------------------------------------------
720
721   Int_t n=1000;
722   Double_t mparam[7];
723   Double_t point1[3],point2[3];
724   Double_t phi,cosphi,sinphi,z;
725   // 0-5 layers, 6 pipe, 7-8 shields 
726   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
727   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
728
729   Int_t ifirst=0,ilast=0;  
730   if(material.Contains("Pipe")) {
731     ifirst=6; ilast=6;
732   } else if(material.Contains("Shields")) {
733     ifirst=7; ilast=8;
734   } else if(material.Contains("Layers")) {
735     ifirst=0; ilast=5;
736   } else {
737     Error("BuildMaterialLUT","Wrong layer name\n");
738   }
739
740   for(Int_t imat=ifirst; imat<=ilast; imat++) {
741     Double_t param[5]={0.,0.,0.,0.,0.};
742     for (Int_t i=0; i<n; i++) {
743       phi = 2.*TMath::Pi()*gRandom->Rndm();
744       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi); 
745       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
746       point1[0] = rmin[imat]*cosphi;
747       point1[1] = rmin[imat]*sinphi;
748       point1[2] = z;
749       point2[0] = rmax[imat]*cosphi;
750       point2[1] = rmax[imat]*sinphi;
751       point2[2] = z;
752       AliTracker::MeanMaterialBudget(point1,point2,mparam);
753       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
754     }
755     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
756     if(imat<=5) {
757       fxOverX0Layer[imat] = param[1];
758       fxTimesRhoLayer[imat] = param[0]*param[4];
759     } else if(imat==6) {
760       fxOverX0Pipe = param[1];
761       fxTimesRhoPipe = param[0]*param[4];
762     } else if(imat==7) {
763       fxOverX0Shield[0] = param[1];
764       fxTimesRhoShield[0] = param[0]*param[4];
765     } else if(imat==8) {
766       fxOverX0Shield[1] = param[1];
767       fxTimesRhoShield[1] = param[0]*param[4];
768     }
769   }
770   /*
771   printf("%s\n",material.Data());
772   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
773   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
774   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
775   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
776   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
777   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
778   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
779   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
780   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
781   */
782   return;
783 }
784
785
786
787
788 //------------------------------------------------------------------------
789 Int_t AliITStrackerHLT::CorrectForTPCtoITSDeadZoneMaterial(AliHLTITSTrack *t) {
790   //--------------------------------------------------------------------
791   // Correction for the material between the TPC and the ITS
792   //--------------------------------------------------------------------
793   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction 
794       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
795       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
796       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
797   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
798       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
799       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
800       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
801   } else {
802     printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
803     return 0;
804   }
805   
806   return 1;
807 }
808
809
810 //------------------------------------------------------------------------
811 Int_t AliITStrackerHLT::CorrectForPipeMaterial(AliHLTITSTrack *t,
812                                                bool InwardDirection) {
813   //-------------------------------------------------------------------
814   // Propagate beyond beam pipe and correct for material
815   // (material budget in different ways according to fUseTGeo value)
816   // Add time if going outward (PropagateTo or PropagateToTGeo)
817   //-------------------------------------------------------------------
818
819   // Define budget mode:
820   // 0: material from AliITSRecoParam (hard coded)
821   // 1: material from TGeo in one step (on the fly)
822   // 2: material from lut
823   // 3: material from TGeo in one step (same for all hypotheses)
824   Int_t mode;
825   switch(fUseTGeo) {
826   case 0:
827     mode=0; 
828     break;    
829   case 1:
830     mode=1;
831     break;    
832   case 2:
833     mode=2;
834     break;
835   case 3:
836     mode=3; 
837     break;
838   case 4:
839     mode=3;
840     break;
841   default:
842     mode=0;
843     break;
844   }
845   
846   Float_t  dir = (InwardDirection ? 1. : -1.);
847   Double_t rToGo= ( InwardDirection ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
848   Double_t xToGo, phi,z;
849
850   if (!t->GetLocalXPhiZat(rToGo,xToGo,phi,z)) return 0;
851
852   Double_t xOverX0,x0,lengthTimesMeanDensity;
853
854   switch(mode) {
855   case 0:
856     xOverX0 = AliITSRecoParam::GetdPipe();
857     x0 = AliITSRecoParam::GetX0Be();
858     lengthTimesMeanDensity = xOverX0*x0;
859     lengthTimesMeanDensity *= dir;
860     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
861     break;
862   case 1:
863     if (!t->PropagateToTGeo(xToGo,1)) return 0;
864     break;
865   case 2:
866     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");  
867     xOverX0 = fxOverX0Pipe;
868     lengthTimesMeanDensity = fxTimesRhoPipe;
869     lengthTimesMeanDensity *= dir;
870     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
871     break;
872   case 3:
873     double xOverX0PipeTrks, xTimesRhoPipeTrks;
874     if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
875     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
876                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
877     xOverX0PipeTrks = TMath::Abs(xOverX0)/angle;
878     xTimesRhoPipeTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
879     xOverX0 = xOverX0PipeTrks;
880     lengthTimesMeanDensity = xTimesRhoPipeTrks;
881     lengthTimesMeanDensity *= dir;
882     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
883     break;
884   }
885   
886   return 1;
887 }
888 //------------------------------------------------------------------------
889 Int_t AliITStrackerHLT::CorrectForShieldMaterial(AliHLTITSTrack *t,
890                                                  Int_t    shieldindex,
891                                                 bool InwardDirection) {
892   //-------------------------------------------------------------------
893   // Propagate beyond SPD or SDD shield and correct for material
894   // (material budget in different ways according to fUseTGeo value)
895   // Add time if going outward (PropagateTo or PropagateToTGeo)
896   //-------------------------------------------------------------------
897
898   // Define budget mode:
899   // 0: material from AliITSRecoParam (hard coded)
900   // 1: material from TGeo in steps of X cm (on the fly)
901   //    X = AliITSRecoParam::GetStepSizeTGeo()
902   // 2: material from lut
903   // 3: material from TGeo in one step (same for all hypotheses)
904   Int_t mode;
905   switch(fUseTGeo) {
906   case 0:
907     mode=0; 
908     break;    
909   case 1:
910     mode=1;
911     break;    
912   case 2:
913     mode=2;
914     break;
915   case 3:
916     mode=3;
917     break;
918   case 4:
919     mode=3;
920     break;
921   default:
922     mode=0;
923     break;
924   }
925
926
927   Float_t  dir = (InwardDirection ? 1. : -1.);
928   Double_t rToGo;
929
930   if (shieldindex==1 ) { // SDDouter
931     rToGo=(InwardDirection ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
932   } else if (shieldindex==0 ) {        // SPDouter
933     rToGo=(InwardDirection ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0)); 
934   } else {
935     Error("CorrectForShieldMaterial"," Wrong shield name\n");
936     return 0;
937   }
938   Double_t xToGo, phi,z;
939
940   if (!t->GetLocalXPhiZat(rToGo,xToGo,phi,z)) return 0;
941
942   Double_t xOverX0,x0,lengthTimesMeanDensity;
943   Int_t nsteps=1;
944
945   switch(mode) {
946   case 0:
947     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
948     x0 = AliITSRecoParam::GetX0shield(shieldindex);
949     lengthTimesMeanDensity = xOverX0*x0;
950     lengthTimesMeanDensity *= dir;
951     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
952     break;
953   case 1:
954     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/fRecoParam->GetStepSizeTGeo())+1;
955     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
956     break;
957   case 2:
958     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");  
959     xOverX0 = fxOverX0Shield[shieldindex];
960     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
961     lengthTimesMeanDensity *= dir;
962     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
963     break;
964   case 3:    
965     if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
966     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
967                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
968     double xOverX0ShieldTrks = TMath::Abs(xOverX0)/angle;
969     double xTimesRhoShieldTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
970     xOverX0 = xOverX0ShieldTrks;
971     lengthTimesMeanDensity = xTimesRhoShieldTrks;
972     lengthTimesMeanDensity *= dir;
973     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
974     break;
975   }
976
977   return 1;
978 }
979 //------------------------------------------------------------------------
980 Int_t AliITStrackerHLT::CorrectForLayerMaterial(AliHLTITSTrack *t,
981                                                Int_t layerindex,
982                                                bool InwardDirection ){
983   //-------------------------------------------------------------------
984   // Propagate beyond layer and correct for material
985   // (material budget in different ways according to fUseTGeo value)
986   // Add time if going outward (PropagateTo or PropagateToTGeo)
987   //-------------------------------------------------------------------
988
989   /*
990     Double_t r=fLayers[layerindex].GetR();
991     Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
992     Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
993     rToGo+= InwardDirection ?-deltar :deltar;
994     Double_t xToGo;
995     if (!t->GetLocalXat(rToGo,xToGo)) return 0;  
996   */
997   
998   if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
999
1000   Double_t lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
1001   if( !InwardDirection ) lengthTimesMeanDensity = -lengthTimesMeanDensity;
1002
1003   return t->CorrectForMeanMaterial(fxOverX0Layer[layerindex],lengthTimesMeanDensity,kTRUE);
1004 }
1005
1006
1007 //------------------------------------------------------------------------
1008 Bool_t AliITStrackerHLT::LocalModuleCoord(Int_t ilayer,Int_t idet,
1009                                        const AliHLTITSTrack *track,
1010                                        Float_t &xloc,Float_t &zloc) const {
1011   //-----------------------------------------------------------------
1012   // Gives position of track in local module ref. frame
1013   //-----------------------------------------------------------------
1014
1015   xloc=0.; 
1016   zloc=0.;
1017
1018   if(idet<0) return kFALSE;
1019
1020   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6 
1021
1022   Int_t lad = Int_t(idet/ndet) + 1;
1023
1024   Int_t det = idet - (lad-1)*ndet + 1;
1025
1026   Double_t xyzGlob[3],xyzLoc[3];
1027
1028   AliHLTITSDetector &detector = fLayers[ilayer].GetDetector(idet);
1029   // take into account the misalignment: xyz at real detector plane
1030   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
1031
1032   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
1033
1034   xloc = (Float_t)xyzLoc[0];
1035   zloc = (Float_t)xyzLoc[2];
1036
1037   return kTRUE;
1038 }
1039
1040 //------------------------------------------------------------------------
1041 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 {
1042   //--------------------------------------------------------------------
1043   // This function computes the rectangular road for this track
1044   //--------------------------------------------------------------------
1045
1046
1047   AliHLTITSDetector &det = fLayers[ilayer].GetDetector(idet);
1048   // take into account the misalignment: propagate track to misaligned detector plane
1049
1050   double y,z,snp,cov[3];
1051   if( !track->GetYZAtPhiX( det.GetPhi(),det.GetRmisal(), y, z, snp, cov))return 0;
1052
1053
1054   Double_t dz=fRecoParam->GetNSigmaRoadZ()*
1055                     TMath::Sqrt(cov[2] + 
1056                                 fRecoParam->GetNSigmaZLayerForRoadZ()*
1057                                 fRecoParam->GetNSigmaZLayerForRoadZ()*
1058                                 fRecoParam->GetSigmaZ2(ilayer));
1059   Double_t dy=fRecoParam->GetNSigmaRoadY()*
1060                     TMath::Sqrt(cov[0] + 
1061                                 fRecoParam->GetNSigmaYLayerForRoadY()*
1062                                 fRecoParam->GetNSigmaYLayerForRoadY()*
1063                                 fRecoParam->GetSigmaY2(ilayer));
1064       
1065   // track at boundary between detectors, enlarge road
1066   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1067   if ( (y-dy < det.GetYmin()+boundaryWidth) || 
1068        (y+dy > det.GetYmax()-boundaryWidth) || 
1069        (z-dz < det.GetZmin()+boundaryWidth) ||
1070        (z+dz > det.GetZmax()-boundaryWidth) ) {
1071     Float_t tgl = TMath::Abs(track->GetTgl());
1072     if (tgl > 1.) tgl=1.;
1073     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1074     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1075     if (snp > fRecoParam->GetMaxSnp()) return kFALSE;
1076     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1077   } // boundary
1078   
1079   // add to the road a term (up to 2-3 mm) to deal with misalignments
1080   dy = TMath::Sqrt(dy*dy + fRecoParam->GetRoadMisal()*fRecoParam->GetRoadMisal());
1081   dz = TMath::Sqrt(dz*dz + fRecoParam->GetRoadMisal()*fRecoParam->GetRoadMisal());
1082
1083   Double_t r = fLayers[ilayer].GetR();
1084   zmin = z - dz; 
1085   zmax = z + dz;
1086   ymin = y + r*det.GetPhi() - dy;
1087   ymax = y + r*det.GetPhi() + dy;
1088
1089   return kTRUE;
1090 }
1091
1092
1093
1094 Int_t AliITStrackerHLT::PropagateBack(AliESDEvent * /*event*/) 
1095 {
1096   // dummy
1097   return 0;
1098 }
1099
1100 Int_t AliITStrackerHLT::RefitInward(AliESDEvent * /*event*/ ) 
1101 {
1102   // dummy
1103   return 0;
1104 }
1105
1106
1107 Bool_t AliITStrackerHLT::GetTrackPoint(Int_t /*index*/, AliTrackPoint& /*p*/) const 
1108 {
1109   // dummy
1110   return 0;
1111 }
1112
1113 Bool_t AliITStrackerHLT::GetTrackPointTrackingError(Int_t /*index*/, 
1114                                                     AliTrackPoint& /*p*/, const AliESDtrack */*t*/) 
1115 {
1116   // dummy
1117   return 0;
1118 }