]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/AliITSUSeed.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added wrapper volumes to CDR setup
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUSeed.cxx
CommitLineData
f8832015 1#include <TString.h>
32d38de2 2#include <TMath.h>
3#include "AliITSUSeed.h"
716ccba7 4#include "AliLog.h"
3e4e3c23 5#include "AliESDtrack.h"
32d38de2 6using namespace TMath;
7
8ClassImp(AliITSUSeed)
9
10//_________________________________________________________________________
11AliITSUSeed::AliITSUSeed()
f8832015 12: fHitsPattern(0)
b8b59e05 13 ,fNChildren(0)
f8832015 14 ,fClID(0)
15 ,fChi2Glo(0)
16 ,fChi2Cl(0)
716ccba7 17 ,fChi2Penalty(0)
9f42ebf7 18 ,fChi2Match(0)
19 ,fChi2ITSSA(0)
c61e50c3 20 ,fParent(0)
08419930 21#ifdef _ITSU_TUNING_MODE_ // this is used only for tuning histo filling
22 ,fOrdBranch(0)
23 ,fOrdCand(0)
24#endif
32d38de2 25{
26 // def c-tor
44785f3e 27 ResetFMatrix();
32d38de2 28}
29
30//_________________________________________________________________________
31AliITSUSeed::~AliITSUSeed()
32{
716ccba7 33 // d-tor
32d38de2 34}
35
36//_________________________________________________________________________
37AliITSUSeed::AliITSUSeed(const AliITSUSeed& src)
c61e50c3 38 :AliExternalTrackParam(src)
f8832015 39 ,fHitsPattern(src.fHitsPattern)
b8b59e05 40 ,fNChildren(src.fNChildren)
c61e50c3 41 ,fClID(src.fClID)
f8832015 42 ,fChi2Glo(src.fChi2Glo)
43 ,fChi2Cl(src.fChi2Cl)
716ccba7 44 ,fChi2Penalty(src.fChi2Penalty)
9f42ebf7 45 ,fChi2Match(src.fChi2Match)
46 ,fChi2ITSSA(src.fChi2ITSSA)
c61e50c3 47 ,fParent(src.fParent)
08419930 48#ifdef _ITSU_TUNING_MODE_ // this is used only for tuning histo filling
49 ,fOrdBranch(src.fOrdBranch)
50 ,fOrdCand(src.fOrdCand)
51#endif
32d38de2 52{
53 // def c-tor
44785f3e 54 for (int i=kNFElem;i--;) fFMatrix[i] = src.fFMatrix[i];
716ccba7 55 for (int i=kNKElem;i--;) fKMatrix[i] = src.fKMatrix[i];
56 for (int i=kNRElem;i--;) fRMatrix[i] = src.fRMatrix[i];
943e1898 57 fResid[0]=src.fResid[0];
58 fResid[1]=src.fResid[1];
59 fCovIYZ[0]=src.fCovIYZ[0];
60 fCovIYZ[1]=src.fCovIYZ[1];
61 fCovIYZ[2]=src.fCovIYZ[2];
62 //
32d38de2 63}
64
65//_________________________________________________________________________
66AliITSUSeed &AliITSUSeed::operator=(const AliITSUSeed& src)
67{
68 // def c-tor
69 if (this == &src) return *this;
943e1898 70 this->~AliITSUSeed();
71 new(this) AliITSUSeed(src);
32d38de2 72 return *this;
73}
f8832015 74
75//_________________________________________________________________________
76void AliITSUSeed::Print(Option_t* opt) const
77{
78 // print seed info
79 int lr,cl = GetLrCluster(lr);
08419930 80 printf("%cLr%d Nchild: %3d Cl:%4d Chi2Glo:%7.2f(%7.2f) Chi2Cl:%7.2f Penalty: %7.2f Mtc:%6.3f Bwd:%6.3f",IsKilled() ? '-':' ',
9f42ebf7 81 lr,GetNChildren(),cl,GetChi2Glo(),GetChi2GloNrm(),GetChi2Cl(), GetChi2Penalty(), GetChi2ITSTPC(), GetChi2ITSSA());
f8832015 82 printf(" |");
716ccba7 83 int lrc=0;
84 const AliITSUSeed *sdc = this;
85 while(1) {
86 if (lrc<lr) printf(".");
87 else {
88 sdc = sdc->GetParent(lrc);
89 if (!sdc) break;
90 printf("%c",sdc->GetClusterID()<0 ? '.': (sdc->IsFake() ? '-':'+'));
91 }
92 lrc++;
93 }
94 printf("|\n");
f8832015 95 TString opts = opt; opts.ToLower();
96 if (opts.Contains("etp")) AliExternalTrackParam::Print();
97 if (opts.Contains("parent") && GetParent()) GetParent()->Print(opt);
98}
3dd9c283 99
3e4e3c23 100//______________________________________________________________________________
9f42ebf7 101void AliITSUSeed::InitFromSeed(const AliExternalTrackParam* seed)
3e4e3c23 102{
103 // init seed from ESD track
104 TObject::Clear();
9f42ebf7 105 AliExternalTrackParam::operator=(*seed);
3e4e3c23 106 ResetFMatrix();
107 fHitsPattern = 0;
108 fClID = 0;
b8b59e05 109 fNChildren = 0;
3e4e3c23 110 fChi2Glo = fChi2Cl = fChi2Penalty = 0;
111 fParent = 0; //!!!
112}
113
114
3dd9c283 115//______________________________________________________________________________
116Float_t AliITSUSeed::GetChi2GloNrm() const
117{
118 int ndf = 2*GetNLayersHit() - 5;
716ccba7 119 return (ndf>0 ? fChi2Glo/ndf : fChi2Glo) + fChi2Penalty;
3dd9c283 120}
121
122
123//______________________________________________________________________________
124Int_t AliITSUSeed::Compare(const TObject* obj) const
125{
126 // compare clusters accodring to specific mode
f9c7eb32 127 const AliITSUSeed* sd = (const AliITSUSeed*)obj;
3dd9c283 128 const Float_t kTol = 1e-5;
129 if (!IsKilled() && sd->IsKilled()) return -1;
130 if ( IsKilled() &&!sd->IsKilled()) return 1;
131 //
716ccba7 132 float chi2This = GetChi2GloNrm();
133 float chi2Other = sd->GetChi2GloNrm();
134
135 if (chi2This+kTol<chi2Other) return -1;
136 else if (chi2This-kTol>chi2Other) return 1;
3dd9c283 137 return 0;
138}
139
140//______________________________________________________________________________
141Bool_t AliITSUSeed::IsEqual(const TObject* obj) const
142{
143 // compare clusters accodring to specific mode
f9c7eb32 144 const AliITSUSeed* sd = (const AliITSUSeed*)obj;
3dd9c283 145 const Float_t kTol = 1e-5;
146 if (IsKilled() != sd->IsKilled()) return kFALSE;
716ccba7 147 return Abs(GetChi2GloNrm() - sd->GetChi2GloNrm())<kTol;
3dd9c283 148}
44785f3e 149
150//______________________________________________________________________________
151Bool_t AliITSUSeed::PropagateToX(Double_t xk, Double_t b)
152{
153 // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
154 Double_t dx=xk-fX;
155 if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
156
157 Double_t crv=GetC(b);
158 if (TMath::Abs(b) < kAlmost0Field) crv=0.;
159 Double_t x2r = crv*dx;
160 Double_t f1=fP[2], f2=f1 + x2r;
161 if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
162 if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
163 if (TMath::Abs(fP[4])< kAlmost0) return kFALSE;
164
165 Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4];
166 Double_t
167 &fC00=fC[0],
168 &fC10=fC[1], &fC11=fC[2],
169 &fC20=fC[3], &fC21=fC[4], &fC22=fC[5],
170 &fC30=fC[6], &fC31=fC[7], &fC32=fC[8], &fC33=fC[9],
171 &fC40=fC[10], &fC41=fC[11], &fC42=fC[12], &fC43=fC[13], &fC44=fC[14];
172
173 Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
174 if (TMath::Abs(r1)<kAlmost0) return kFALSE;
175 if (TMath::Abs(r2)<kAlmost0) return kFALSE;
176
177 fX=xk;
178 double dy2dx = (f1+f2)/(r1+r2);
179 fP0 += dx*dy2dx;
180 if (TMath::Abs(x2r)<0.05) {
181 fP1 += dx*(r2 + f2*dy2dx)*fP3; // Many thanks to P.Hristov !
182 fP2 += x2r;
183 }
184 else {
185 // for small dx/R the linear apporximation of the arc by the segment is OK,
186 // but at large dx/R the error is very large and leads to incorrect Z propagation
187 // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
188 // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
189 // Similarly, the rotation angle in linear in dx only for dx<<R
190 double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
191 double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
192 fP1 += rot/crv*fP3;
193 fP2 = TMath::Sin(rot + TMath::ASin(fP2));
194 }
195
196 //f = F - 1
197 double r1i = 1./r1;
198 double r2i = 1./r2;
199 double tg1 = f1*r1i;
200 double tg2 = f2*r2i;
201 double v0 = 1. + dy2dx*tg2;
202 double v1 = (r1i+r2i)*(dy2dx*(tg1+tg2)+2);
203 double v2 = (r1i+r2i)*v0;
204 //
205 double f24 = dx*crv/fP4;
206 double f02 = dx*v1;
207 double f04 = dx*v2*f24;
208 double f12 = dx*fP3* (f2*v1+dy2dx-tg2);
209 double f13 = dx*r2*v0;
210 double f14 = dx*f24*fP3*(f2*v2+dy2dx-tg2);
211 //
212 //b = C*ft
213 Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
214 Double_t b02=f24*fC40;
215 Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
216 Double_t b12=f24*fC41;
217 Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
218 Double_t b22=f24*fC42;
219 Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
220 Double_t b42=f24*fC44;
221 Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
222 Double_t b32=f24*fC43;
223
224 //a = f*b = f*C*ft
225 Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
226 Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
227 Double_t a22=f24*b42;
228
229 //F*C*Ft = C + (b + bt + a)
230 fC00 += b00 + b00 + a00;
231 fC10 += b10 + b01 + a01;
232 fC20 += b20 + b02 + a02;
233 fC30 += b30;
234 fC40 += b40;
235 fC11 += b11 + b11 + a11;
236 fC21 += b21 + b12 + a12;
237 fC31 += b31;
238 fC41 += b41;
239 fC22 += b22 + b22 + a22;
240 fC32 += b32;
241 fC42 += b42;
242 //
243 // update stored transformation matrix F = Fnew*Fold
244 fFMatrix[kF04] += f04 + f24*fFMatrix[kF02];
245 fFMatrix[kF14] += f14 + f24*fFMatrix[kF12];
246 fFMatrix[kF02] += f02;
247 fFMatrix[kF12] += f12;
248 fFMatrix[kF13] += f13;
249 fFMatrix[kF24] += f24;
250 //
251 CheckCovariance();
252
253 return kTRUE;
254}
943e1898 255
e7d83d38 256//__________________________________________________________________
257Int_t AliITSUSeed::GetClusterIndex(Int_t ind) const
258{
259 // get ind-th cluster index
260 int ncl = 0;
261 const AliITSUSeed* seed = this;
262 while(seed) {
c425b166 263 if ( seed->HasCluster() && (ncl++==ind) ) return seed->GetLrClusterID();//GetClusterID();
e7d83d38 264 seed = (AliITSUSeed*)seed->GetParent();
265 }
266 return -1;
267 //
268}
269
716ccba7 270//______________________________________________________________________________
271Bool_t AliITSUSeed::RotateToAlpha(Double_t alpha)
272{
273 // Transform this track to the local coord. system rotated
274 // by angle "alpha" (rad) with respect to the global coord. system.
275 //
276 if (TMath::Abs(fP[2]) >= kAlmost1) {
277 AliError(Form("Precondition is not satisfied: |sin(phi)|>1 ! %f",fP[2]));
278 return kFALSE;
279 }
280 //
281 if (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
282 else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
283 //
284 Double_t &fP0=fP[0];
285 Double_t &fP2=fP[2];
286 Double_t &fC00=fC[0];
287 Double_t &fC10=fC[1];
288 Double_t &fC20=fC[3];
289 Double_t &fC21=fC[4];
290 Double_t &fC22=fC[5];
291 Double_t &fC30=fC[6];
292 Double_t &fC32=fC[8];
293 Double_t &fC40=fC[10];
294 Double_t &fC42=fC[12];
295 //
296 Double_t x=fX;
297 Double_t ca=TMath::Cos(alpha-fAlpha), sa=TMath::Sin(alpha-fAlpha);
298 Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); // Improve precision
299 // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
300 // direction in local frame is along the X axis
301 if ((cf*ca+sf*sa)<0) {
302 AliDebug(1,Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
303 return kFALSE;
304 }
305 //
306 Double_t tmp=sf*ca - cf*sa;
307
308 if (TMath::Abs(tmp) >= kAlmost1) {
309 if (TMath::Abs(tmp) > 1.+ Double_t(FLT_EPSILON))
310 AliWarning(Form("Rotation failed ! %.10e",tmp));
311 return kFALSE;
312 }
313 fAlpha = alpha;
314 fX = x*ca + fP0*sa;
315 fP0= -x*sa + fP0*ca;
316 fP2= tmp;
317
318 if (TMath::Abs(cf)<kAlmost0) {
319 AliError(Form("Too small cosine value %f",cf));
320 cf = kAlmost0;
321 }
322
323 Double_t rr=(ca+sf/cf*sa);
324
325 fC00 *= (ca*ca);
326 fC10 *= ca;
327 fC20 *= ca*rr;
328 fC21 *= rr;
329 fC22 *= rr*rr;
330 fC30 *= ca;
331 fC32 *= rr;
332 fC40 *= ca;
333 fC42 *= rr;
334 //
335 fRMatrix[kR00] = ca;
336 fRMatrix[kR22] = rr;
337 //
338 CheckCovariance();
339
340 return kTRUE;
341}
342
943e1898 343//______________________________________________________________________________
344Bool_t AliITSUSeed::GetTrackingXAtXAlpha(double xOther, double alpOther, double bz, double &xdst)
345{
346 // calculate X and Y in the tracking frame of the track, corresponding to other X,Alpha tracking
347 double ca=TMath::Cos(alpOther-fAlpha), sa=TMath::Sin(alpOther-fAlpha);
348 double &y=fP[0], &sf=fP[2], cf=Sqrt((1.-sf)*(1.+sf));
349 double eta = xOther - fX*ca - y*sa;
350 double xi = sf*ca - cf*sa;
716ccba7 351 if (Abs(xi)>= kAlmost1) return kFALSE;
943e1898 352 double nu = xi + GetC(bz)*eta;
716ccba7 353 if (Abs(nu)>= kAlmost1) return kFALSE;
943e1898 354 xdst = xOther*ca - sa*( y*ca-fX*sa + eta*(xi+nu)/(Sqrt((1.-xi)*(1.+xi)) + Sqrt((1.-nu)*(1.+nu))) );
355 return kTRUE;
356}
357
358//____________________________________________________________________
359Double_t AliITSUSeed::GetPredictedChi2(Double_t p[2],Double_t cov[3])
360{
361 // Estimate the chi2 of the space point "p" with the cov. matrix "cov"
362 // Store info needed for update and smoothing
363 Double_t sdd = fC[0] + cov[0];
364 Double_t sdz = fC[1] + cov[1];
365 Double_t szz = fC[2] + cov[2];
366 Double_t det = sdd*szz - sdz*sdz;
367 if (TMath::Abs(det) < kAlmost0) return kVeryBig;
368 det = 1./det;
369 fCovIYZ[0] = szz*det;
370 fCovIYZ[1] = -sdz*det;
371 fCovIYZ[2] = sdd*det;
372 double &dy = fResid[0] = p[0] - fP[0];
373 double &dz = fResid[1] = p[1] - fP[1];
374 //
375 return dy*(dy*fCovIYZ[0]+dz*fCovIYZ[1]) + dz*(dy*fCovIYZ[1]+dz*(fCovIYZ[2]));
376 //
377}
378
379//____________________________________________________________________
380Bool_t AliITSUSeed::Update()
381{
382 // Update the track parameters with the measurement stored during GetPredictedChi2
383 //
716ccba7 384 Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4],
385 &fC00=fC[kS00],
386 &fC10=fC[kS10], &fC11=fC[kS11],
387 &fC20=fC[kS20], &fC21=fC[kS21], &fC22=fC[kS22],
388 &fC30=fC[kS30], &fC31=fC[kS31], &fC32=fC[kS32], &fC33=fC[kS33],
389 &fC40=fC[kS40], &fC41=fC[kS41], &fC42=fC[kS42], &fC43=fC[kS43], &fC44=fC[kS44];
943e1898 390 //
391 double &r00=fCovIYZ[0],&r01=fCovIYZ[1],&r11=fCovIYZ[2];
392 double &dy=fResid[0], &dz=fResid[1];
393 //
716ccba7 394 // store info needed for smoothing in the fKMatrix
395 double &k00 = fKMatrix[kK00] = fC00*r00+fC10*r01;
396 double &k01 = fKMatrix[kK01] = fC00*r01+fC10*r11;
397 double &k10 = fKMatrix[kK10] = fC10*r00+fC11*r01;
398 double &k11 = fKMatrix[kK11] = fC10*r01+fC11*r11;
399 double &k20 = fKMatrix[kK20] = fC20*r00+fC21*r01;
400 double &k21 = fKMatrix[kK21] = fC20*r01+fC21*r11;
401 double &k30 = fKMatrix[kK30] = fC30*r00+fC31*r01;
402 double &k31 = fKMatrix[kK31] = fC30*r01+fC31*r11;
403 double &k40 = fKMatrix[kK40] = fC40*r00+fC41*r01;
404 double &k41 = fKMatrix[kK41] = fC40*r01+fC41*r11;
943e1898 405 //
406 Double_t sf=fP2 + k20*dy + k21*dz;
407 if (TMath::Abs(sf) > kAlmost1) return kFALSE;
408
409 fP0 += k00*dy + k01*dz;
410 fP1 += k10*dy + k11*dz;
411 fP2 = sf;
412 fP3 += k30*dy + k31*dz;
413 fP4 += k40*dy + k41*dz;
414 //
415 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
416 Double_t c12=fC21, c13=fC31, c14=fC41;
417
418 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
419 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
420 fC40-=k00*c04+k01*c14;
421
422 fC11-=k10*c01+k11*fC11;
423 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
424 fC41-=k10*c04+k11*c14;
425
426 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
427 fC42-=k20*c04+k21*c14;
428
429 fC33-=k30*c03+k31*c13;
430 fC43-=k30*c04+k31*c14;
431
432 fC44-=k40*c04+k41*c14;
433 //
943e1898 434 CheckCovariance();
716ccba7 435 //
436 return kTRUE;
437}
438
439
440//____________________________________________________________________
441Bool_t AliITSUSeed::Smooth(Double_t vecL[5],Double_t matL[15])
442{
443 // Prepare MBF smoothing auxiliary params for smoothing at prev. point:
444 // \hat{l_N} = 0
445 // \hat{L_N} = 0
446 // \tilde{l_j} = -H_j^T N_{j}^{-1} z_j + B_{j}^T \hat{l_j}
447 // \tilde{L_j} = H_j^T N_{j}^{-1} H_j + B_j^T \hat{L_j} B_j
448 // \hat{l_j} = F_j^T \tilde{l_{j+1}}
449 // \hat{L_j} = F_j^T \tilde{L_{j+1}} F_j
450 //
451 // P_{j/N} = P_{j/j} - P_{j/j} \hat{L_j} P_{j/j}
452 // \hat{x_{j/N}} = \hat{x_{j/j}} - P_{j/j} \hat{l_j}
453 //
454 // N^-1 = fCovIYZ
455 // z = fResid
456 // B = I - K H
457 // H = {{1,0,0,0,0},{0,1,0,0,0}}
458 //
459 // calc. \tilde{l_j}
460 //
461 if (GetClusterID()<0) return kTRUE;
462 //
463
464 double
465 &k00=fKMatrix[kK00],&k01=fKMatrix[kK01],
466 &k10=fKMatrix[kK10],&k11=fKMatrix[kK11],
467 &k20=fKMatrix[kK20],&k21=fKMatrix[kK21],
468 &k30=fKMatrix[kK30],&k31=fKMatrix[kK31],
469 &k40=fKMatrix[kK40],&k41=fKMatrix[kK41];
470 double
471 &l00=matL[kS00],
472 &l10=matL[kS10], &l11=matL[kS11],
473 &l20=matL[kS20], &l21=matL[kS21], &l22=matL[kS22],
474 &l30=matL[kS30], &l31=matL[kS31], &l32=matL[kS32], &l33=matL[kS33],
475 &l40=matL[kS40], &l41=matL[kS41], &l42=matL[kS42], &l43=matL[kS43], &l44=matL[kS44];
476 //
477 // calculate correction
478 double corrVec[5]={0},corrMat[15]={0};
479 corrVec[0] = fC[kS00]*vecL[0] + fC[kS10]*vecL[1] + fC[kS20]*vecL[2] + fC[kS30]*vecL[3] + fC[kS40]*vecL[4];
480 corrVec[1] = fC[kS10]*vecL[0] + fC[kS11]*vecL[1] + fC[kS21]*vecL[2] + fC[kS31]*vecL[3] + fC[kS41]*vecL[4];
481 corrVec[2] = fC[kS20]*vecL[0] + fC[kS21]*vecL[1] + fC[kS22]*vecL[2] + fC[kS32]*vecL[3] + fC[kS42]*vecL[4];
482 corrVec[3] = fC[kS30]*vecL[0] + fC[kS31]*vecL[1] + fC[kS32]*vecL[2] + fC[kS33]*vecL[3] + fC[kS43]*vecL[4];
483 corrVec[4] = fC[kS40]*vecL[0] + fC[kS41]*vecL[1] + fC[kS42]*vecL[2] + fC[kS43]*vecL[3] + fC[kS44]*vecL[4];
484 //
485 double *crm = ProdABA(fC,matL);
486 for (int i=0;i<15;i++) corrMat[i] = crm[i];
487
488 double vcL0 = vecL[0], vcL1 = vecL[1];
489 vecL[0] -= k00*vcL0+k10*vcL1+k20*vecL[2]+k30*vecL[3]+k40*vecL[4] + fCovIYZ[0]*fResid[0] + fCovIYZ[1]*fResid[1];
490 vecL[1] -= k01*vcL0+k11*vcL1+k21*vecL[2]+k31*vecL[3]+k41*vecL[4] + fCovIYZ[1]*fResid[0] + fCovIYZ[2]*fResid[1];
491
492 /*
493 double vcL0 = vecL[0], vcL1 = vecL[1];
494 vecL[0] -= k00*vcL0+k10*vcL1+fKMatrix[kK20]*vecL[2]+k30*vecL[3]+k40*vecL[4] + fCovIYZ[0]*fResid[0] + fCovIYZ[1]*fResid[1];
495 vecL[1] -= k01*vcL0+fKMatrix[kK11]*vcL1+k21*vecL[2]+k31*vecL[3]+k41*vecL[4] + fCovIYZ[1]*fResid[0] + fCovIYZ[2]*fResid[1];
496 vecL[3] += fFMatrix[kF13]*vecL[1];
497 vecL[4] = fFMatrix[kF04]*vecL[0] + fFMatrix[kF14]*vecL[1] + fFMatrix[kF24]*vecL[2] + fFMatrix[kF44]*vecL[4];
498 vecL[2] += fFMatrix[kF02]*vecL[0] + fFMatrix[kF12]*vecL[1];
499 //
500 */
501 // and \hat{l_j} in one go
502
503 // L = H^T * sg * H + (I-KH)^T * L * (I - KH)
504 double v00 = k00*l00+k10*l10+k20*l20+k30*l30+k40*l40;
505 double v10 = k00*l10+k10*l11+k20*l21+k30*l31+k40*l41;
506 double v20 = k00*l20+k10*l21+k20*l22+k30*l32+k40*l42;
507 double v30 = k00*l30+k10*l31+k20*l32+k30*l33+k40*l43;
508 double v40 = k00*l40+k10*l41+k20*l42+k30*l43+k40*l44;
509 //
510 double v01 = k01*l00+k11*l10+k21*l20+k31*l30+k41*l40;
511 double v11 = k01*l10+k11*l11+k21*l21+k31*l31+k41*l41;
512 double v21 = k01*l20+k11*l21+k21*l22+k31*l32+k41*l42;
513 double v31 = k01*l30+k11*l31+k21*l32+k31*l33+k41*l43;
514 double v41 = k01*l40+k11*l41+k21*l42+k31*l43+k41*l44;
515 //
516 // (H^T * K^T * L * K * H) - (L * K * H) - (H^T * K^T * L) + (H^T*N^-1*H)
517 l00 += k00*v00 + k10*v10 + k20*v20 + k30*v30 + k40*v40 - v00 - v00 + fCovIYZ[0];
518 l10 += k01*v00 + k11*v10 + k21*v20 + k31*v30 + k41*v40 - v01 - v10 + fCovIYZ[1];
519 l11 += k01*v01 + k11*v11 + k21*v21 + k31*v31 + k41*v41 - v11 - v11 + fCovIYZ[2];
520 //
521 l20 -= v20;
522 l21 -= v21;
523 l30 -= v30;
524 l31 -= v31;
525 l40 -= v40;
526 l41 -= v41;
527 //
528 printf("CorrMt:\n");
529 printf("%+e\n%+e %+e\n%+e %+e %+e\n%+e %+e %+e %+e\n%+e %+e %+e %+e %+e\n",
530 corrMat[kS00],corrMat[kS10],corrMat[kS11],corrMat[kS20],corrMat[kS21],corrMat[kS22],
531 corrMat[kS30],corrMat[kS31],corrMat[kS32],corrMat[kS33],
532 corrMat[kS40],corrMat[kS41],corrMat[kS42],corrMat[kS43],corrMat[kS44]);
533
534 printf("SMcorr: %+e %+e %+e %+e %+e\n",corrVec[0],corrVec[1],corrVec[2],corrVec[3],corrVec[4]);
535
536 printf("State : "); this->AliExternalTrackParam::Print("");
537 //
538 printf("\nBefore transport back (RotElems: %+e %+e)\n",fRMatrix[kR00],fRMatrix[kR22]);
539 printf("Res: %+e %+e | Err: %+e %+e %+e\n",fResid[0],fResid[1],fCovIYZ[0],fCovIYZ[1],fCovIYZ[2]);
540 printf("Lr%d VecL: ",GetLayerID()); for (int i=0;i<5;i++) printf("%+e ",vecL[i]); printf("\n");
541 //
542 printf("%+e\n%+e %+e\n%+e %+e %+e\n%+e %+e %+e %+e\n%+e %+e %+e %+e %+e\n",
543 matL[kS00],matL[kS10],matL[kS11],matL[kS20],matL[kS21],matL[kS22],
544 matL[kS30],matL[kS31],matL[kS32],matL[kS33],matL[kS40],matL[kS41],matL[kS42],matL[kS43],matL[kS44]);
545 //
546 printf("F: "); for (int i=0;i<kNFElem;i++) printf("%+e ",fFMatrix[i]); printf("\n");
547 printf("K: "); for (int i=0;i<kNKElem;i++) printf("%+e ",fKMatrix[i]); printf("\n");
548 //
549 // apply rotation matrix (diagonal)
550 vecL[0] *= fRMatrix[kR00];
551 vecL[2] *= fRMatrix[kR22];
552 //
553 l00 *= fRMatrix[kR00]*fRMatrix[kR00];
554 l10 *= fRMatrix[kR00];
555 l20 *= fRMatrix[kR22]*fRMatrix[kR00];
556 l21 *= fRMatrix[kR22];
557 l22 *= fRMatrix[kR22]*fRMatrix[kR22];
558 l30 *= fRMatrix[kR00];
559 l32 *= fRMatrix[kR22];
560 l40 *= fRMatrix[kR00];
561 l42 *= fRMatrix[kR22];
562 //
563 // Apply translation matrix F^T. Note, that fFMatrix keeps non-trivial elems of F-1 = f, except the e-loss coeff f44
564 // We need F^T*L* F = L + (L*f) + (L*f)^T + f^T * (L*f)
565 //
566 double
567 &f02=fFMatrix[kF02],&f04=fFMatrix[kF04],
568 &f12=fFMatrix[kF12],&f13=fFMatrix[kF13],&f14=fFMatrix[kF14],
569 &f24=fFMatrix[kF24],
570 f44 =fFMatrix[kF44];
571 //
572 vecL[4] = f04*vecL[0]+f14*vecL[1]+f24*vecL[2]+f44*vecL[4];
573 vecL[3] += f13*vecL[1];
574 vecL[2] += f02*vecL[0]+f12*vecL[1];
575 //
576 f44 -= 1.0; // !!!!!
577 //
578 //b = L*f
579 Double_t b02=l00*f02+l10*f12, b03=l10*f13, b04=l00*f04+l10*f14+l20*f24+l40*f44;
580 Double_t b12=l10*f02+l11*f12, b13=l11*f13, b14=l10*f04+l11*f14+l21*f24+l41*f44;
581 Double_t b22=l20*f02+l21*f12, b23=l21*f13, b24=l20*f04+l21*f14+l22*f24+l42*f44;
582 Double_t b32=l30*f02+l31*f12, b33=l31*f13, b34=l30*f04+l31*f14+l32*f24+l43*f44;
583 Double_t b42=l40*f02+l41*f12, b43=l41*f13, b44=l40*f04+l41*f14+l42*f24+l44*f44;
584 //
585 //a = f^T * b = f^T * L * f, profit from symmetry
586 Double_t a22=f02*b02+f12*b12, a33=f13*b13, a44=f04*b04+f14*b14+f24*b24+f44*b44,
587 a32=f13*b12, //= a23=f02*b03+f12*b13,
588 a42=f02*b04+f12*b14, //f04*b02+f14*b12+f24*b22+f44*b42 = a24
589 a43=f13*b14; //f04*b03+f14*b13+f24*b23+f44*b43 = a34
590 //
591 // F^T*L* F = L + (b + b^T + a)
592 l44 += b44 + b44 + a44;
593 l43 += b43 + b34 + a43;
594 l42 += b42 + b24 + a42;
595 l41 += b14;
596 l40 += b04;
597 l33 += b33 + b33 + a33;
598 l32 += b32 + b23 + a32;
599 l31 += b13;
600 l30 += b03;
601 l22 += b22 + b23 + a22;
602 l21 += b12;
603 l20 += b02;
604 //
605 printf("After transport back\n");
606 printf("Lr%d VecL: ",GetLayerID()); for (int i=0;i<5;i++) printf("%+e ",vecL[i]); printf("\n");
607 //
608 printf("%+e\n%+e %+e\n%+e %+e %+e\n%+e %+e %+e %+e\n%+e %+e %+e %+e %+e\n",
609 matL[kS00],matL[kS10],matL[kS11],matL[kS20],matL[kS21],matL[kS22],
610 matL[kS30],matL[kS31],matL[kS32],matL[kS33],matL[kS40],matL[kS41],matL[kS42],matL[kS43],matL[kS44]);
611
612 return kTRUE;
613}
614
615//____________________________________________________________________
616Double_t* AliITSUSeed::ProdABA(const double a[15],const double b[15]) const
617{
618 // product of symmetric matrices A*B*A
619 //
620 const Short_t knd[5][5] = {
621 {kS00,kS10,kS20,kS30,kS40},
622 {kS10,kS11,kS21,kS31,kS41},
623 {kS20,kS21,kS22,kS32,kS42},
624 {kS30,kS31,kS32,kS33,kS43},
625 {kS40,kS41,kS42,kS43,kS44}
626 };
627 //
628 static double aba[15];
629 // 1) ba = B*A
630 double ba[5][5];
631 for (int i=5;i--;) for (int j=5;j--;) {
632 ba[i][j] = 0;
633 for (int k=5;k--;) ba[i][j] += b[knd[i][k]]*a[knd[k][j]];
634 }
635 //
636 // 2) A * ba, lower triangle only
637 for (int i=5;i--;) for (int j=i+1;j--;) {
638 aba[knd[i][j]] = 0;
639 for (int k=5;k--;) aba[knd[i][j]] += a[knd[i][k]]*ba[k][j];
640 }
641 //
642 return &aba[0];
643}
943e1898 644
38997c3c 645//____________________________________________________________________
646Bool_t AliITSUSeed::ContainsFake() const
647{
648 // check if the full branch containes a fake cluster
649 const AliITSUSeed* seed = this;
650 while(seed) {
651 if ( seed->IsFake() ) return kTRUE;
652 seed = (AliITSUSeed*)seed->GetParent();
653 }
654 return kFALSE;
655}
656
9f42ebf7 657//__________________________________________________________________
658Int_t AliITSUSeed::FetchClusterInfo(Int_t *clIDarr) const
659{
660 // fill cl.id's in the array. The clusters of layer L will be set at slots
661 // clID[2L] (and clID[2L+1] if there is an extra cluster).
662 Int_t lr,ncl=0;
663 const AliITSUSeed* seed = this;
664 do {
665 int clID = seed->GetLrCluster(lr);
666 if (clID>=0) {
667 lr<<=1;
668 clIDarr[ clIDarr[lr]<0 ? lr : lr+1 ] = clID;
669 ncl++;
670 }
671 } while ((seed=(AliITSUSeed*)seed->GetParent()));
672 return ncl;
673}
674
716ccba7 675/*
676//____________________________________________________________________
677Bool_t AliITSUSeed::Smooth(Double_t vecL[5],Double_t matL[15])
678{
679 // Prepare MBF smoothing auxiliary params for smoothing at prev. point:
680 // \hat{l_N} = 0
681 // \hat{L_N} = 0
682 // \tilde{l_j} = -H_j^T N_{j}^{-1} z_j + B_{j}^T \hat{l_j}
683 // \tilde{L_j} = H_j^T N_{j}^{-1} H_j + B_j^T \hat{L_j} B_j
684 // \hat{l_j} = F_j^T \tilde{l_{j+1}}
685 // \hat{L_j} = F_j^T \tilde{L_{j+1}} F_j
686 //
687 // P_{j/N} = P_{j/j} - P_{j/j} \hat{L_j} P_{j/j}
688 // \hat{x_{j/N}} = \hat{x_{j/j}} - P_{j/j} \hat{l_j}
689 //
690 // N^-1 = fCovIYZ
691 // z = fResid
692 // B = I - K H
693 // H = {{1,0,0,0,0},{0,1,0,0,0}}
694 //
695 // calc. \tilde{l_j} and \hat{l_j} in one go
696 //
697 if (GetClusterID()<0) return kTRUE;
698 //
699 double
700 &k00=fKMatrix[kK00],&k01=fKMatrix[kK01],
701 &k10=fKMatrix[kK10],&k11=fKMatrix[kK11],
702 &k20=fKMatrix[kK20],&k21=fKMatrix[kK21],
703 &k30=fKMatrix[kK30],&k31=fKMatrix[kK31],
704 &k40=fKMatrix[kK40],&k41=fKMatrix[kK41];
705 double
706 &matL00=matL[kS00],
707 &matL10=matL[kS01], &matL11=matL[kS11],
708 &matL20=matL[kS20], &matL21=matL[kS21], &matL22=matL[kS22],
709 &matL30=matL[kS30], &matL31=matL[kS31], &matL32=matL[kS32], &matL33=matL[kS33],
710 &matL40=matL[kS40], &matL41=matL[kS41], &matL42=matL[kS42], &matL43=matL[kS43], &matL44=matL[kS44];
711 //
712 double vcL0 = vecL[0], vcL1 = vecL[1];
713 vecL[0] -= k00*vcL0+k10*vcL1+fKMatrix[kK20]*vecL[2]+k30*vecL[3]+k40*vecL[4] + fCovIYZ[0]*fResid[0] + fCovIYZ[1]*fResid[1];
714 vecL[1] -= k01*vcL0+fKMatrix[kK11]*vcL1+k21*vecL[2]+k31*vecL[3]+k41*vecL[4] + fCovIYZ[1]*fResid[0] + fCovIYZ[2]*fResid[1];
715 vecL[3] += fFMatrix[kF13]*vecL[1];
716 vecL[4] = fFMatrix[kF04]*vecL[0] + fFMatrix[kF14]*vecL[1] + fFMatrix[kF24]*vecL[2] + fFMatrix[kF44]*vecL[4];
717 vecL[2] += fFMatrix[kF02]*vecL[0] + fFMatrix[kF12]*vecL[1];
718 //
719
720 // L = H^T * sg * H + (I-KH)^T * L * (I - KH)
721 double v00 = k00*matL00+k10*matL10+k20*matL20+k30*matL30+k40*matL40;
722 double v10 = k00*matL10+k10*matL11+k20*matL21+k30*matL31+k40*matL41;
723 double v20 = k00*matL20+k10*matL12+k20*matL22+k30*matL32+k40*matL42;
724 double v30 = k00*matL30+k10*matL13+k20*matL23+k30*matL33+k40*matL43;
725 double v40 = k00*matL40+k10*matL14+k20*matL24+k30*matL34+k40*matL44;
726 //
727 double v01 = k01*matL00+k11*matL10+k21*matL20+k31*matL30+k41*matL40;
728 double v11 = k01*matL01+k11*matL11+k21*matL21+k31*matL31+k41*matL41;
729 double v21 = k01*matL02+k11*matL12+k21*matL22+k31*matL32+k41*matL42;
730 double v31 = k01*matL03+k11*matL13+k21*matL23+k31*matL33+k41*matL43;
731 double v41 = k01*matL04+k11*matL14+k21*matL24+k31*matL34+k41*matL44;
732 //
733 double t00 = k00*matL00+k10*matL01+k20*matL02+k30*matL03+k40*matL04;
734 double t10 = k00*matL10+k10*matL11+k20*matL12+k30*matL13+k40*matL14;
735 double t20 = k00*matL20+k10*matL21+k20*matL22+k30*matL23+k40*matL24;
736 double t30 = k00*matL30+k10*matL31+k20*matL32+k30*matL33+k40*matL34;
737 double t40 = k00*matL40+k10*matL41+k20*matL42+k30*matL43+k40*matL44;
738 //
739 double t01 = k01*matL00+k11*matL01+k21*matL02+k31*matL03+k41*matL04;
740 double t11 = k01*matL10+k11*matL11+k21*matL12+k31*matL13+k41*matL14;
741 double t21 = k01*matL20+k11*matL21+k21*matL22+k31*matL23+k41*matL24;
742 double t31 = k01*matL30+k11*matL31+k21*matL32+k31*matL33+k41*matL34;
743 double t41 = k01*matL40+k11*matL41+k21*matL42+k31*matL43+k41*matL44;
744 //
745 // (H^T * K^T * L * K * H) - (L * K * H) - (H^T * K^T * L) + (H^T*N^-1*H)
746 matL00 += k00*v00+k10*v10+k20*v20*k30*v30+k40*v40 - t00 - v00 + fCovIYZ[0];
747 matL01 += k01*v00+k11*v10+k21*v20*k31*v30+k41*v40 - t01 - v10 + fCovIYZ[1];
748 matL10 += k00*v01+k10*v11+k20*v21*k30*v31+k40*v41 - t10 - v01 + fCovIYZ[1];
749 matL11 += k01*v01+k11*v11+k21*v21*k31*v31+k41*v41 - t11 - v11 + fCovIYZ[2];
750 //
751 matL20 -= t20;
752 matL21 -= t21;
753 matL30 -= t30;
754 matL31 -= t31;
755 matL40 -= t40;
756 matL41 -= t41;
757 //
758 matL02 -= v20;
759 matL03 -= v30;
760 matL04 -= v40;
761 matL12 -= v21;
762 matL13 -= v31;
763 matL14 -= v41;
764 //
765 printf("Lr%d VecL: ",GetLayerID()); for (int i=0;i<5;i++) printf("%+e ",vecL[i]); printf("\n");
766 printf("F: "); for (int i=0;i<kNFElem;i++) printf("%+e ",fFMatrix[i]); printf("\n");
767 printf("K: "); for (int i=0;i<kNKElem;i++) printf("%+e ",fKMatrix[i]); printf("\n");
768 //
769 for (int j=0;j<5;j++) {
770 for (int i=0;i<5;i++) printf("%+e ",matL[j][i]); printf("\n");
771 }
772 //
943e1898 773 return kTRUE;
774}
716ccba7 775
776 */
e7d83d38 777