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