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