]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/ITS/tracking/AliITStrackerHLT.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
initialisation bug fixed
[u/mrichter/AliRoot.git] / HLT / ITS / tracking / AliITStrackerHLT.cxx
CommitLineData
2f399afc 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 "AliITSDetTypeRec.h"
52#include "AliITSRecPoint.h"
53#include "AliITSgeomTGeo.h"
54#include "AliITSReconstructor.h"
55#include "AliITSClusterParam.h"
56#include "AliITSsegmentation.h"
57#include "AliITSCalibration.h"
58#include "AliITSV0Finder.h"
59#include "AliITStrackerHLT.h"
60//#include "AliHLTTPCCATrackParam.h"
ca02aa66 61//#include "AliHLTVertexer.h"
2f399afc 62
63
133288b9 64ClassImp(AliITStrackerHLT)
2f399afc 65
6edb0fb5 66
67
2f399afc 68Bool_t AliITStrackerHLT::TransportToX( AliExternalTrackParam *t, double x ) const
69{
70 return t->PropagateTo( x, t->GetBz() );
71}
72
73Bool_t AliITStrackerHLT::TransportToPhiX( AliExternalTrackParam *t, double phi, double x ) const
74{
75 return t->Propagate( phi, x, t->GetBz() );
76}
77
78
79//------------------------------------------------------------------------
80Int_t AliITStrackerHLT::UpdateMI(AliHLTITSTrack* track, const AliITSRecPoint* cl,Double_t /*chi2*/,Int_t index) const
81{
82 //
83 // Update ITS track
84 //
85
86 if (cl->GetQ()<=0) return 0; // ingore the "virtual" clusters
87
88 Int_t layer = (index & 0xf0000000) >> 28;
89
90 // Take into account the mis-alignment (bring track to cluster plane)
91
92 Double_t xTrOrig=track->GetX();
93 if (!TransportToX( track, xTrOrig + cl->GetX() ) ) return 0;
94
95 Double_t err2Y, err2Z;
96
97 GetClusterErrors2( layer, cl, track, err2Y, err2Z );
98
99 Double_t p[2]={ cl->GetY(), cl->GetZ()};
100 Double_t cov[3]={err2Y, 0., err2Z};
101
102 Int_t updated = 1;
103 //if( layer!=2 && layer!=3 )
104 updated = track->AliExternalTrackParam::Update(p,cov);
105
106 int n = track->GetNumberOfClusters();
107 track->SetClusterIndex(n,index);
108 track->SetNumberOfClusters(n+1);
109
110 return updated;
111}
112
113
114AliHLTITSLayer AliITStrackerHLT::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
115
6edb0fb5 116AliITStrackerHLT::AliITStrackerHLT()
117 :AliTracker(),
87cc43e1 118 fRecoParam(0),
6edb0fb5 119 fEsd(0),
120 fUseTGeo(3),
121 fxOverX0Pipe(-1.),
122 fxTimesRhoPipe(-1.),
123 fDebugStreamer(0),
124 fITSChannelStatus(0),
125 fkDetTypeRec(0),
126 fTracks()
2f399afc 127{
128 //Default constructor
129 Int_t i;
130 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
131 for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
132 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
133}
134//------------------------------------------------------------------------
6edb0fb5 135AliITStrackerHLT::AliITStrackerHLT(const Char_t *geom)
136: AliTracker(),
87cc43e1 137 fRecoParam(0),
6edb0fb5 138 fEsd(0),
139 fUseTGeo(3),
140 fxOverX0Pipe(-1.),
141 fxTimesRhoPipe(-1.),
142 fDebugStreamer(0),
143 fITSChannelStatus(0),
144 fkDetTypeRec(0),
145 fTracks()
2f399afc 146{
147 //--------------------------------------------------------------------
148 //This is the AliITStrackerHLT constructor
149 //--------------------------------------------------------------------
150 if (geom) {
151 AliWarning("\"geom\" is actually a dummy argument !");
152 }
153
87cc43e1 154 fRecoParam = AliITSReconstructor::GetRecoParam();
155 if( !fRecoParam ) fRecoParam = AliITSRecoParam::GetLowFluxParam();
156
2f399afc 157 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
158 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
159 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
160
161 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
162 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
163 Double_t poff=TMath::ATan2(y,x);
164 Double_t zoff=z;
165 Double_t r=TMath::Sqrt(x*x + y*y);
166
167 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
168 r += TMath::Sqrt(x*x + y*y);
169 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
170 r += TMath::Sqrt(x*x + y*y);
171 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
172 r += TMath::Sqrt(x*x + y*y);
173 r*=0.25;
174
175 new (fgLayers+i-1) AliHLTITSLayer(r,poff,zoff,nlad,ndet);
176
177 for (Int_t j=1; j<nlad+1; j++) {
178 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
179 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
180 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
181 m.Multiply(tm);
182 Double_t txyz[3]={0.};
183 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
184 m.LocalToMaster(txyz,xyz);
185 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
186 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
187
188 if (phi<0) phi+=TMath::TwoPi();
189 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
190
191 AliHLTITSDetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
192 new(&det) AliHLTITSDetector(r,phi);
193 // compute the real radius (with misalignment)
194 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
195 mmisal.Multiply(tm);
196 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
197 mmisal.LocalToMaster(txyz,xyz);
198 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
199 det.SetRmisal(rmisal);
200
201 } // end loop on detectors
202 } // end loop on ladders
203 } // end loop on layers
204
205
87cc43e1 206 Double_t xyzVtx[]={ fRecoParam->GetXVdef(),
207 fRecoParam->GetYVdef(),
208 fRecoParam->GetZVdef()};
209 Double_t ersVtx[]={ fRecoParam->GetSigmaXVdef(),
210 fRecoParam->GetSigmaYVdef(),
211 fRecoParam->GetSigmaZVdef()};
2f399afc 212
2f399afc 213 SetVertex(xyzVtx,ersVtx);
214
215 // store positions of centre of SPD modules (in z)
216 Double_t tr[3];
217 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
218 fSPDdetzcentre[0] = tr[2];
219 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
220 fSPDdetzcentre[1] = tr[2];
221 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
222 fSPDdetzcentre[2] = tr[2];
223 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
224 fSPDdetzcentre[3] = tr[2];
225
87cc43e1 226 //fUseTGeo = fRecoParam->GetUseTGeoInTracker();
227 //if(fRecoParam->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
228 //AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
229 //fUseTGeo = 3;
230 //}
2f399afc 231
232 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
233 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
234
235 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
236
237}
238//------------------------------------------------------------------------
6edb0fb5 239AliITStrackerHLT::AliITStrackerHLT(const AliITStrackerHLT &tracker)
240:AliTracker(tracker),
87cc43e1 241 fRecoParam( tracker.fRecoParam),
6edb0fb5 242 fEsd(tracker.fEsd),
243 fUseTGeo(tracker.fUseTGeo),
244 fxOverX0Pipe(tracker.fxOverX0Pipe),
245 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
246 fDebugStreamer(tracker.fDebugStreamer),
247 fITSChannelStatus(tracker.fITSChannelStatus),
248 fkDetTypeRec(tracker.fkDetTypeRec),
249 fTracks()
2f399afc 250{
251 //Copy constructor
252 Int_t i;
253 for(i=0;i<4;i++) {
254 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
255 }
256 for(i=0;i<6;i++) {
257 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
258 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
259 }
260 for(i=0;i<2;i++) {
261 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
262 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
263 }
264}
265//------------------------------------------------------------------------
266AliITStrackerHLT & AliITStrackerHLT::operator=(const AliITStrackerHLT &tracker){
267 //Assignment operator
268 this->~AliITStrackerHLT();
269 new(this) AliITStrackerHLT(tracker);
270 return *this;
271}
272//------------------------------------------------------------------------
273AliITStrackerHLT::~AliITStrackerHLT()
274{
275 //
276 //destructor
277 //
278 if (fDebugStreamer) {
279 //fDebugStreamer->Close();
280 delete fDebugStreamer;
281 }
282 if(fITSChannelStatus) delete fITSChannelStatus;
283}
284
285
286//------------------------------------------------------------------------
287void AliITStrackerHLT::ReadBadFromDetTypeRec() {
288 //--------------------------------------------------------------------
289 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
290 //i.e. from OCDB
291 //--------------------------------------------------------------------
292
87cc43e1 293 if(!fRecoParam->GetUseBadZonesFromOCDB()) return;
2f399afc 294
295 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
296
297 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
298
299 // ITS channels map
300 if(fITSChannelStatus) delete fITSChannelStatus;
301 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
302
303 // ITS detectors and chips
304 Int_t i=0,j=0,k=0,ndet=0;
305 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
306 Int_t nBadDetsPerLayer=0;
307 ndet=AliITSgeomTGeo::GetNDetectors(i);
308 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
309 for (k=1; k<ndet+1; k++) {
310 AliHLTITSDetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
311 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
312 if(det.IsBad()) {nBadDetsPerLayer++;}
313 } // end loop on detectors
314 } // end loop on ladders
315 Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
316 } // end loop on layers
317
318 return;
319}
6edb0fb5 320
321
322//------------------------------------------------------------------------
323void AliITStrackerHLT::LoadClusters( std::vector<AliITSRecPoint> clusters)
324{
325 //--------------------------------------------------------------------
326 //This function loads ITS clusters
327 //--------------------------------------------------------------------
87cc43e1 328
6edb0fb5 329 //SignDeltas(clusters,GetZ());
87cc43e1 330 //std::cout<<"CA ITS: NClusters "<<clusters.size()<<std::endl;
6edb0fb5 331 for( unsigned int icl=0; icl<clusters.size(); icl++ ){
87cc43e1 332 std::cout<<"CA ITS: icl "<<icl<<std::endl;
333
6edb0fb5 334 AliITSRecPoint &cl = clusters[icl];
87cc43e1 335 //std::cout<<"CA ITS: icl "<<icl<<": "
336 //<<cl.GetX()<<" "<<cl.GetY()<<" "<<cl.GetZ()<<" "<<cl.GetDetectorIndex()<<" "<<cl.GetLayer()
337 //<<" "<<cl.GetVolumeId()
338 //<<std::endl;
339
340
6edb0fb5 341 Int_t i=cl.GetLayer();
342 //Int_t ndet=fgLayers[i].GetNdetectors();
87cc43e1 343 //Int_t detector=cl.GetDetectorIndex();
344 if (!cl.Misalign()) AliWarning("Can't misalign this cluster !"); //SG!!!
345 fgLayers[i].InsertCluster(new AliITSRecPoint(cl));
6edb0fb5 346 }
347
348 for( int i=0; i<AliITSgeomTGeo::GetNLayers(); i++ ){
349 for( int detector = 0; detector<fgLayers[i].GetNdetectors(); detector++ ){
6edb0fb5 350 // add dead zone "virtual" cluster in SPD, if there is a cluster within
351 // zwindow cm from the dead zone
87cc43e1 352
353 fRecoParam->GetAddVirtualClustersInDeadZone();
354
355 if (i<2 && fRecoParam->GetAddVirtualClustersInDeadZone()) {
356 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*fRecoParam->GetXPassDeadZoneHits()) {
6edb0fb5 357 Int_t lab[4] = {0,0,0,detector};
358 Int_t info[3] = {0,0,i};
359 Float_t q = 0.; // this identifies virtual clusters
360 Float_t hit[5] = {xdead,
361 0.,
87cc43e1 362 fRecoParam->GetSigmaXDeadZoneHit2(),
363 fRecoParam->GetSigmaZDeadZoneHit2(),
6edb0fb5 364 q};
365 Bool_t local = kTRUE;
87cc43e1 366 Double_t zwindow = fRecoParam->GetZWindowDeadZone();
6edb0fb5 367 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
368 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
369 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
370 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
371 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
372 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
373 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
374 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
375 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
376 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
377 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
378 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
379 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
380 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
381 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
382 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
383 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
384 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
385 }
386 } // "virtual" clusters in SPD
387 }
388 fgLayers[i].ResetRoad(); //road defined by the cluster density
389 fgLayers[i].SortClusters();
390 }
391}
392
393
2f399afc 394//------------------------------------------------------------------------
395Int_t AliITStrackerHLT::LoadClusters(TTree *cTree) {
396 //--------------------------------------------------------------------
397 //This function loads ITS clusters
398 //--------------------------------------------------------------------
399 TBranch *branch=cTree->GetBranch("ITSRecPoints");
400 if (!branch) {
401 Error("LoadClusters"," can't get the branch !\n");
402 return 1;
403 }
404
405 static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
406 branch->SetAddress(&clusters);
407
408 Int_t i=0,j=0,ndet=0;
409 Int_t detector=0;
410 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
411 ndet=fgLayers[i].GetNdetectors();
412 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
413 for (; j<jmax; j++) {
414 if (!cTree->GetEvent(j)) continue;
415 Int_t ncl=clusters->GetEntriesFast();
416 SignDeltas(clusters,GetZ());
417
418 while (ncl--) {
419 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
420 detector=c->GetDetectorIndex();
421
422 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
423
424 fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
425 }
426 clusters->Delete();
427 // add dead zone "virtual" cluster in SPD, if there is a cluster within
428 // zwindow cm from the dead zone
87cc43e1 429 if (i<2 && fRecoParam->GetAddVirtualClustersInDeadZone()) {
430 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*fRecoParam->GetXPassDeadZoneHits()) {
2f399afc 431 Int_t lab[4] = {0,0,0,detector};
432 Int_t info[3] = {0,0,i};
433 Float_t q = 0.; // this identifies virtual clusters
434 Float_t hit[5] = {xdead,
435 0.,
87cc43e1 436 fRecoParam->GetSigmaXDeadZoneHit2(),
437 fRecoParam->GetSigmaZDeadZoneHit2(),
2f399afc 438 q};
439 Bool_t local = kTRUE;
87cc43e1 440 Double_t zwindow = fRecoParam->GetZWindowDeadZone();
2f399afc 441 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
442 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
443 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
444 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
445 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
446 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
447 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
448 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
449 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
450 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
451 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
452 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
453 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
454 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
455 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
456 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
457 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
458 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
459 }
460 } // "virtual" clusters in SPD
461
462 }
463 //
464 fgLayers[i].ResetRoad(); //road defined by the cluster density
465 fgLayers[i].SortClusters();
466 }
467
468 dummy.Clear();
469
470 return 0;
471}
472//------------------------------------------------------------------------
473void AliITStrackerHLT::UnloadClusters() {
474 //--------------------------------------------------------------------
475 //This function unloads ITS clusters
476 //--------------------------------------------------------------------
477 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
478}
479
480
481
482//------------------------------------------------------------------------
483Int_t AliITStrackerHLT::CorrectForTPCtoITSDeadZoneMaterial(AliHLTITSTrack *t) {
484 //--------------------------------------------------------------------
485 // Correction for the material between the TPC and the ITS
486 //--------------------------------------------------------------------
487 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
488 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
489 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
490 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
491 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
492 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
493 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
494 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
495 } else {
496 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
497 return 0;
498 }
499
500 return 1;
501}
502
503#include "TStopwatch.h"
504
6edb0fb5 505void AliITStrackerHLT::Reconstruct( std::vector<AliExternalTrackParam> tracksTPC )
506{
507 //--------------------------------------------------------------------
508 // This functions reconstructs ITS tracks
509 // The clusters must be already loaded !
510 //--------------------------------------------------------------------
511
512 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
513 fTracks.clear();
514
515 for( unsigned int itr=0; itr<tracksTPC.size(); itr++ ){
516
517 AliHLTITSTrack tMI( tracksTPC[itr] );
518 AliHLTITSTrack *t = &tMI;
519 t->SetTPCtrackId( itr );
520 t->SetMass(pimass);
521 t->SetExpQ(0);
522
523 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) continue;
524
525 //Int_t tpcLabel=t->GetLabel(); //save the TPC track label
526
527 FollowProlongationTree(t);
528 int nclu=0;
529 for(Int_t i=0; i<6; i++) {
530 if( t->GetClusterIndex(i)>=0 ) nclu++;
531 }
532 //cout<<"N assigned ITS clusters = "<<nclu<<std::endl;
533 if( nclu>0 ){
534 t->SetLabel(-1);//tpcLabel);
535 t->SetFakeRatio(1.);
536 CookLabel(t,0.); //For comparison only
537 //cout<<"label = "<<t->GetLabel()<<" / "<<tpcLabel<<endl;
538 TransportToX(t, 0 );
539 //cout<<"\n fill track : parameters at "<<t->GetX()<<": "<< TMath::Sqrt(TMath::Abs(t->GetSigmaY2()))<<" "<< TMath::Sqrt(TMath::Abs(t->GetSigmaY2()))<<endl;
540 //t->Print();
541 fTracks.push_back( *t );
542 }
543 }
544
545 //AliHLTVertexer vertexer;
546 //vertexer.SetESD( event );
547 //vertexer.FindPrimaryVertex();
548 //vertexer.FindV0s();
549}
550
551
552
2f399afc 553//------------------------------------------------------------------------
554Int_t AliITStrackerHLT::Clusters2Tracks(AliESDEvent *event) {
555 //--------------------------------------------------------------------
556 // This functions reconstructs ITS tracks
557 // The clusters must be already loaded !
558 //--------------------------------------------------------------------
87cc43e1 559
2f399afc 560 TStopwatch timer;
561
562 fEsd = event; // store pointer to the esd
6edb0fb5 563 std::vector<AliExternalTrackParam> tracksTPC;
2f399afc 564
6edb0fb5 565 for( int itr=0; itr<event->GetNumberOfTracks(); itr++ ){
2f399afc 566
6edb0fb5 567 AliESDtrack *esdTrack = event->GetTrack(itr);
2f399afc 568
6edb0fb5 569 if ((esdTrack->GetStatus()&AliESDtrack::kTPCin)==0) continue;
570 if (esdTrack->GetStatus()&AliESDtrack::kTPCout) continue;
571 if (esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
572 if (esdTrack->GetKinkIndex(0)>0) continue; //kink daughter
573
574 AliHLTITSTrack t(*esdTrack);
575 t.SetTPCtrackId( itr );
576 tracksTPC.push_back( t );
577 }
2f399afc 578
6edb0fb5 579 Reconstruct( tracksTPC );
580
581 for( unsigned int itr=0; itr<fTracks.size(); itr++ ){
582 AliHLTITSTrack &t = fTracks[itr];
583 UpdateESDtrack(event->GetTrack(t.TPCtrackId()), &t, AliESDtrack::kITSin);
584 }
585
ca02aa66 586 //AliHLTVertexer vertexer;
587 //vertexer.SetESD( event );
588 //vertexer.FindPrimaryVertex();
589 //vertexer.FindV0s();
2f399afc 590
591 timer.Stop();
592 static double totalTime = 0;
593 static int nEvnts = 0;
594 totalTime+=timer.CpuTime();
595 nEvnts++;
596 std::cout<<"\n\n ITS tracker time = "<<totalTime/nEvnts<<" [s/ev] for "<<nEvnts<<" events\n\n "<<std::endl;
597 return 0;
598}
599
600
601//------------------------------------------------------------------------
602Int_t AliITStrackerHLT::PropagateBack(AliESDEvent * /*event*/) {
603 return 0;
604}
605
606//------------------------------------------------------------------------
607Int_t AliITStrackerHLT::RefitInward(AliESDEvent * /*event*/ ) {
608 return 0;
609}
610//------------------------------------------------------------------------
611AliCluster *AliITStrackerHLT::GetCluster(Int_t index) const {
612 //--------------------------------------------------------------------
613 // Return pointer to a given cluster
614 //--------------------------------------------------------------------
615 Int_t l=(index & 0xf0000000) >> 28;
616 Int_t c=(index & 0x0fffffff) >> 00;
617 return fgLayers[l].GetCluster(c);
618}
619//------------------------------------------------------------------------
620Bool_t AliITStrackerHLT::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
621 //--------------------------------------------------------------------
622 // Get track space point with index i
623 //--------------------------------------------------------------------
624
625 Int_t l=(index & 0xf0000000) >> 28;
626 Int_t c=(index & 0x0fffffff) >> 00;
627 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
628 Int_t idet = cl->GetDetectorIndex();
629
630 Float_t xyz[3];
631 Float_t cov[6];
632 cl->GetGlobalXYZ(xyz);
633 cl->GetGlobalCov(cov);
634 p.SetXYZ(xyz, cov);
635 p.SetCharge(cl->GetQ());
636 p.SetDriftTime(cl->GetDriftTime());
637 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
638 switch (l) {
639 case 0:
640 iLayer = AliGeomManager::kSPD1;
641 break;
642 case 1:
643 iLayer = AliGeomManager::kSPD2;
644 break;
645 case 2:
646 iLayer = AliGeomManager::kSDD1;
647 break;
648 case 3:
649 iLayer = AliGeomManager::kSDD2;
650 break;
651 case 4:
652 iLayer = AliGeomManager::kSSD1;
653 break;
654 case 5:
655 iLayer = AliGeomManager::kSSD2;
656 break;
657 default:
658 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
659 break;
660 };
661 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
662 p.SetVolumeID((UShort_t)volid);
663 return kTRUE;
664}
665//------------------------------------------------------------------------
666Bool_t AliITStrackerHLT::GetTrackPointTrackingError(Int_t index,
667 AliTrackPoint& p, const AliESDtrack *t) {
668 //--------------------------------------------------------------------
669 // Get track space point with index i
670 // (assign error estimated during the tracking)
671 //--------------------------------------------------------------------
672
673 Int_t l=(index & 0xf0000000) >> 28;
674 Int_t c=(index & 0x0fffffff) >> 00;
675 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
676 Int_t idet = cl->GetDetectorIndex();
677
678 const AliHLTITSDetector &det=fgLayers[l].GetDetector(idet);
679
680 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
681 Float_t detxy[2];
682 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
683 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
684 Double_t alpha = t->GetAlpha();
685 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
686 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
687 phi += alpha-det.GetPhi();
688 Float_t tgphi = TMath::Tan(phi);
689
690 Float_t tgl = t->GetTgl(); // tgl about const along track
691 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
692
693 Float_t errlocalx,errlocalz;
694 Bool_t addMisalErr=kFALSE;
695 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errlocalx,errlocalz,addMisalErr);
696
697 Float_t xyz[3];
698 Float_t cov[6];
699 cl->GetGlobalXYZ(xyz);
700 // cl->GetGlobalCov(cov);
701 Float_t pos[3] = {0.,0.,0.};
702 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errlocalx*errlocalx,errlocalz*errlocalz,0);
703 tmpcl.GetGlobalCov(cov);
704
705 p.SetXYZ(xyz, cov);
706 p.SetCharge(cl->GetQ());
707 p.SetDriftTime(cl->GetDriftTime());
708
709 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
710 switch (l) {
711 case 0:
712 iLayer = AliGeomManager::kSPD1;
713 break;
714 case 1:
715 iLayer = AliGeomManager::kSPD2;
716 break;
717 case 2:
718 iLayer = AliGeomManager::kSDD1;
719 break;
720 case 3:
721 iLayer = AliGeomManager::kSDD2;
722 break;
723 case 4:
724 iLayer = AliGeomManager::kSSD1;
725 break;
726 case 5:
727 iLayer = AliGeomManager::kSSD2;
728 break;
729 default:
730 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
731 break;
732 };
733 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
734
735 p.SetVolumeID((UShort_t)volid);
736 return kTRUE;
737}
738
739
740//------------------------------------------------------------------------
741void AliITStrackerHLT::FollowProlongationTree(AliHLTITSTrack * track )
742{
743 //cout<<endl;
744 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
745 //cout<<"\nLayer "<<ilayer<<endl;
746
747 AliHLTITSLayer &layer=fgLayers[ilayer];
748
749 //cout<<" shield material.. "<<endl;
750
751 // material between SSD and SDD, SDD and SPD
752 if (ilayer==3 && !CorrectForShieldMaterial(track,"SDD","inward")) continue;
753 if (ilayer==1 && !CorrectForShieldMaterial(track,"SPD","inward")) continue;
754
755 int idet;
756
757 {
758 // propagate to the layer radius
759
760 Double_t r = layer.GetR(), phi,z;
761 Double_t xToGo;
762 //cout<<" propagate to layer R= "<<r<<" .."<<endl;
763 if( !track->GetLocalXat(r,xToGo) ) continue;
764 if( !TransportToX(track, xToGo) ) continue;
765
766 // detector number
767
768 if (!track->GetPhiZat(r,phi,z)) continue;
769 idet=layer.FindDetectorIndex(phi,z);
770 //cout<<" detector number = "<<idet<<endl;
771 }
772
773
774 //cout<<" correct for the layer material .. "<<endl;
775
776 // correct for the layer material
777 {
778 Double_t trackGlobXYZ1[3];
779 if (!track->GetXYZ(trackGlobXYZ1)) continue;
780 CorrectForLayerMaterial(track,ilayer,trackGlobXYZ1,"inward");
781 }
782
783 // track outside layer acceptance in z
784
785 if( idet<0 ) continue;
786
787 // propagate to the intersection with the detector plane
788
789 //cout<<" propagate to the intersection with the detector .. "<<endl;
790
791 const AliHLTITSDetector &det=layer.GetDetector( idet );
792 if (!TransportToPhiX( track, det.GetPhi(), det.GetR() ) ) continue;
793
794 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
795
796 // road in global (rphi,z) [i.e. in tracking ref. system]
797
798 Double_t zmin,zmax,ymin,ymax;
799
800 //cout<<" ComputeRoad .. "<<endl;
801 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
802
803 //cout<<" Road: y["<<ymin<<","<<ymax<<"], z["<<zmin<<","<<zmax<<"] "<<endl;
804
805 // select clusters in road
806
807 //cout<<" SelectClusters .. "<<endl;
808 layer.SelectClusters(zmin,zmax,ymin,ymax);
809
810 // Define criteria for track-cluster association
811
812 Double_t msz = track->GetSigmaZ2() +
87cc43e1 813 fRecoParam->GetNSigmaZLayerForRoadZ()*
814 fRecoParam->GetNSigmaZLayerForRoadZ()*
815 fRecoParam->GetSigmaZ2(ilayer);
2f399afc 816
817 Double_t msy = track->GetSigmaY2() +
87cc43e1 818 fRecoParam->GetNSigmaYLayerForRoadY()*
819 fRecoParam->GetNSigmaYLayerForRoadY()*
820 fRecoParam->GetSigmaY2(ilayer);
2f399afc 821
87cc43e1 822 msz *= fRecoParam->GetNSigma2RoadZNonC();
823 msy *= fRecoParam->GetNSigma2RoadYNonC();
2f399afc 824
825 msz = 1./msz; // 1/RoadZ^2
826 msy = 1./msy; // 1/RoadY^2
827
828 //
829 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
830 //
831
832 const AliITSRecPoint *cl=0;
833 Int_t clidx=-1;
87cc43e1 834 Double_t chi2trkcl=fRecoParam->GetMaxChi2(); // init with big value
2f399afc 835 Bool_t deadzoneSPD=kFALSE;
836
837 // check if the road contains a dead zone
838 Bool_t noClusters = !layer.GetNextCluster(clidx,kTRUE);
839
840 Double_t dz=0.5*(zmax-zmin);
841 Double_t dy=0.5*(ymax-ymin);
842
843 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,noClusters);
844
845 // create a prolongation without clusters (check also if there are no clusters in the road)
846
847 if (dead==1) { // dead zone at z=0,+-7cm in SPD
848 deadzoneSPD=kTRUE;
849 }
850
851 clidx=-1;
852
853 //cout<<" loop over clusters in the road .. "<<endl;
854
855 // loop over clusters in the road
856 const AliITSRecPoint *bestCluster=0;
857 double bestChi2 = 1.e10;
858 AliHLTITSTrack bestTrack( *track );
859 int bestIdx = -1;
860 while ((cl=layer.GetNextCluster(clidx))!=0) {
861 //cout<<" cluster: "<<cl->GetX()<<" "<<cl->GetY()<<" "<<cl->GetZ()<<endl;
862 AliHLTITSTrack t(*track);
863 if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
864
865 Int_t idetc=cl->GetDetectorIndex();
866
867 //cout<<" cluster detector: "<<idetc<<endl;
868
869 if ( idet !=idetc ) { // new cluster's detector
870 const AliHLTITSDetector &detc=layer.GetDetector(idetc);
871 if (!TransportToPhiX( track, detc.GetPhi(),detc.GetR()) ) continue;
872 t = *track;
873 idet = idetc;
874 }
875
876 // take into account misalignment (bring track to real detector plane)
877
878 if (!TransportToX( &t, t.GetX() + cl->GetX() ) ) continue;
879 double chi2 = ( (t.GetZ()-cl->GetZ())*(t.GetZ()-cl->GetZ())*msz +
880 (t.GetY()-cl->GetY())*(t.GetY()-cl->GetY())*msy );
881 //cout<<" chi2="<<chi2<<endl;
882 if ( chi2 < bestChi2 ){
883 bestChi2 = chi2;
884 bestCluster = cl;
885 bestTrack = t;
886 bestIdx = clidx;
887 continue;
888 }
889 }
890
891 //cout<<" best chi2= "<<bestChi2<<endl;
892
893 if( bestCluster && bestChi2 <=1. ){
894
895 //cout<<" cluster found "<<endl;
896 *track = bestTrack;
897
898
899 // calculate track-clusters chi2
900 chi2trkcl = GetPredictedChi2MI(track,bestCluster,ilayer);
901 //cout<<" track-clusters chi2 = "<<chi2trkcl<<endl;
87cc43e1 902 //cout<<" max chi2 = "<<fRecoParam->GetMaxChi2s(ilayer)<<endl;
2f399afc 903
904 // chi2 cut
87cc43e1 905 if (chi2trkcl < fRecoParam->GetMaxChi2s(ilayer)) {
2f399afc 906 if (bestCluster->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
907 //cout<<"set index.."<<endl;
908 //cout<<"set index ok"<<endl;
909 if (bestCluster->GetQ()!=0) { // real cluster
910 //cout<<" UpdateMI ... "<<endl;
911 if (!UpdateMI(track,bestCluster,chi2trkcl,(ilayer<<28)+bestIdx)) {
912 continue;
913 }
914 }
915 }
916 }
917 //cout<<" goto next layer "<<endl;
918
919 }
920}
921
922
923//------------------------------------------------------------------------
924AliHLTITSLayer & AliITStrackerHLT::GetLayer(Int_t layer) const
925{
926 //--------------------------------------------------------------------
927 //
928 //
929 return fgLayers[layer];
930}
931
932
933
2f399afc 934//------------------------------------------------------------------------
6edb0fb5 935void AliITStrackerHLT::CookLabel(AliHLTITSTrack *track,Float_t wrong) const
936{
937 // get MC label for the track
938
939 Int_t mcLabel = -1;
940
941 vector<int> labels;
942 Int_t nClusters = track->GetNumberOfClusters();
943
944 for (Int_t i=0; i<nClusters; i++){
945 Int_t cindex = track->GetClusterIndex(i);
946 //Int_t l=(cindex & 0xf0000000) >> 28;
947 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
948 if ( cl->GetLabel(0) >= 0 ) labels.push_back(cl->GetLabel(0)) ;
949 if ( cl->GetLabel(1) >= 0 ) labels.push_back(cl->GetLabel(1)) ;
950 if ( cl->GetLabel(2) >= 0 ) labels.push_back(cl->GetLabel(2)) ;
951 }
952 std::sort( labels.begin(), labels.end() );
953 labels.push_back( -1 ); // put -1 to the end
954 int labelMax = -1, labelCur = -1, nLabelsMax = 0, nLabelsCurr = 0;
955 for ( unsigned int iLab = 0; iLab < labels.size(); iLab++ ) {
956 if ( labels[iLab] != labelCur ) {
957 if ( labelCur >= 0 && nLabelsMax< nLabelsCurr ) {
958 nLabelsMax = nLabelsCurr;
959 labelMax = labelCur;
960 }
961 labelCur = labels[iLab];
962 nLabelsCurr = 0;
963 }
964 nLabelsCurr++;
965 }
966
967 if( labelMax>=0 && nLabelsMax <= wrong * nClusters ) labelMax = -labelMax;
968
969 mcLabel = labelMax;
970
971 track->SetLabel( mcLabel );
2f399afc 972}
973
974
975
976void AliITStrackerHLT::GetClusterErrors2( Int_t layer, const AliITSRecPoint *cluster, AliHLTITSTrack* track, double &err2Y, double &err2Z ) const
977{
978 Float_t erry,errz;
979 Float_t theta = track->GetTgl();
980 Float_t phi = track->GetSnp();
981 phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
982 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
983 err2Y = erry*erry;
984 err2Z = errz*errz;
985}
986
987
988//------------------------------------------------------------------------
989Double_t AliITStrackerHLT::GetPredictedChi2MI(AliHLTITSTrack* track, const AliITSRecPoint *cluster,Int_t layer)
990{
991 //
992 // Compute predicted chi2
993 //
994
995 AliHLTITSTrack t(*track);
996 if (!t.Propagate( t.GetX()+cluster->GetX())) return 1000.;
997
998 Double_t err2Y,err2Z;
999 GetClusterErrors2( layer, cluster, &t, err2Y, err2Z );
1000
1001 Double_t chi2 = t.GetPredictedChi2(cluster->GetY(),cluster->GetZ(),err2Y,err2Z);
1002
1003 Float_t ny,nz;
1004 Float_t theta = track->GetTgl();
1005 Float_t phi = track->GetSnp();
1006 phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
1007 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
1008 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
1009 if (delta>1){
1010 chi2+=0.5*TMath::Min(delta/2,2.);
1011 chi2+=2.*cluster->GetDeltaProbability();
1012 }
1013 return chi2;
1014}
1015
1016
6edb0fb5 1017
1018
2f399afc 1019//------------------------------------------------------------------------
1020void AliITStrackerHLT::SignDeltas(const TObjArray *clusterArray, Float_t vz)
1021{
1022 //
1023 // Clusters from delta electrons?
1024 //
1025 Int_t entries = clusterArray->GetEntriesFast();
1026 if (entries<4) return;
1027 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
1028 Int_t layer = cluster->GetLayer();
1029 if (layer>1) return;
1030 Int_t index[10000];
1031 Int_t ncandidates=0;
1032 Float_t r = (layer>0)? 7:4;
1033 //
1034 for (Int_t i=0;i<entries;i++){
1035 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
1036 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
1037 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
1038 index[ncandidates] = i; //candidate to belong to delta electron track
1039 ncandidates++;
1040 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
1041 cl0->SetDeltaProbability(1);
1042 }
1043 }
1044 //
1045 //
1046 //
1047 for (Int_t i=0;i<ncandidates;i++){
1048 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
1049 if (cl0->GetDeltaProbability()>0.8) continue;
1050 //
1051 Int_t ncl = 0;
1052 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
1053 sumy=sumz=sumy2=sumyz=sumw=0.0;
1054 for (Int_t j=0;j<ncandidates;j++){
1055 if (i==j) continue;
1056 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
1057 //
1058 Float_t dz = cl0->GetZ()-cl1->GetZ();
1059 Float_t dy = cl0->GetY()-cl1->GetY();
1060 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
1061 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
1062 y[ncl] = cl1->GetY();
1063 z[ncl] = cl1->GetZ();
1064 sumy+= y[ncl]*weight;
1065 sumz+= z[ncl]*weight;
1066 sumy2+=y[ncl]*y[ncl]*weight;
1067 sumyz+=y[ncl]*z[ncl]*weight;
1068 sumw+=weight;
1069 ncl++;
1070 }
1071 }
1072 if (ncl<4) continue;
1073 Float_t det = sumw*sumy2 - sumy*sumy;
1074 Float_t delta=1000;
1075 if (TMath::Abs(det)>0.01){
1076 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
1077 Float_t k = (sumyz*sumw - sumy*sumz)/det;
1078 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
1079 }
1080 else{
1081 Float_t z0 = sumyz/sumy;
1082 delta = TMath::Abs(cl0->GetZ()-z0);
1083 }
1084 if ( delta<0.05) {
1085 cl0->SetDeltaProbability(1-20.*delta);
1086 }
1087 }
1088}
1089//------------------------------------------------------------------------
6edb0fb5 1090void AliITStrackerHLT::UpdateESDtrack(AliESDtrack *tESD, AliHLTITSTrack* track, ULong_t flags) const
2f399afc 1091{
1092 //
1093 // Update ESD track
1094 //
6edb0fb5 1095 tESD->UpdateTrackParams(track,flags);
1096 AliHLTITSTrack * oldtrack = (AliHLTITSTrack*)(tESD->GetITStrack());
2f399afc 1097 if (oldtrack) delete oldtrack;
6edb0fb5 1098 tESD->SetITStrack(new AliHLTITSTrack(*track));
2f399afc 1099}
1100
1101
1102
1103
1104//------------------------------------------------------------------------
1105void AliITStrackerHLT::BuildMaterialLUT(TString material) {
1106 //--------------------------------------------------------------------
1107 // Fill a look-up table with mean material
1108 //--------------------------------------------------------------------
1109
1110 Int_t n=1000;
1111 Double_t mparam[7];
1112 Double_t point1[3],point2[3];
1113 Double_t phi,cosphi,sinphi,z;
1114 // 0-5 layers, 6 pipe, 7-8 shields
1115 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
1116 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
1117
1118 Int_t ifirst=0,ilast=0;
1119 if(material.Contains("Pipe")) {
1120 ifirst=6; ilast=6;
1121 } else if(material.Contains("Shields")) {
1122 ifirst=7; ilast=8;
1123 } else if(material.Contains("Layers")) {
1124 ifirst=0; ilast=5;
1125 } else {
1126 Error("BuildMaterialLUT","Wrong layer name\n");
1127 }
1128
1129 for(Int_t imat=ifirst; imat<=ilast; imat++) {
1130 Double_t param[5]={0.,0.,0.,0.,0.};
1131 for (Int_t i=0; i<n; i++) {
1132 phi = 2.*TMath::Pi()*gRandom->Rndm();
1133 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
1134 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
1135 point1[0] = rmin[imat]*cosphi;
1136 point1[1] = rmin[imat]*sinphi;
1137 point1[2] = z;
1138 point2[0] = rmax[imat]*cosphi;
1139 point2[1] = rmax[imat]*sinphi;
1140 point2[2] = z;
1141 AliTracker::MeanMaterialBudget(point1,point2,mparam);
1142 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
1143 }
1144 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
1145 if(imat<=5) {
1146 fxOverX0Layer[imat] = param[1];
1147 fxTimesRhoLayer[imat] = param[0]*param[4];
1148 } else if(imat==6) {
1149 fxOverX0Pipe = param[1];
1150 fxTimesRhoPipe = param[0]*param[4];
1151 } else if(imat==7) {
1152 fxOverX0Shield[0] = param[1];
1153 fxTimesRhoShield[0] = param[0]*param[4];
1154 } else if(imat==8) {
1155 fxOverX0Shield[1] = param[1];
1156 fxTimesRhoShield[1] = param[0]*param[4];
1157 }
1158 }
1159 /*
1160 printf("%s\n",material.Data());
1161 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
1162 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
1163 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
1164 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
1165 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
1166 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
1167 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
1168 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
1169 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
1170 */
1171 return;
1172}
1173
1174//------------------------------------------------------------------------
1175Int_t AliITStrackerHLT::CorrectForPipeMaterial(AliHLTITSTrack *t,
1176 TString direction) {
1177 //-------------------------------------------------------------------
1178 // Propagate beyond beam pipe and correct for material
1179 // (material budget in different ways according to fUseTGeo value)
1180 // Add time if going outward (PropagateTo or PropagateToTGeo)
1181 //-------------------------------------------------------------------
1182
1183 // Define budget mode:
1184 // 0: material from AliITSRecoParam (hard coded)
1185 // 1: material from TGeo in one step (on the fly)
1186 // 2: material from lut
1187 // 3: material from TGeo in one step (same for all hypotheses)
1188 Int_t mode;
1189 switch(fUseTGeo) {
1190 case 0:
1191 mode=0;
1192 break;
1193 case 1:
1194 mode=1;
1195 break;
1196 case 2:
1197 mode=2;
1198 break;
1199 case 3:
1200 mode=3;
1201 break;
1202 case 4:
1203 mode=3;
1204 break;
1205 default:
1206 mode=0;
1207 break;
1208 }
1209
1210 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
1211 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
1212 Double_t xToGo;
1213 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1214
1215 Double_t xOverX0,x0,lengthTimesMeanDensity;
1216
1217 switch(mode) {
1218 case 0:
1219 xOverX0 = AliITSRecoParam::GetdPipe();
1220 x0 = AliITSRecoParam::GetX0Be();
1221 lengthTimesMeanDensity = xOverX0*x0;
1222 lengthTimesMeanDensity *= dir;
1223 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1224 break;
1225 case 1:
1226 if (!t->PropagateToTGeo(xToGo,1)) return 0;
1227 break;
1228 case 2:
1229 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
1230 xOverX0 = fxOverX0Pipe;
1231 lengthTimesMeanDensity = fxTimesRhoPipe;
1232 lengthTimesMeanDensity *= dir;
1233 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1234 break;
1235 case 3:
1236 double xOverX0PipeTrks, xTimesRhoPipeTrks;
1237 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
1238 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1239 ((1.-t->GetSnp())*(1.+t->GetSnp())));
1240 xOverX0PipeTrks = TMath::Abs(xOverX0)/angle;
1241 xTimesRhoPipeTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
1242 xOverX0 = xOverX0PipeTrks;
1243 lengthTimesMeanDensity = xTimesRhoPipeTrks;
1244 lengthTimesMeanDensity *= dir;
1245 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1246 break;
1247 }
1248
1249 return 1;
1250}
1251//------------------------------------------------------------------------
1252Int_t AliITStrackerHLT::CorrectForShieldMaterial(AliHLTITSTrack *t,
1253 TString shield,
1254 TString direction) {
1255 //-------------------------------------------------------------------
1256 // Propagate beyond SPD or SDD shield and correct for material
1257 // (material budget in different ways according to fUseTGeo value)
1258 // Add time if going outward (PropagateTo or PropagateToTGeo)
1259 //-------------------------------------------------------------------
1260
1261 // Define budget mode:
1262 // 0: material from AliITSRecoParam (hard coded)
1263 // 1: material from TGeo in steps of X cm (on the fly)
1264 // X = AliITSRecoParam::GetStepSizeTGeo()
1265 // 2: material from lut
1266 // 3: material from TGeo in one step (same for all hypotheses)
1267 Int_t mode;
1268 switch(fUseTGeo) {
1269 case 0:
1270 mode=0;
1271 break;
1272 case 1:
1273 mode=1;
1274 break;
1275 case 2:
1276 mode=2;
1277 break;
1278 case 3:
1279 mode=3;
1280 break;
1281 case 4:
1282 mode=3;
1283 break;
1284 default:
1285 mode=0;
1286 break;
1287 }
1288
1289
1290 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
1291 Double_t rToGo;
1292 Int_t shieldindex=0;
1293 if (shield.Contains("SDD")) { // SDDouter
1294 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
1295 shieldindex=1;
1296 } else if (shield.Contains("SPD")) { // SPDouter
1297 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
1298 shieldindex=0;
1299 } else {
1300 Error("CorrectForShieldMaterial"," Wrong shield name\n");
1301 return 0;
1302 }
1303 Double_t xToGo;
1304 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1305
1306 Double_t xOverX0,x0,lengthTimesMeanDensity;
1307 Int_t nsteps=1;
1308
1309 switch(mode) {
1310 case 0:
1311 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
1312 x0 = AliITSRecoParam::GetX0shield(shieldindex);
1313 lengthTimesMeanDensity = xOverX0*x0;
1314 lengthTimesMeanDensity *= dir;
1315 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1316 break;
1317 case 1:
87cc43e1 1318 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/fRecoParam->GetStepSizeTGeo())+1;
2f399afc 1319 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
1320 break;
1321 case 2:
1322 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
1323 xOverX0 = fxOverX0Shield[shieldindex];
1324 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
1325 lengthTimesMeanDensity *= dir;
1326 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1327 break;
1328 case 3:
1329 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
1330 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1331 ((1.-t->GetSnp())*(1.+t->GetSnp())));
1332 double xOverX0ShieldTrks = TMath::Abs(xOverX0)/angle;
1333 double xTimesRhoShieldTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
1334 xOverX0 = xOverX0ShieldTrks;
1335 lengthTimesMeanDensity = xTimesRhoShieldTrks;
1336 lengthTimesMeanDensity *= dir;
1337 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1338 break;
1339 }
1340
1341 return 1;
1342}
1343//------------------------------------------------------------------------
1344Int_t AliITStrackerHLT::CorrectForLayerMaterial(AliHLTITSTrack *t,
1345 Int_t layerindex,
1346 Double_t oldGlobXYZ[3],
1347 TString direction) {
1348 //-------------------------------------------------------------------
1349 // Propagate beyond layer and correct for material
1350 // (material budget in different ways according to fUseTGeo value)
1351 // Add time if going outward (PropagateTo or PropagateToTGeo)
1352 //-------------------------------------------------------------------
1353
1354 // Define budget mode:
1355 // 0: material from AliITSRecoParam (hard coded)
1356 // 1: material from TGeo in stepsof X cm (on the fly)
1357 // X = AliITSRecoParam::GetStepSizeTGeo()
1358 // 2: material from lut
1359 // 3: material from TGeo in one step (same for all hypotheses)
1360 Int_t mode;
1361 switch(fUseTGeo) {
1362 case 0:
1363 mode=0;
1364 break;
1365 case 1:
1366 mode=1;
1367 break;
1368 case 2:
1369 mode=2;
1370 break;
1371 case 3:
1372 mode=3;
1373 break;
1374 case 4:
1375 mode=3;
1376 break;
1377 default:
1378 mode=0;
1379 break;
1380 }
1381
1382
1383 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
1384
1385 Double_t r=fgLayers[layerindex].GetR();
1386 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
1387
1388 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
1389 Double_t xToGo;
1390 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1391
1392
1393 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
1394 Int_t nsteps=1;
1395
1396 // back before material (no correction)
1397 Double_t rOld,xOld;
1398 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
1399 if (!t->GetLocalXat(rOld,xOld)) return 0;
1400 if (!TransportToX( t, xOld)) return 0;
1401
1402 switch(mode) {
1403 case 0:
1404 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
1405 lengthTimesMeanDensity = xOverX0*x0;
1406 lengthTimesMeanDensity *= dir;
1407 // Bring the track beyond the material
1408 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1409 break;
1410 case 1:
87cc43e1 1411 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/fRecoParam->GetStepSizeTGeo())+1;
2f399afc 1412 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
1413 break;
1414 case 2:
1415 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
1416 xOverX0 = fxOverX0Layer[layerindex];
1417 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
1418 lengthTimesMeanDensity *= dir;
1419 // Bring the track beyond the material
1420 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1421 break;
1422 case 3:
87cc43e1 1423 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/fRecoParam->GetStepSizeTGeo())+1;
2f399afc 1424 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
1425 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1426 ((1.-t->GetSnp())*(1.+t->GetSnp())));
1427 double xOverX0LayerTrks = TMath::Abs(xOverX0)/angle;
1428 double xTimesRhoLayerTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
1429 xOverX0 = xOverX0LayerTrks;
1430 lengthTimesMeanDensity = xTimesRhoLayerTrks;
1431 lengthTimesMeanDensity *= dir;
1432 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1433 break;
1434 }
1435
1436
1437 return 1;
1438}
1439
1440
1441//------------------------------------------------------------------------
1442Int_t AliITStrackerHLT::CheckDeadZone(AliHLTITSTrack *track,
1443 Int_t ilayer,Int_t idet,
1444 Double_t dz,Double_t dy,
1445 Bool_t noClusters) const {
1446 //-----------------------------------------------------------------
1447 // This method is used to decide whether to allow a prolongation
1448 // without clusters, because there is a dead zone in the road.
1449 // In this case the return value is > 0:
1450 // return 1: dead zone at z=0,+-7cm in SPD
1451 // return 2: all road is "bad" (dead or noisy) from the OCDB
1452 // return 3: something "bad" (dead or noisy) from the OCDB
1453 //-----------------------------------------------------------------
1454
1455 // check dead zones at z=0,+-7cm in the SPD
87cc43e1 1456 if (ilayer<2 && !fRecoParam->GetAddVirtualClustersInDeadZone()) {
2f399afc 1457 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
1458 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
1459 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
1460 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
1461 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
1462 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
1463 for (Int_t i=0; i<3; i++)
1464 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
1465 AliDebug(2,Form("crack SPD %d",ilayer));
1466 return 1;
1467 }
1468 }
1469
1470 // check bad zones from OCDB
87cc43e1 1471 if (!fRecoParam->GetUseBadZonesFromOCDB()) return 0;
2f399afc 1472
1473 if (idet<0) return 0;
1474
1475 AliHLTITSDetector &det=fgLayers[ilayer].GetDetector(idet);
1476
1477 Int_t detType=-1;
1478 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
1479 if (ilayer==0 || ilayer==1) { // ---------- SPD
1480 detType = 0;
1481 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
1482 detType = 1;
1483 detSizeFactorX *= 2.;
1484 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
1485 detType = 2;
1486 }
1487 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
1488 if (detType==2) segm->SetLayer(ilayer+1);
1489 Float_t detSizeX = detSizeFactorX*segm->Dx();
1490 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
1491
1492 // check if the road overlaps with bad chips
1493 Float_t xloc,zloc;
1494 LocalModuleCoord(ilayer,idet,track,xloc,zloc);
1495 Float_t zlocmin = zloc-dz;
1496 Float_t zlocmax = zloc+dz;
1497 Float_t xlocmin = xloc-dy;
1498 Float_t xlocmax = xloc+dy;
1499 Int_t chipsInRoad[100];
1500
1501 // check if road goes out of detector
1502 Bool_t touchNeighbourDet=kFALSE;
1503 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.5*detSizeX; touchNeighbourDet=kTRUE;}
1504 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.5*detSizeX; touchNeighbourDet=kTRUE;}
1505 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.5*detSizeZ; touchNeighbourDet=kTRUE;}
1506 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.5*detSizeZ; touchNeighbourDet=kTRUE;}
1507 AliDebug(2,Form("layer %d det %d zmim zmax %f %f xmin xmax %f %f %f %f",ilayer,idet,zlocmin,zlocmax,xlocmin,xlocmax,detSizeZ,detSizeX));
1508
1509 // check if this detector is bad
1510 if (det.IsBad()) {
1511 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
1512 if(!touchNeighbourDet) {
1513 return 2; // all detectors in road are bad
1514 } else {
1515 return 3; // at least one is bad
1516 }
1517 }
1518
1519 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
1520 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
1521 if (!nChipsInRoad) return 0;
1522
1523 Bool_t anyBad=kFALSE,anyGood=kFALSE;
1524 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
1525 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
1526 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
1527 if (det.IsChipBad(chipsInRoad[iCh])) {
1528 anyBad=kTRUE;
1529 } else {
1530 anyGood=kTRUE;
1531 }
1532 }
1533
1534 if (!anyGood) {
1535 if(!touchNeighbourDet) {
1536 AliDebug(2,"all bad in road");
1537 return 2; // all chips in road are bad
1538 } else {
1539 return 3; // at least a bad chip in road
1540 }
1541 }
1542
1543 if (anyBad) {
1544 AliDebug(2,"at least a bad in road");
1545 return 3; // at least a bad chip in road
1546 }
1547
1548
87cc43e1 1549 if (!fRecoParam->GetUseSingleBadChannelsFromOCDB()
2f399afc 1550 || ilayer==4 || ilayer==5 // SSD
1551 || !noClusters) return 0;
1552
1553 // There are no clusters in road: check if there is at least
1554 // a bad SPD pixel or SDD anode
1555
1556 Int_t idetInITS=idet;
1557 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
1558
1559 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
1560 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
1561 return 3;
1562 }
87cc43e1 1563 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > fRecoParam->GetMinFractionOfBadInRoad()) return 3;
2f399afc 1564
1565 return 0;
1566}
1567//------------------------------------------------------------------------
1568Bool_t AliITStrackerHLT::LocalModuleCoord(Int_t ilayer,Int_t idet,
1569 const AliHLTITSTrack *track,
1570 Float_t &xloc,Float_t &zloc) const {
1571 //-----------------------------------------------------------------
1572 // Gives position of track in local module ref. frame
1573 //-----------------------------------------------------------------
1574
1575 xloc=0.;
1576 zloc=0.;
1577
1578 if(idet<0) return kFALSE;
1579
1580 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
1581
1582 Int_t lad = Int_t(idet/ndet) + 1;
1583
1584 Int_t det = idet - (lad-1)*ndet + 1;
1585
1586 Double_t xyzGlob[3],xyzLoc[3];
1587
1588 AliHLTITSDetector &detector = fgLayers[ilayer].GetDetector(idet);
1589 // take into account the misalignment: xyz at real detector plane
1590 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
1591
1592 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
1593
1594 xloc = (Float_t)xyzLoc[0];
1595 zloc = (Float_t)xyzLoc[2];
1596
1597 return kTRUE;
1598}
1599
1600//------------------------------------------------------------------------
1601Bool_t AliITStrackerHLT::ComputeRoad(AliHLTITSTrack* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1602 //--------------------------------------------------------------------
1603 // This function computes the rectangular road for this track
1604 //--------------------------------------------------------------------
1605
1606
1607 AliHLTITSDetector &det = fgLayers[ilayer].GetDetector(idet);
1608 // take into account the misalignment: propagate track to misaligned detector plane
1609 if (!TransportToPhiX( track, det.GetPhi(),det.GetRmisal() ) ) return kFALSE;
1610
87cc43e1 1611 Double_t dz=fRecoParam->GetNSigmaRoadZ()*
2f399afc 1612 TMath::Sqrt(track->GetSigmaZ2() +
87cc43e1 1613 fRecoParam->GetNSigmaZLayerForRoadZ()*
1614 fRecoParam->GetNSigmaZLayerForRoadZ()*
1615 fRecoParam->GetSigmaZ2(ilayer));
1616 Double_t dy=fRecoParam->GetNSigmaRoadY()*
2f399afc 1617 TMath::Sqrt(track->GetSigmaY2() +
87cc43e1 1618 fRecoParam->GetNSigmaYLayerForRoadY()*
1619 fRecoParam->GetNSigmaYLayerForRoadY()*
1620 fRecoParam->GetSigmaY2(ilayer));
2f399afc 1621
1622 // track at boundary between detectors, enlarge road
1623 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1624 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1625 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1626 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1627 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1628 Float_t tgl = TMath::Abs(track->GetTgl());
1629 if (tgl > 1.) tgl=1.;
1630 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1631 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1632 Float_t snp = TMath::Abs(track->GetSnp());
87cc43e1 1633 if (snp > fRecoParam->GetMaxSnp()) return kFALSE;
2f399afc 1634 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1635 } // boundary
1636
1637 // add to the road a term (up to 2-3 mm) to deal with misalignments
87cc43e1 1638 dy = TMath::Sqrt(dy*dy + fRecoParam->GetRoadMisal()*fRecoParam->GetRoadMisal());
1639 dz = TMath::Sqrt(dz*dz + fRecoParam->GetRoadMisal()*fRecoParam->GetRoadMisal());
2f399afc 1640
1641 Double_t r = fgLayers[ilayer].GetR();
1642 zmin = track->GetZ() - dz;
1643 zmax = track->GetZ() + dz;
1644 ymin = track->GetY() + r*det.GetPhi() - dy;
1645 ymax = track->GetY() + r*det.GetPhi() + dy;
1646
1647 // bring track back to idead detector plane
1648 if (!TransportToPhiX( track, det.GetPhi(),det.GetR())) return kFALSE;
1649
1650 return kTRUE;
1651}
28df84b5 1652