]>
Commit | Line | Data |
---|---|---|
4811a3f4 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id: AliTrackerBase.cxx 38069 2009-12-24 16:56:18Z belikov $ */ | |
17 | ||
18 | //------------------------------------------------------------------------- | |
19 | // Implementation of the AliTrackerBase class | |
075f4221 | 20 | // that is the base for the AliTracker class |
21 | // Origin: Marian.Ivanov@cern.ch | |
4811a3f4 | 22 | //------------------------------------------------------------------------- |
23 | #include <TClass.h> | |
24 | #include <TMath.h> | |
25 | #include <TGeoManager.h> | |
26 | ||
27 | #include "AliLog.h" | |
28 | #include "AliTrackerBase.h" | |
29 | #include "AliExternalTrackParam.h" | |
075f4221 | 30 | #include "AliTrackPointArray.h" |
363db7c3 | 31 | #include "TVectorD.h" |
4811a3f4 | 32 | |
33 | extern TGeoManager *gGeoManager; | |
34 | ||
35 | ClassImp(AliTrackerBase) | |
36 | ||
37 | AliTrackerBase::AliTrackerBase(): | |
38 | TObject(), | |
39 | fX(0), | |
40 | fY(0), | |
41 | fZ(0), | |
42 | fSigmaX(0.005), | |
43 | fSigmaY(0.005), | |
44 | fSigmaZ(0.010) | |
45 | { | |
46 | //-------------------------------------------------------------------- | |
47 | // The default constructor. | |
48 | //-------------------------------------------------------------------- | |
49 | if (!TGeoGlobalMagField::Instance()->GetField()) | |
50 | AliWarning("Field map is not set."); | |
51 | } | |
52 | ||
53 | //__________________________________________________________________________ | |
54 | AliTrackerBase::AliTrackerBase(const AliTrackerBase &atr): | |
55 | TObject(atr), | |
56 | fX(atr.fX), | |
57 | fY(atr.fY), | |
58 | fZ(atr.fZ), | |
59 | fSigmaX(atr.fSigmaX), | |
60 | fSigmaY(atr.fSigmaY), | |
61 | fSigmaZ(atr.fSigmaZ) | |
62 | { | |
63 | //-------------------------------------------------------------------- | |
64 | // The default constructor. | |
65 | //-------------------------------------------------------------------- | |
66 | if (!TGeoGlobalMagField::Instance()->GetField()) | |
67 | AliWarning("Field map is not set."); | |
68 | } | |
69 | ||
70 | //__________________________________________________________________________ | |
71 | Double_t AliTrackerBase::GetBz() | |
72 | { | |
73 | AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); | |
93ba142d | 74 | if (!fld) { |
75 | AliFatalClass("Field is not loaded"); | |
76 | //if (!fld) | |
77 | return 0.5*kAlmost0Field; | |
78 | } | |
4811a3f4 | 79 | Double_t bz = fld->SolenoidField(); |
80 | return TMath::Sign(0.5*kAlmost0Field,bz) + bz; | |
81 | } | |
82 | ||
83 | //__________________________________________________________________________ | |
84 | Double_t AliTrackerBase::GetBz(const Double_t *r) { | |
85 | //------------------------------------------------------------------ | |
86 | // Returns Bz (kG) at the point "r" . | |
87 | //------------------------------------------------------------------ | |
88 | AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); | |
93ba142d | 89 | if (!fld) { |
90 | AliFatalClass("Field is not loaded"); | |
91 | // if (!fld) | |
92 | return 0.5*kAlmost0Field; | |
93 | } | |
4811a3f4 | 94 | Double_t bz = fld->GetBz(r); |
95 | return TMath::Sign(0.5*kAlmost0Field,bz) + bz; | |
96 | } | |
97 | ||
98 | //__________________________________________________________________________ | |
99 | void AliTrackerBase::GetBxByBz(const Double_t r[3], Double_t b[3]) { | |
100 | //------------------------------------------------------------------ | |
101 | // Returns Bx, By and Bz (kG) at the point "r" . | |
102 | //------------------------------------------------------------------ | |
103 | AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); | |
104 | if (!fld) { | |
93ba142d | 105 | AliFatalClass("Field is not loaded"); |
106 | // b[0] = b[1] = 0.; | |
107 | // b[2] = 0.5*kAlmost0Field; | |
108 | return; | |
4811a3f4 | 109 | } |
110 | ||
111 | if (fld->IsUniform()) { | |
112 | b[0] = b[1] = 0.; | |
113 | b[2] = fld->SolenoidField(); | |
114 | } else { | |
115 | fld->Field(r,b); | |
116 | } | |
117 | b[2] = (TMath::Sign(0.5*kAlmost0Field,b[2]) + b[2]); | |
118 | return; | |
119 | } | |
120 | ||
121 | Double_t AliTrackerBase::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam) | |
122 | { | |
123 | // | |
124 | // Calculate mean material budget and material properties between | |
125 | // the points "start" and "end". | |
126 | // | |
127 | // "mparam" - parameters used for the energy and multiple scattering | |
128 | // corrections: | |
129 | // | |
130 | // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3] | |
131 | // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional] | |
132 | // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional] | |
133 | // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional] | |
134 | // mparam[4] - length: sum(x_i) [cm] | |
135 | // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional] | |
136 | // mparam[6] - number of boundary crosses | |
137 | // | |
138 | // Origin: Marian Ivanov, Marian.Ivanov@cern.ch | |
139 | // | |
140 | // Corrections and improvements by | |
141 | // Andrea Dainese, Andrea.Dainese@lnl.infn.it, | |
142 | // Andrei Gheata, Andrei.Gheata@cern.ch | |
143 | // | |
144 | ||
145 | mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0; | |
146 | mparam[4]=0; mparam[5]=0; mparam[6]=0; | |
147 | // | |
148 | Double_t bparam[6]; // total parameters | |
149 | Double_t lparam[6]; // local parameters | |
150 | ||
151 | for (Int_t i=0;i<6;i++) bparam[i]=0; | |
152 | ||
153 | if (!gGeoManager) { | |
2ad9a021 | 154 | AliFatalClass("No TGeo\n"); |
4811a3f4 | 155 | return 0.; |
156 | } | |
157 | // | |
158 | Double_t length; | |
159 | Double_t dir[3]; | |
160 | length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+ | |
161 | (end[1]-start[1])*(end[1]-start[1])+ | |
162 | (end[2]-start[2])*(end[2]-start[2])); | |
163 | mparam[4]=length; | |
164 | if (length<TGeoShape::Tolerance()) return 0.0; | |
165 | Double_t invlen = 1./length; | |
166 | dir[0] = (end[0]-start[0])*invlen; | |
167 | dir[1] = (end[1]-start[1])*invlen; | |
168 | dir[2] = (end[2]-start[2])*invlen; | |
169 | ||
170 | // Initialize start point and direction | |
171 | TGeoNode *currentnode = 0; | |
172 | TGeoNode *startnode = gGeoManager->InitTrack(start, dir); | |
173 | if (!startnode) { | |
174 | AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f", | |
175 | start[0],start[1],start[2])); | |
176 | return 0.0; | |
177 | } | |
178 | TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial(); | |
179 | lparam[0] = material->GetDensity(); | |
180 | lparam[1] = material->GetRadLen(); | |
181 | lparam[2] = material->GetA(); | |
182 | lparam[3] = material->GetZ(); | |
183 | lparam[4] = length; | |
184 | lparam[5] = lparam[3]/lparam[2]; | |
185 | if (material->IsMixture()) { | |
186 | TGeoMixture * mixture = (TGeoMixture*)material; | |
187 | lparam[5] =0; | |
188 | Double_t sum =0; | |
189 | for (Int_t iel=0;iel<mixture->GetNelements();iel++){ | |
190 | sum += mixture->GetWmixt()[iel]; | |
191 | lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel]; | |
192 | } | |
193 | lparam[5]/=sum; | |
194 | } | |
195 | ||
196 | // Locate next boundary within length without computing safety. | |
197 | // Propagate either with length (if no boundary found) or just cross boundary | |
198 | gGeoManager->FindNextBoundaryAndStep(length, kFALSE); | |
199 | Double_t step = 0.0; // Step made | |
200 | Double_t snext = gGeoManager->GetStep(); | |
201 | // If no boundary within proposed length, return current density | |
202 | if (!gGeoManager->IsOnBoundary()) { | |
203 | mparam[0] = lparam[0]; | |
204 | mparam[1] = lparam[4]/lparam[1]; | |
205 | mparam[2] = lparam[2]; | |
206 | mparam[3] = lparam[3]; | |
207 | mparam[4] = lparam[4]; | |
208 | return lparam[0]; | |
209 | } | |
210 | // Try to cross the boundary and see what is next | |
211 | Int_t nzero = 0; | |
212 | while (length>TGeoShape::Tolerance()) { | |
213 | currentnode = gGeoManager->GetCurrentNode(); | |
214 | if (snext<2.*TGeoShape::Tolerance()) nzero++; | |
215 | else nzero = 0; | |
216 | if (nzero>3) { | |
217 | // This means navigation has problems on one boundary | |
218 | // Try to cross by making a small step | |
219 | AliErrorClass("Cannot cross boundary\n"); | |
220 | mparam[0] = bparam[0]/step; | |
221 | mparam[1] = bparam[1]; | |
222 | mparam[2] = bparam[2]/step; | |
223 | mparam[3] = bparam[3]/step; | |
224 | mparam[5] = bparam[5]/step; | |
225 | mparam[4] = step; | |
226 | mparam[0] = 0.; // if crash of navigation take mean density 0 | |
227 | mparam[1] = 1000000; // and infinite rad length | |
228 | return bparam[0]/step; | |
229 | } | |
230 | mparam[6]+=1.; | |
231 | step += snext; | |
232 | bparam[1] += snext/lparam[1]; | |
233 | bparam[2] += snext*lparam[2]; | |
234 | bparam[3] += snext*lparam[3]; | |
235 | bparam[5] += snext*lparam[5]; | |
236 | bparam[0] += snext*lparam[0]; | |
237 | ||
238 | if (snext>=length) break; | |
239 | if (!currentnode) break; | |
240 | length -= snext; | |
241 | material = currentnode->GetVolume()->GetMedium()->GetMaterial(); | |
242 | lparam[0] = material->GetDensity(); | |
243 | lparam[1] = material->GetRadLen(); | |
244 | lparam[2] = material->GetA(); | |
245 | lparam[3] = material->GetZ(); | |
246 | lparam[5] = lparam[3]/lparam[2]; | |
247 | if (material->IsMixture()) { | |
248 | TGeoMixture * mixture = (TGeoMixture*)material; | |
249 | lparam[5]=0; | |
250 | Double_t sum =0; | |
251 | for (Int_t iel=0;iel<mixture->GetNelements();iel++){ | |
252 | sum+= mixture->GetWmixt()[iel]; | |
253 | lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel]; | |
254 | } | |
255 | lparam[5]/=sum; | |
256 | } | |
257 | gGeoManager->FindNextBoundaryAndStep(length, kFALSE); | |
258 | snext = gGeoManager->GetStep(); | |
259 | } | |
260 | mparam[0] = bparam[0]/step; | |
261 | mparam[1] = bparam[1]; | |
262 | mparam[2] = bparam[2]/step; | |
263 | mparam[3] = bparam[3]/step; | |
264 | mparam[5] = bparam[5]/step; | |
265 | return bparam[0]/step; | |
266 | } | |
267 | ||
268 | ||
269 | Bool_t | |
270 | AliTrackerBase::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo, | |
cd5c09f1 | 271 | Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Int_t sign, Bool_t addTimeStep, Bool_t correctMaterialBudget){ |
4811a3f4 | 272 | //---------------------------------------------------------------- |
273 | // | |
274 | // Propagates the track to the plane X=xk (cm) using the magnetic field map | |
275 | // and correcting for the crossed material. | |
276 | // | |
1d26da6d | 277 | // mass - mass used in propagation - used for energy loss correction (if <0 then q = 2) |
4811a3f4 | 278 | // maxStep - maximal step for propagation |
279 | // | |
280 | // Origin: Marian Ivanov, Marian.Ivanov@cern.ch | |
281 | // | |
282 | //---------------------------------------------------------------- | |
283 | const Double_t kEpsilon = 0.00001; | |
284 | Double_t xpos = track->GetX(); | |
13545423 | 285 | Int_t dir = (xpos<xToGo) ? 1:-1; |
4811a3f4 | 286 | // |
287 | while ( (xToGo-xpos)*dir > kEpsilon){ | |
288 | Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); | |
289 | Double_t x = xpos+step; | |
290 | Double_t xyz0[3],xyz1[3],param[7]; | |
291 | track->GetXYZ(xyz0); //starting global position | |
292 | ||
293 | Double_t bz=GetBz(xyz0); // getting the local Bz | |
4811a3f4 | 294 | if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE; // no prolongation |
295 | xyz1[2]+=kEpsilon; // waiting for bug correction in geo | |
cd5c09f1 | 296 | |
7391b629 | 297 | if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE; |
4811a3f4 | 298 | if (!track->PropagateTo(x,bz)) return kFALSE; |
299 | ||
cd5c09f1 | 300 | if (correctMaterialBudget){ |
301 | MeanMaterialBudget(xyz0,xyz1,param); | |
302 | Double_t xrho=param[0]*param[4], xx0=param[1]; | |
303 | if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed | |
304 | else { // determine automatically the sign from direction | |
305 | if (dir>0) xrho = -xrho; // outward should be negative | |
306 | } | |
307 | // | |
308 | if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE; | |
13545423 | 309 | } |
cd5c09f1 | 310 | |
4811a3f4 | 311 | if (rotateTo){ |
dc981549 | 312 | track->GetXYZ(xyz1); // global position |
313 | Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]); | |
7391b629 | 314 | if (maxSnp>0) { |
315 | if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE; | |
cd5c09f1 | 316 | |
7391b629 | 317 | // |
318 | Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha()); | |
319 | Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf)); | |
320 | Double_t sinNew = sf*ca - cf*sa; | |
cd5c09f1 | 321 | if (TMath::Abs(sinNew) >= maxSnp) return kFALSE; |
322 | ||
7391b629 | 323 | } |
dc981549 | 324 | if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE; |
cd5c09f1 | 325 | |
4811a3f4 | 326 | } |
327 | xpos = track->GetX(); | |
13545423 | 328 | if (addTimeStep && track->IsStartedTimeIntegral()) { |
dc981549 | 329 | if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted |
13545423 | 330 | Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2]; |
331 | Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ); | |
332 | if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction | |
333 | else { // determine automatically the sign from direction | |
334 | if (dir<0) d = -d; | |
335 | } | |
336 | track->AddTimeStep(d); | |
337 | } | |
4811a3f4 | 338 | } |
339 | return kTRUE; | |
340 | } | |
341 | ||
69749404 | 342 | Int_t AliTrackerBase::PropagateTrackTo2(AliExternalTrackParam *track, Double_t xToGo, |
343 | Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Int_t sign, Bool_t addTimeStep, Bool_t correctMaterialBudget){ | |
344 | //---------------------------------------------------------------- | |
345 | // | |
346 | // Propagates the track to the plane X=xk (cm) using the magnetic field map | |
347 | // and correcting for the crossed material. | |
348 | // | |
349 | // mass - mass used in propagation - used for energy loss correction | |
350 | // maxStep - maximal step for propagation | |
351 | // | |
352 | // Origin: Marian Ivanov, Marian.Ivanov@cern.ch | |
353 | // | |
354 | //---------------------------------------------------------------- | |
355 | const Double_t kEpsilon = 0.00001; | |
356 | Double_t xpos = track->GetX(); | |
357 | Int_t dir = (xpos<xToGo) ? 1:-1; | |
358 | // | |
359 | while ( (xToGo-xpos)*dir > kEpsilon){ | |
360 | Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); | |
361 | Double_t x = xpos+step; | |
362 | Double_t xyz0[3],xyz1[3],param[7]; | |
363 | track->GetXYZ(xyz0); //starting global position | |
364 | ||
365 | Double_t bz=GetBz(xyz0); // getting the local Bz | |
366 | if (!track->GetXYZAt(x,bz,xyz1)) return -1; // no prolongation | |
367 | xyz1[2]+=kEpsilon; // waiting for bug correction in geo | |
368 | ||
369 | if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return -2; | |
370 | if (!track->PropagateTo(x,bz)) return -3; | |
371 | ||
372 | if (correctMaterialBudget){ | |
373 | MeanMaterialBudget(xyz0,xyz1,param); | |
374 | Double_t xrho=param[0]*param[4], xx0=param[1]; | |
375 | if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed | |
376 | else { // determine automatically the sign from direction | |
377 | if (dir>0) xrho = -xrho; // outward should be negative | |
378 | } | |
379 | // | |
380 | if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return -4; | |
381 | } | |
382 | ||
383 | if (rotateTo){ | |
384 | track->GetXYZ(xyz1); // global position | |
385 | Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]); | |
386 | if (maxSnp>0) { | |
387 | if (TMath::Abs(track->GetSnp()) >= maxSnp) return -5; | |
388 | ||
389 | // | |
390 | Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha()); | |
391 | Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf)); | |
392 | Double_t sinNew = sf*ca - cf*sa; | |
393 | if (TMath::Abs(sinNew) >= maxSnp) return -6; | |
394 | ||
395 | } | |
396 | if (!track->AliExternalTrackParam::Rotate(alphan)) return -7; | |
397 | ||
398 | } | |
399 | xpos = track->GetX(); | |
400 | if (addTimeStep && track->IsStartedTimeIntegral()) { | |
401 | if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted | |
402 | Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2]; | |
403 | Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ); | |
404 | if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction | |
405 | else { // determine automatically the sign from direction | |
406 | if (dir<0) d = -d; | |
407 | } | |
408 | track->AddTimeStep(d); | |
409 | } | |
410 | } | |
411 | return 1; | |
412 | } | |
413 | ||
4811a3f4 | 414 | Bool_t |
415 | AliTrackerBase::PropagateTrackToBxByBz(AliExternalTrackParam *track, | |
13545423 | 416 | Double_t xToGo,Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Int_t sign, Bool_t addTimeStep){ |
4811a3f4 | 417 | //---------------------------------------------------------------- |
418 | // | |
419 | // Propagates the track to the plane X=xk (cm) | |
420 | // taking into account all the three components of the magnetic field | |
421 | // and correcting for the crossed material. | |
422 | // | |
1d26da6d | 423 | // mass - mass used in propagation - used for energy loss correction (if <0 then q=2) |
4811a3f4 | 424 | // maxStep - maximal step for propagation |
425 | // | |
426 | // Origin: Marian Ivanov, Marian.Ivanov@cern.ch | |
427 | // | |
428 | //---------------------------------------------------------------- | |
429 | const Double_t kEpsilon = 0.00001; | |
430 | Double_t xpos = track->GetX(); | |
13545423 | 431 | Int_t dir = (xpos<xToGo) ? 1:-1; |
4811a3f4 | 432 | // |
433 | while ( (xToGo-xpos)*dir > kEpsilon){ | |
434 | Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); | |
435 | Double_t x = xpos+step; | |
436 | Double_t xyz0[3],xyz1[3],param[7]; | |
437 | track->GetXYZ(xyz0); //starting global position | |
438 | ||
439 | Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz | |
440 | ||
441 | if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE; // no prolongation | |
442 | xyz1[2]+=kEpsilon; // waiting for bug correction in geo | |
443 | ||
7391b629 | 444 | if (maxSnp>0 && TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE; |
4811a3f4 | 445 | if (!track->PropagateToBxByBz(x,b)) return kFALSE; |
446 | ||
13545423 | 447 | MeanMaterialBudget(xyz0,xyz1,param); |
448 | Double_t xrho=param[0]*param[4], xx0=param[1]; | |
449 | if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed | |
450 | else { // determine automatically the sign from direction | |
451 | if (dir>0) xrho = -xrho; // outward should be negative | |
452 | } | |
453 | // | |
4811a3f4 | 454 | if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE; |
455 | if (rotateTo){ | |
dc981549 | 456 | track->GetXYZ(xyz1); // global position |
457 | Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]); | |
7391b629 | 458 | if (maxSnp>0) { |
459 | if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE; | |
460 | Double_t ca=TMath::Cos(alphan-track->GetAlpha()), sa=TMath::Sin(alphan-track->GetAlpha()); | |
461 | Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf)); | |
462 | Double_t sinNew = sf*ca - cf*sa; | |
463 | if (TMath::Abs(sinNew) >= maxSnp) return kFALSE; | |
464 | } | |
dc981549 | 465 | if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE; |
4811a3f4 | 466 | } |
13545423 | 467 | xpos = track->GetX(); |
468 | if (addTimeStep && track->IsStartedTimeIntegral()) { | |
dc981549 | 469 | if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted |
13545423 | 470 | Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2]; |
471 | Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ); | |
472 | if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction | |
473 | else { // determine automatically the sign from direction | |
474 | if (dir<0) d = -d; | |
475 | } | |
476 | track->AddTimeStep(d); | |
477 | } | |
4811a3f4 | 478 | } |
479 | return kTRUE; | |
480 | } | |
481 | ||
482 | Double_t AliTrackerBase::GetTrackPredictedChi2(AliExternalTrackParam *track, | |
483 | Double_t mass, Double_t step, | |
484 | const AliExternalTrackParam *backup) { | |
485 | // | |
486 | // This function brings the "track" with particle "mass" [GeV] | |
487 | // to the same local coord. system and the same reference plane as | |
488 | // of the "backup", doing it in "steps" [cm]. | |
489 | // Then, it calculates the 5D predicted Chi2 for these two tracks | |
490 | // | |
491 | Double_t chi2=kVeryBig; | |
492 | Double_t alpha=backup->GetAlpha(); | |
493 | if (!track->Rotate(alpha)) return chi2; | |
494 | ||
495 | Double_t xb=backup->GetX(); | |
496 | Double_t sign=(xb < track->GetX()) ? 1. : -1.; | |
497 | if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2; | |
498 | ||
499 | chi2=track->GetPredictedChi2(backup); | |
500 | ||
501 | return chi2; | |
502 | } | |
075f4221 | 503 | |
504 | ||
505 | ||
506 | ||
507 | Double_t AliTrackerBase::MakeC(Double_t x1,Double_t y1, | |
508 | Double_t x2,Double_t y2, | |
509 | Double_t x3,Double_t y3) | |
510 | { | |
511 | //----------------------------------------------------------------- | |
512 | // Initial approximation of the track curvature | |
513 | //----------------------------------------------------------------- | |
514 | x3 -=x1; | |
515 | x2 -=x1; | |
516 | y3 -=y1; | |
517 | y2 -=y1; | |
518 | // | |
519 | Double_t det = x3*y2-x2*y3; | |
520 | if (TMath::Abs(det)<1e-10) { | |
521 | return 0; | |
522 | } | |
523 | // | |
524 | Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det; | |
525 | Double_t x0 = x3*0.5-y3*u; | |
526 | Double_t y0 = y3*0.5+x3*u; | |
527 | Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0); | |
528 | if (det>0) c2*=-1; | |
529 | return c2; | |
530 | } | |
531 | ||
532 | Double_t AliTrackerBase::MakeSnp(Double_t x1,Double_t y1, | |
533 | Double_t x2,Double_t y2, | |
534 | Double_t x3,Double_t y3) | |
535 | { | |
536 | //----------------------------------------------------------------- | |
537 | // Initial approximation of the track snp | |
538 | //----------------------------------------------------------------- | |
539 | x3 -=x1; | |
540 | x2 -=x1; | |
541 | y3 -=y1; | |
542 | y2 -=y1; | |
543 | // | |
544 | Double_t det = x3*y2-x2*y3; | |
545 | if (TMath::Abs(det)<1e-10) { | |
546 | return 0; | |
547 | } | |
548 | // | |
549 | Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det; | |
550 | Double_t x0 = x3*0.5-y3*u; | |
551 | Double_t y0 = y3*0.5+x3*u; | |
552 | Double_t c2 = 1./TMath::Sqrt(x0*x0+y0*y0); | |
553 | x0*=c2; | |
554 | x0=TMath::Abs(x0); | |
555 | if (y2*x2<0.) x0*=-1; | |
556 | return x0; | |
557 | } | |
558 | ||
559 | Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1, | |
560 | Double_t x2,Double_t y2, | |
561 | Double_t z1,Double_t z2, Double_t c) | |
562 | { | |
563 | //----------------------------------------------------------------- | |
564 | // Initial approximation of the tangent of the track dip angle | |
565 | //----------------------------------------------------------------- | |
566 | // | |
e94b22cc | 567 | const Double_t kEpsilon =0.00001; |
075f4221 | 568 | x2-=x1; |
569 | y2-=y1; | |
570 | z2-=z1; | |
571 | Double_t d = TMath::Sqrt(x2*x2+y2*y2); // distance straight line | |
572 | if (TMath::Abs(d*c*0.5)>1) return 0; | |
573 | Double_t angle2 = TMath::ASin(d*c*0.5); | |
e94b22cc | 574 | if (TMath::Abs(angle2)>kEpsilon) { |
575 | angle2 = z2*TMath::Abs(c/(angle2*2.)); | |
576 | }else{ | |
577 | angle2=z2/d; | |
578 | } | |
075f4221 | 579 | return angle2; |
580 | } | |
581 | ||
582 | ||
583 | Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1, | |
584 | Double_t x2,Double_t y2, | |
585 | Double_t z1,Double_t z2) | |
586 | { | |
587 | //----------------------------------------------------------------- | |
588 | // Initial approximation of the tangent of the track dip angle | |
589 | //----------------------------------------------------------------- | |
590 | return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); | |
591 | } | |
592 | ||
593 | ||
594 | AliExternalTrackParam * AliTrackerBase::MakeSeed( AliTrackPoint &point0, AliTrackPoint &point1, AliTrackPoint &point2){ | |
595 | // | |
596 | // Make Seed - AliExternalTrackParam from input 3 points | |
597 | // returning seed in local frame of point0 | |
598 | // | |
599 | Double_t xyz0[3]={0,0,0}; | |
600 | Double_t xyz1[3]={0,0,0}; | |
601 | Double_t xyz2[3]={0,0,0}; | |
602 | Double_t alpha=point0.GetAngle(); | |
603 | Double_t xyz[3]={point0.GetX(),point0.GetY(),point0.GetZ()}; | |
604 | Double_t bxyz[3]; GetBxByBz(xyz,bxyz); | |
605 | Double_t bz = bxyz[2]; | |
606 | // | |
607 | // get points in frame of point 0 | |
608 | // | |
609 | AliTrackPoint p0r = point0.Rotate(alpha); | |
610 | AliTrackPoint p1r = point1.Rotate(alpha); | |
611 | AliTrackPoint p2r = point2.Rotate(alpha); | |
612 | xyz0[0]=p0r.GetX(); | |
613 | xyz0[1]=p0r.GetY(); | |
614 | xyz0[2]=p0r.GetZ(); | |
615 | xyz1[0]=p1r.GetX(); | |
616 | xyz1[1]=p1r.GetY(); | |
617 | xyz1[2]=p1r.GetZ(); | |
618 | xyz2[0]=p2r.GetX(); | |
619 | xyz2[1]=p2r.GetY(); | |
620 | xyz2[2]=p2r.GetZ(); | |
621 | // | |
622 | // make covariance estimate | |
623 | // | |
624 | Double_t covar[15]; | |
625 | Double_t param[5]={0,0,0,0,0}; | |
626 | for (Int_t m=0; m<15; m++) covar[m]=0; | |
627 | // | |
628 | // calculate intitial param | |
629 | param[0]=xyz0[1]; | |
630 | param[1]=xyz0[2]; | |
631 | param[2]=MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]); | |
632 | param[4]=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]); | |
633 | param[3]=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2],param[4]); | |
634 | ||
635 | //covariance matrix - only diagonal elements | |
636 | //Double_t dist=p0r.GetX()-p2r.GetX(); | |
637 | Double_t deltaP=0; | |
638 | covar[0]= p0r.GetCov()[3]; | |
639 | covar[2]= p0r.GetCov()[5]; | |
640 | //sigma snp | |
641 | deltaP= (MakeSnp(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[2]); | |
642 | covar[5]+= deltaP*deltaP; | |
643 | deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[2]); | |
644 | covar[5]+= deltaP*deltaP; | |
645 | deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p1r.GetCov()[3]))-param[2]); | |
646 | covar[5]+= deltaP*deltaP; | |
647 | //sigma tgl | |
648 | // | |
649 | deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2]+TMath::Sqrt(p1r.GetCov()[5]),xyz1[2],param[4])-param[3]; | |
650 | covar[9]+= deltaP*deltaP; | |
651 | deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2]+TMath::Sqrt(p1r.GetCov()[5]),param[4])-param[3]; | |
652 | covar[9]+= deltaP*deltaP; | |
653 | // | |
654 | ||
655 | deltaP=MakeC(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[4]; | |
656 | covar[14]+= deltaP*deltaP; | |
657 | deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[4]; | |
658 | covar[14]+= deltaP*deltaP; | |
659 | deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p2r.GetCov()[3]))-param[4]; | |
660 | covar[14]+= deltaP*deltaP; | |
661 | ||
662 | covar[14]/=(bz*kB2C)*(bz*kB2C); | |
663 | param[4]/=(bz*kB2C); // transform to 1/pt | |
664 | AliExternalTrackParam * trackParam = new AliExternalTrackParam(xyz0[0],alpha,param, covar); | |
665 | if (0) { | |
666 | // consistency check -to put warnings here | |
667 | // small disagrement once Track extrapolation used | |
668 | // nice agreement in seeds with MC track parameters - problem in extrapoloation - to be fixed | |
669 | // to check later | |
670 | Double_t y1,y2,z1,z2; | |
671 | trackParam->GetYAt(xyz1[0],bz,y1); | |
672 | trackParam->GetZAt(xyz1[0],bz,z1); | |
673 | trackParam->GetYAt(xyz2[0],bz,y2); | |
674 | trackParam->GetZAt(xyz2[0],bz,z2); | |
675 | if (TMath::Abs(y1-xyz1[1])> TMath::Sqrt(p1r.GetCov()[3]*5)){ | |
676 | AliWarningClass("Seeding problem y1\n"); | |
677 | } | |
678 | if (TMath::Abs(y2-xyz2[1])> TMath::Sqrt(p2r.GetCov()[3]*5)){ | |
679 | AliWarningClass("Seeding problem y2\n"); | |
680 | } | |
681 | if (TMath::Abs(z1-xyz1[2])> TMath::Sqrt(p1r.GetCov()[5]*5)){ | |
682 | AliWarningClass("Seeding problem z1\n"); | |
683 | } | |
684 | } | |
685 | return trackParam; | |
686 | } | |
687 | ||
688 | Double_t AliTrackerBase::FitTrack(AliExternalTrackParam * trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep){ | |
689 | // | |
690 | // refit the track - trackParam using the points in point array | |
691 | // | |
692 | const Double_t kMaxSnp=0.99; | |
693 | if (!trackParam) return 0; | |
694 | Int_t npoints=pointArray->GetNPoints(); | |
695 | AliTrackPoint point,point2; | |
696 | Double_t pointPos[2]={0,0}; | |
697 | Double_t pointCov[3]={0,0,0}; | |
698 | // choose coordinate frame | |
699 | // in standard way the coordinate frame should be changed point by point | |
700 | // Some problems with rotation observed | |
701 | // rotate method of AliExternalTrackParam should be revisited | |
702 | pointArray->GetPoint(point,0); | |
703 | pointArray->GetPoint(point2,npoints-1); | |
704 | Double_t alpha=TMath::ATan2(point.GetY()-point2.GetY(), point.GetX()-point2.GetX()); | |
705 | ||
706 | for (Int_t ipoint=npoints-1; ipoint>0; ipoint-=1){ | |
707 | pointArray->GetPoint(point,ipoint); | |
708 | AliTrackPoint pr = point.Rotate(alpha); | |
709 | trackParam->Rotate(alpha); | |
710 | Bool_t status = PropagateTrackTo(trackParam,pr.GetX(),mass,maxStep,kFALSE,kMaxSnp); | |
711 | if(!status){ | |
712 | AliWarningClass("Problem to propagate\n"); | |
713 | break; | |
714 | } | |
715 | if (TMath::Abs(trackParam->GetSnp())>kMaxSnp){ | |
716 | AliWarningClass("sin(phi) > kMaxSnp \n"); | |
717 | break; | |
718 | } | |
719 | pointPos[0]=pr.GetY();//local y | |
720 | pointPos[1]=pr.GetZ();//local z | |
721 | pointCov[0]=pr.GetCov()[3];//simay^2 | |
722 | pointCov[1]=pr.GetCov()[4];//sigmayz | |
723 | pointCov[2]=pr.GetCov()[5];//sigmaz^2 | |
724 | trackParam->Update(pointPos,pointCov); | |
725 | } | |
726 | return 0; | |
727 | } | |
363db7c3 | 728 | |
729 | ||
730 | ||
731 | void AliTrackerBase::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){ | |
732 | // | |
733 | // Update track 1 with track 2 | |
734 | // | |
735 | // | |
736 | // | |
737 | TMatrixD vecXk(5,1); // X vector | |
738 | TMatrixD covXk(5,5); // X covariance | |
739 | TMatrixD matHk(5,5); // vector to mesurement | |
740 | TMatrixD measR(5,5); // measurement error | |
741 | TMatrixD vecZk(5,1); // measurement | |
742 | // | |
743 | TMatrixD vecYk(5,1); // Innovation or measurement residual | |
744 | TMatrixD matHkT(5,5); | |
745 | TMatrixD matSk(5,5); // Innovation (or residual) covariance | |
746 | TMatrixD matKk(5,5); // Optimal Kalman gain | |
747 | TMatrixD mat1(5,5); // update covariance matrix | |
748 | TMatrixD covXk2(5,5); // | |
749 | TMatrixD covOut(5,5); | |
750 | // | |
751 | Double_t *param1=(Double_t*) track1.GetParameter(); | |
752 | Double_t *covar1=(Double_t*) track1.GetCovariance(); | |
753 | Double_t *param2=(Double_t*) track2.GetParameter(); | |
754 | Double_t *covar2=(Double_t*) track2.GetCovariance(); | |
755 | // | |
756 | // copy data to the matrix | |
757 | for (Int_t ipar=0; ipar<5; ipar++){ | |
758 | for (Int_t jpar=0; jpar<5; jpar++){ | |
759 | covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)]; | |
760 | measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)]; | |
761 | matHk(ipar,jpar)=0; | |
762 | mat1(ipar,jpar)=0; | |
763 | } | |
764 | vecXk(ipar,0) = param1[ipar]; | |
765 | vecZk(ipar,0) = param2[ipar]; | |
766 | matHk(ipar,ipar)=1; | |
767 | mat1(ipar,ipar)=1; | |
768 | } | |
769 | // | |
770 | // | |
771 | // | |
772 | // | |
773 | // | |
774 | vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual | |
775 | matHkT=matHk.T(); matHk.T(); | |
776 | matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance | |
777 | matSk.Invert(); | |
778 | matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain | |
779 | vecXk += matKk*vecYk; // updated vector | |
780 | covXk2 = (mat1-(matKk*matHk)); | |
781 | covOut = covXk2*covXk; | |
782 | // | |
783 | // | |
784 | // | |
785 | // copy from matrix to parameters | |
786 | if (0) { | |
787 | vecXk.Print(); | |
788 | vecZk.Print(); | |
789 | // | |
790 | measR.Print(); | |
791 | covXk.Print(); | |
792 | covOut.Print(); | |
793 | // | |
794 | track1.Print(); | |
795 | track2.Print(); | |
796 | } | |
797 | ||
798 | for (Int_t ipar=0; ipar<5; ipar++){ | |
799 | param1[ipar]= vecXk(ipar,0) ; | |
800 | for (Int_t jpar=0; jpar<5; jpar++){ | |
801 | covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar); | |
802 | } | |
803 | } | |
804 | } | |
805 | ||
806 |