FairRoot/PandaRoot
PndFTSCATrackParamVector.h
Go to the documentation of this file.
1 //-*- Mode: C++ -*-
2 // $Id: PndFTSCATrackParamVector.h,v 1.3 2011/01/31 17:18:32 fisyak Exp $
3 // ************************************************************************
4 // This file is property of and copyright by the ALICE HLT Project *
5 // ALICE Experiment at CERN, All rights reserved. *
6 // See cxx source for full Copyright notice *
7 // *
8 //*************************************************************************
9 
10 
11 #ifndef PNDFTSCATRACKPARAMVECTOR_H
12 #define PNDFTSCATRACKPARAMVECTOR_H
13 
14 class PndFTSCATrackParam;
15 
16 #ifdef PANDA_FTS
17 
18 #include "PndFTSCADef.h"
19 #include "PndFTSVector.h"
20 #include "PndFTSCAMath.h"
21 #include "L1XYMeasurementInfo.h"
22 #include "L1MaterialInfo.h"
23 #include "L1Field.h"
24 
25 #include "FairField.h"
26 #include "FairRunAna.h"
27 // #include "CAX1X2MeasurementInfo.h"
28 #include "FTSCAStation.h"
29 
30 class PndFTSCAParam;
31 class FTSCAHit;
32 class FTSCAHitV;
33 class FTSCATarget;
34 
42 class PndFTSCATrackParamVector
43 {
44  friend std::istream &operator>>( std::istream &, PndFTSCATrackParamVector & );
45  friend std::ostream &operator<<( std::ostream &, const PndFTSCATrackParamVector & );
46  public:
47  PndFTSCATrackParamVector():
48  fX(0.f),fY(0.f),fTx(0.f),fTy(0.f),fQP(0.f),fZ(0.f),
49  fC00(0.f),
50  fC10(0.f), fC11(0.f),
51  fC20(0.f), fC21(0.f), fC22(0.f),
52  fC30(0.f), fC31(0.f), fC32(0.f), fC33(0.f),
53  fC40(0.f), fC41(0.f), fC42(0.f), fC43(0.f), fC44(0.f),
54  fChi2(0.f), fNDF(0.f)
55  { fZB[0] = fZB[1] = 10e10; fDirection = true; }
56 
57  void SetTrackParam(const PndFTSCATrackParamVector &param, const float_m &m = float_m( true ) )
58  {
59  for(int i=0; i<5; i++) Par(i)(m) = param.Par(i);
60  for(int i=0; i<15; i++) Cov(i)(m) = param.Cov(i);
61  fX(m) = param.X();
62  fChi2(m) = param.Chi2();
63  fZ(m) = param.Z();
64  fNDF(static_cast<int_m>(m)) = param.NDF();
65  fAlpha(static_cast<int_m>(m)) = param.Angle();
66  }
67 
68  void SetTrackParamOne(int iV, const PndFTSCATrackParamVector &param, int iVa )
69  {
70  for(int i=0; i<5; i++) Par(i)[iV] = param.Par(i)[iVa];
71  for(int i=0; i<15; i++) Cov(i)[iV] = param.Cov(i)[iVa];
72  fX[iV] = param.X()[iVa];
73  fZ[iV] = param.Z()[iVa];
74  fChi2[iV] = param.Chi2()[iVa];
75  fNDF[iV] = param.NDF()[iVa];
76  fAlpha[iV] = param.Angle()[iVa];
77  }
78 
79  void ConvertTrackParamToVector( PndFTSCATrackParam t0[float_v::Size], int nTracksV );
80 
81  void PrintCovMat()
82  {
83  cout << "CovMat " << endl;
84  //for ( int i = 0; i < 15; i++ )
85  cout << Cov(0) << endl;
86  cout << Cov(1) << " "<< Cov(2) << endl;
87  cout << Cov(3) << " "<< Cov(4) <<" "<< Cov(5) << endl;
88  cout << Cov(6) << " "<< Cov(7) <<" "<< Cov(8) <<" "<< Cov(9) << endl;
89  cout << Cov(10) <<" "<< Cov(11) <<" "<< Cov(12) <<" "<< Cov(13) <<" "<< Cov(14) << endl;
90  /*cout << Cov(0)[0] << endl;
91  cout << Cov(1)[0] << " "<< Cov(2)[0] << endl;
92  cout << Cov(3)[0] << " "<< Cov(4)[0] <<" "<< Cov(5)[0] << endl;
93  cout << Cov(6)[0] << " "<< Cov(7)[0] <<" "<< Cov(8)[0] <<" "<< Cov(9)[0] << endl;
94  cout << Cov(10)[0] <<" "<< Cov(11)[0] <<" "<< Cov(12)[0] <<" "<< Cov(13)[0] <<" "<< Cov(14)[0];*/
95  }
96 
97  void InitCovMatrix( float_v d2QMom = 0.f );
98  void InitByTarget( const FTSCATarget& target );
99  void InitByHit( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_v& dQP );
100 
101  void InitDirection( float_v r0, float_v r1, float_v r2 ) { // initialize direction parameters according to a given tangent vector r = { r0, r1, r2 }
102  float_m mask = (r0 > 0.f);
103  float_v tx (Vc::Zero);
104  tx(mask) = r1/r0;
105  float_v ty(Vc::Zero);
106  ty(mask) = r2/r0;
107  SetTx( tx );
108  SetTy( ty );
109  }
110 
111 
112  float_v X() const { return fX; }
113  float_v Y() const { return fY; }
114  float_v Z() const { return fZ; }
115  float_v Tx() const { return fTx; }
116  float_v Ty() const { return fTy; }
117  float_v QP() const { return fQP;}
118 
119  float_v Bx() const { return fBx; }
120  float_v By() const { return fBy; }
121 
122 
123  float_v Chi2() const { return fChi2; }
124  int_v NDF() const { return fNDF; }
125 
126  float_v Angle() const { return fAlpha; }
127 
128  float_v X0() const { return Z(); }
129  float_v X1() const { return X(); }
130  float_v X2() const { return Y(); }
131  float_v Tx1() const { return Tx(); }
132  float_v Tx2() const { return Ty(); }
133  float_v QMomentum() const { return QP(); } // used for triplets comparison
134  float_v SinPhi() const { // dx/dxz
135  const float_v tx12 = Tx1()*Tx1();
136  return sqrt( tx12/(1.f + tx12) );
137  }
138 
139  const float_v& Par( int i ) const {
140  switch ( i ) {
141  case 0: return fX;
142  case 1: return fY;
143  case 2: return fTx;
144  case 3: return fTy;
145  case 4: return fQP;
146  case 5: return fZ;
147  };
148  assert(0); return fX;
149  }
150 
151  const float_v& Cov( int i ) const {
152  switch ( i ) {
153  case 0: return fC00;
154  case 1: return fC10;
155  case 2: return fC11;
156  case 3: return fC20;
157  case 4: return fC21;
158  case 5: return fC22;
159  case 6: return fC30;
160  case 7: return fC31;
161  case 8: return fC32;
162  case 9: return fC33;
163  case 10: return fC40;
164  case 11: return fC41;
165  case 12: return fC42;
166  case 13: return fC43;
167  case 14: return fC44;
168  };
169  assert(0); return fC00;
170  }
171 
172  float_v Err2X1() const { return Err2X(); }
173  float_v Err2X2() const { return Err2Y(); }
174 
175  float_v Err2X() const { return fC00; }
176  float_v Err2Y() const { return fC11; }
177  float_v Err2QP() const { return fC44; }
178  float_v Err2QMomentum() const { return Err2QP(); }
179 
180 
181 
182  void SetX( const float_v &v ) { fX = v; }
183  void SetY( const float_v &v ) { fY = v; }
184  void SetZ( const float_v &v ) { fZ = v; }
185  void SetTx( const float_v &v ) { fTx = v; }
186  void SetTy( const float_v &v ) { fTy = v; }
187  void SetQP( const float_v &v ) { fQP = v; }
188  void SetQMomentum( const float_v &v ) { fQP = v; }
189  void SetBx( const float_v &v ) { fBx = v; }
190  void SetBy( const float_v &v ) { fBy = v; }
191 
192  void SetChi2( const float_v &v ) { fChi2 = v; }
193  void SetNDF( int v ) { fNDF = v; }
194  void SetNDF( const int_v &v ) { fNDF = v; }
195 
196  void SetAngle( const float_v& v ) { fAlpha = v; }
197 
198  void SetPar( int i, float_v v ) {
199  switch ( i ) {
200  case 0: fX = v; break;
201  case 1: fY = v; break;
202  case 2: fTx = v; break;
203  case 3: fTy = v; break;
204  case 4: fQP = v; break;
205  }
206  }
207  void SetCov( int i, float_v v ) {
208  switch ( i ) {
209  case 0: fC00 = v; break;
210  case 1: fC10 = v; break;
211  case 2: fC11 = v; break;
212  case 3: fC20 = v; break;
213  case 4: fC21 = v; break;
214  case 5: fC22 = v; break;
215  case 6: fC30 = v; break;
216  case 7: fC31 = v; break;
217  case 8: fC32 = v; break;
218  case 9: fC33 = v; break;
219  case 10: fC40 = v; break;
220  case 11: fC41 = v; break;
221  case 12: fC42 = v; break;
222  case 13: fC43 = v; break;
223  case 14: fC44 = v; break;
224  }
225  }
226 
227  void SetDirection( bool b ) { fDirection = b; }
228 
229  float_v& Par( int i ) {
230  switch ( i ) {
231  case 0: return fX;
232  case 1: return fY;
233  case 2: return fTx;
234  case 3: return fTy;
235  case 4: return fQP;
236  };
237  assert(0); return fX;
238  }
239 
240  float_v& Cov( int i ) {
241  switch ( i ) {
242  case 0: return fC00;
243  case 1: return fC10;
244  case 2: return fC11;
245  case 3: return fC20;
246  case 4: return fC21;
247  case 5: return fC22;
248  case 6: return fC30;
249  case 7: return fC31;
250  case 8: return fC32;
251  case 9: return fC33;
252  case 10: return fC40;
253  case 11: return fC41;
254  case 12: return fC42;
255  case 13: return fC43;
256  case 14: return fC44;
257  };
258  assert(0); return fC00;
259  }
260 
261  void SetX( const float_v &v, const float_m &m ) { fX( m ) = v; }
262  void SetY( const float_v &v, const float_m &m ) { fY( m ) = v; }
263  void SetZ( const float_v &v, const float_m &m ) { fZ( m ) = v; }
264  void SetTx( const float_v &v, const float_m &m ) { fTx( m ) = v; }
265  void SetTy( const float_v &v, const float_m &m ) { fTy( m ) = v; }
266  void SetQP( const float_v &v, const float_m &m ) { fQP( m ) = v; }
267  void SetQMomentum( const float_v &v, const float_m &m ) { fQP( m ) = v; }
268 
269  void SetChi2( const float_v &v, const float_m &m ) { fChi2(m) = v; }
270  void SetNDF( const int_v &v, const int_m &m ) { fNDF(m) = v; }
271 
272  void SetField( int i, const CAFieldValue& b, const float_v& zb ) { fB[i] = b; fZB[i] = zb; }
273 
274  void UpdateFieldValues(const FTSCAHitV& hit, const PndFTSCAParam& param, L1FieldRegion& f, const float_m& mask);
275  void UpdateFieldValues(const FTSCAHitV& hit, int_v& iVrt, float_v& zVirtualStation, const PndFTSCAParam& param, L1FieldRegion& f, const float_m& mask);
276 
277  void SetCovX12( float_v v00, float_v v10, float_v v11 ) { fC00 = v00; fC10 = v10; fC11 = v11; }
278 
279 // float_m Transport(const int_v& ista, const int_v& iVSta, const PndFTSCAParam& param, const float_m&mask = float_m( true ) );
280  //float_m Transport( const int_v& ista, const PndFTSCAParam& param, const float_m &mask = float_m( true ) );
281 
282  float_m Transport( const FTSCAHitV& hit, const PndFTSCAParam& p, float_v& qp0, const float_m &mask = float_m( true ) );
283  float_m Transport( const FTSCAHitV& hit, const PndFTSCAParam& p, const float_m &mask = float_m( true ) );
284 
285  float_m TransportByLine( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_m &mask = float_m( true ) );
286 
287  float_m Filter( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_m &mask = float_m( true ), const float_v& chi2Cut = 10e10f );
288 
289  float_m Transport( const FTSCAHit& hit, const PndFTSCAParam& p, const float_m &mask = float_m( true ) );
290  float_m Filter( const FTSCAHit& hit, const PndFTSCAParam& param, const float_m &mask = float_m( true ), const float_v& chi2Cut = 10e10f );
291 
292  float_m AddTarget( const FTSCATarget& t, const float_m &mask = float_m( true ) );
293 
294  private:
295 // float_m Transport( const FTSCAHitV& hit, const L1FieldRegion& F, const PndFTSCAParam& p, const float_m &mask = float_m( true ) ); // dbg TODO delme
296 
297  float_m TransportToX0WithMaterial( const float_v &x0, const L1FieldRegion &F, const L1MaterialInfo &material, float_v &qp0, const float_m &mask = float_m( true ) );
298  float_m TransportToX0( const float_v &x0, const L1FieldRegion &F, const float_v &qp0, const float_m &mask = float_m( true ) );
299  float_m RK4TransportToX0( const float_v &x0_out, const L1FieldRegion &F, const float_v &qp0, const float_m &mask = float_m(true) );
300  float_m TransportToX0Line( const float_v &x0, const float_m &mask = float_m( true ) );
301  float_m PassMaterial( const L1MaterialInfo &info, const float_v &qp0, const float_m &mask = float_m( true ) );
302  void EnergyLossCorrection( const float_v& mass2, const float_v& radThick, float_v &qp0, float_v direction, const float_m &mask);
303  float_v ApproximateBetheBloch( const float_v &bg2 );
304 
305  float_m Filter( const float_v &y0, const float_v &z0, const float_v &r, const FTSCAStripInfoVector &info, float_v err2, const float_m &active, const float_v& chi2Cut = 10e10f );
306 
307  float_m FilterVtx( const float_v &xV, const float_v& yV, const L1XYMeasurementInfo &info, float_v& extrDx, float_v& extrDy, float_v J[], const float_m &active = float_m( true ) );
308 
309  float_m TransportJXY0ToX0( const float_v &x0, const L1FieldRegion &F, float_v& extrDx, float_v& extrDy, float_v &J04, float_v &J14, const float_m &active = float_m( true ) );
310 
311 
312  float_v fX,fY,fTx,fTy,fQP,fZ,
313  fC00,
314  fC10, fC11,
315  fC20, fC21, fC22,
316  fC30, fC31, fC32, fC33,
317  fC40, fC41, fC42, fC43, fC44;
318 
319  float_v fChi2; // the chi^2 value
320  int_v fNDF; // the Number of Degrees of Freedom
321 
322  float_v fBx;
323  float_v fBy;
324  CAFieldValue fB[2]; // field at two previous points, which track passed
325  float_v fZB[2]; // fZ position of the filld point
326  bool fDirection;
327  float_v fAlpha; // coor system
328  float_v fDelta;
329 };
330 
331 #include "debug.h"
332 
333 inline float_m PndFTSCATrackParamVector::TransportToX0Line( const float_v &x0_out, const float_m &mask )
334 {
335  float_v dz = (x0_out - fZ);
336  //cout<<"dz "<<dz<<endl;
337 
338  fX(mask) += fTx*dz;
339  fY(mask) += fTy*dz;
340  fZ(mask) += dz;
341 
342  const float_v dzC32_in = dz * fC32;
343 
344  fC21(mask) += dzC32_in;
345  fC10(mask) += dz * ( fC21 + fC30 );
346 
347  const float_v C20_in = fC20;
348 
349  fC20(mask) += dz * fC22;
350  fC00(mask) += dz * ( fC20 + C20_in );
351 
352  const float_v C31_in = fC31;
353 
354  fC31(mask) += dz * fC33;
355  fC11(mask) += dz * ( fC31 + C31_in );
356  fC30(mask) += dzC32_in;
357 
358  fC40(mask) += dz * fC42;
359  fC41(mask) += dz * fC43;
360 
361  return mask;
362 }
363 
364 inline float_m PndFTSCATrackParamVector::Filter( const float_v &x0, const float_v &y0, const float_v &r, const FTSCAStripInfoVector &info, float_v err2, const float_m &active, const float_v& /*chi2Cut*/ ) //[R.K. 9/2018] unused
365  // Adds the tube measurement with the Kalman filter
366  // @beta is angle between the strip and z-axis, clockwise. The wire equation is: {x,y,z} - {x0,y0,z0} = t*e_s, where ort e_s = { 0, -sinB, cos B }
367 {
368  // linearize track in current point, which must be x == x0
369  // distance between wire and track are : r = - ({x,y,z} - {x0,y0,z0}) * e_o / |e_o|, where ort e_o = |[e_t x e_s]|
370 
371  float_m success = active;
372 
373  success &= ( err2 > 1.e-8f );
374 
375  //const float_v& h0 = info.cos;
376  //const float_v& h1 = info.sin;
377  float_v h0(Vc::Zero);
378  float_v h1(Vc::Zero);
379  h0(active) = info.cos;
380  h1(active) = info.sin;
381  /*foreach_bit(unsigned short iV, active) {
382  h0[iV] = info.cos;
383  h1[iV] = info.sin;
384  }*/
385 
386 
387  const float_v tx = h0*fTx + h1*fTy;
388 
389  float_v zeta;
390 
391 // float_v Sign_1 = float_v(0.0);
392 
393  float_v rCorrection = sqrt(1.f+tx*tx);
394 
395 
396  float_v Diff = h0*(fX - x0) + h1*(fY - y0);
397  float_v zeta_1, zeta_2;
398  zeta_1 = h0*(fX - x0) + h1*(fY - y0) - r*rCorrection;
399  zeta_2 = h0*(fX - x0) + h1*(fY - y0) + r*rCorrection;
400 
401  for (int iDiff=0; iDiff < 4; iDiff++)
402  {
403  float Diff_f = (float)Diff[iDiff];
404  if (Diff_f> 0.f )
405  {
406  //zeta = h0*(fX - x0) + h1*(fY - y0) - r*rCorrection;
407  zeta[iDiff] = zeta_1[iDiff];
408  }
409  else
410  {
411  zeta[iDiff] = zeta_2[iDiff];
412  }
413 
414  } // iDiff
415 
416  //const float_v zeta2 = h0*(fX - x0) + h1*(fY - y0) + r*rCorrection;
417 
418  if (err2[0] < 0.01f) err2[0] = err2[0] + r[0]*r[0];
419  err2 *= rCorrection*rCorrection;
420 
421  float_v wi, zetawi, HCH;
422  float_v F0, F1, F2, F3, F4;
423  float_v K1, K2, K3, K4;
424 
425  if (fC00[0] < 0.f) fC00[0]=0.f;
426 
427  // F = CH'
428  F0 = h0*fC00 + h1*fC10;
429  F1 = h0*fC10 + h1*fC11;
430 
431  //HCH = ( F0*h0 + F1*h1 );
432  HCH = (h0*h0*fC00 + 2.f*h0*h1*fC10) + h1*h1*fC11;
433 
434  F2 = h0*fC20 + h1*fC21;
435  F3 = h0*fC30 + h1*fC31;
436  F4 = h0*fC40 + h1*fC41;
437 
438  float_v dChi2(Vc::Zero);
439 #if 0 // use mask
440  cout<<"we don't have to be here for now \n";
441  const float_m mask = success && (HCH > err2 * 16.f);
442  HCH(mask) += err2;
443  wi(mask) = 1.f/HCH ;
444  zetawi(mask) = zeta *wi;
445  dChi2(mask) += (zeta * zetawi);
446  fChi2 += dChi2;
447 #else
448  wi = 1.f/( err2 + HCH );
449  zetawi = zeta *wi;
450  dChi2(success) += zeta * zetawi;
451  //success &= dChi2 < chi2Cut;
452  fChi2 += dChi2;
453  //fChi2(success) += zeta * zetawi;
454 #endif // 0
455 
456  //success &= fChi2 < chi2Cut;
457 
458  K1 = F1*wi;
459  K2 = F2*wi;
460  K3 = F3*wi;
461  K4 = F4*wi;
462 
463  fX(success) -= F0*zetawi;
464  fY(success) -= F1*zetawi;
465  fTx(success) -= F2*zetawi;
466  fTy(success) -= F3*zetawi;
467  fQP(success) -= F4*zetawi;
468 
469  fC00(success) -= F0*F0*wi;
470  fC10(success) -= K1*F0;
471  fC11(success) -= K1*F1;
472  fC20(success) -= K2*F0;
473  fC21(success) -= K2*F1;
474  fC22(success) -= K2*F2;
475  fC30(success) -= K3*F0;
476  fC31(success) -= K3*F1;
477  fC32(success) -= K3*F2;
478  fC33(success) -= K3*F3;
479  fC40(success) -= K4*F0;
480  fC41(success) -= K4*F1;
481  fC42(success) -= K4*F2;
482  fC43(success) -= K4*F3;
483  fC44(success) -= K4*F4;
484 
485  fNDF( static_cast<int_m>(success) ) += 1;
486 
487  return success;
488 }
489 
490 inline float_m PndFTSCATrackParamVector::TransportToX0WithMaterial( const float_v &x0, const L1FieldRegion &F, const L1MaterialInfo &material, float_v &qp0, const float_m &mask )
491 { // TODO material
492  float_m active = mask;
493  active &= TransportToX0( x0, F, qp0, active );
494  //active &= TransportToX0Line(x0, active);
495  //active &= RK4TransportToX0( x0, F, qp0, active );
496  active &= PassMaterial( material, qp0, active );
497  const float_v mass2 = 0.13957f*0.13957f;
498  float_v direction = -1.f;
499  if(fDirection)
500  direction = 1.f;
501  //std::cout << "qp " << qp0[0] << " " << " fQP " << fQP[0]<<" ";
502  EnergyLossCorrection(mass2, material.RadThick, qp0, direction, active);
503  //std::cout << qp0[0] << " " << fQP[0] << std::endl;
504  return active;
505 }
506 
507 inline float_m PndFTSCATrackParamVector::RK4TransportToX0( const float_v &x0_out, const L1FieldRegion &/*F*/, const float_v &/*qp0*/, const float_m &mask ) //[R.K. 9/2018] 2 unused
508 {
509  // Forth-order Runge-Kutta method for solution of the equation
510  // of motion of a particle with parameter qp = Q /P
511  // in the magnetic field B()
512  //
513  // | x | tx
514  // | y | ty
515  // d/dz = | tx| = ft * ( ty * ( B(3)+tx*b(1) ) - (1+tx**2)*B(2) )
516  // | ty| ft * (-tx * ( B(3)+ty+b(2) - (1+ty**2)*B(1) ) ,
517  //
518  // where ft = c_light*qp*sqrt ( 1 + tx**2 + ty**2 ) .
519  //
520  // In the following for RK step :
521  //
522  // x=x[0], y=x[1], tx=x[3], ty=x[4].
523  //
524  //========================================================================
525 
526 // const float ZERO = 0.0, ONE = 1.;
527 
528  //RK4ORDER
529  float_v xOut[5];
530  float_v Fmat[25];
531  const float_v fC = 0.000299792458f;
532 
533  float_v coef[4];
534  coef[0] = 0.f; coef[1] = 0.5f;
535  coef[2] = 0.5f; coef[3] = 1.f;
536 
537  float_v xIn[4];
538  xIn[0] = fX; xIn[2] = fTx;
539  xIn[1] = fY; xIn[3] = fTy; xIn[4] = fQP;
540 
541  float_v Ax[4], Ay[4];
542  float_v dAx_dtx[4], dAy_dtx[4], dAx_dty[4], dAy_dty[4];
543  float_v k[4][4];
544 
545  float_v h = x0_out - fZ;
546  //cout<<"Step h: "<<h<<endl;
547  float_v hC = h * fC;
548  float_v hCqp = h * fC * xIn[4];
549  float_v x0[4];
550 
551  float_v x[4];
552  x[0] = fX; x[2] = fTx;
553  x[1] = fY; x[3] = fTy; //x[4] = fQP;
554 
555  for (unsigned int iStep = 0; iStep < 4; iStep++) { // 1
556  if (iStep > 0) {
557 
558  x[0] = xIn[0] + coef[iStep] * k[0][iStep - 1];
559  x[1] = xIn[1] + coef[iStep] * k[1][iStep - 1];
560  x[2] = xIn[2] + coef[iStep] * k[2][iStep - 1];
561  x[3] = xIn[3] + coef[iStep] * k[3][iStep - 1];
562  }
563 
564  float_v bx(0.f), by(0.f), bz(0.f);
565  double p[3];
566  double b[3];
567  for (int xxx = 0; xxx < float_v::Size; xxx++)
568  {
569 // if (fabs(x[0][xxx])>0.001)
570  {
571  //double p[3] = {fX[ctr], fY[ctr], fZ[ctr]};
572  //double b[3] = {0., 0., 0.};
573  p[0] = x[0][xxx];
574  p[1] = x[1][xxx];
575  p[2] = fZ[xxx] + coef[iStep][xxx] * h[xxx];
576  b[0] = 0;
577  b[1] = 0;
578  b[2] = 0;
579  FairRunAna::Instance()->GetField()->Field(p,b);//CbmKF::Instance()->GetMagneticField()->GetFieldValue(p,b);
580  bx[xxx] = b[0];
581  by[xxx] = b[1];
582  bz[xxx] = b[2];
583  }
584  }
585  //fField->GetFieldValue(x[0], x[1], zIn + coef[iStep] * h, Bx, By, Bz);
586  float_v tx = x[2];
587  float_v ty = x[3];
588  float_v txtx = tx * tx;
589  float_v tyty = ty * ty;
590  float_v txty = tx * ty;
591  float_v txtxtyty1= 1.f + txtx + tyty;
592  float_v t1 = sqrt(txtxtyty1);
593  float_v t2 = 1.f / txtxtyty1;
594 
595  Ax[iStep] = ( txty * bx + ty * bz - ( 1.f + txtx ) * by ) * t1;
596  Ay[iStep] = (-txty * by - tx * bz + ( 1.f + tyty ) * bx ) * t1;
597 
598  dAx_dtx[iStep] = Ax[iStep] * tx * t2 + ( ty * bx - 2.f * tx * by ) * t1;
599  dAx_dty[iStep] = Ax[iStep] * ty * t2 + ( tx * bx + bz ) * t1;
600  dAy_dtx[iStep] = Ay[iStep] * tx * t2 + (-ty * by - bz ) * t1;
601  dAy_dty[iStep] = Ay[iStep] * ty * t2 + (-tx * by + 2.f * ty * bx ) * t1;
602 
603  k[0][iStep] = tx * h;
604  k[1][iStep] = ty * h;
605  k[2][iStep] = Ax[iStep] * hCqp;
606  k[3][iStep] = Ay[iStep] * hCqp;
607 
608  } // 1
609 
610  //in + k[0]/6.f + k[1]/3.f + k[2]/3.f + k[3]/6.f;
611  //for (unsigned int i = 0; i < 4; i++) { xOut[i] = CalcOut(xIn[i], k[i]); }
612  xOut[0] = xIn[0] + k[0][0]/6.f + k[0][1]/3.f + k[0][2]/3.f + k[0][3]/6.f;
613  xOut[1] = xIn[1] + k[1][0]/6.f + k[1][1]/3.f + k[1][2]/3.f + k[1][3]/6.f;
614  xOut[2] = xIn[2] + k[2][0]/6.f + k[2][1]/3.f + k[2][2]/3.f + k[2][3]/6.f;
615  xOut[3] = xIn[3] + k[3][0]/6.f + k[3][1]/3.f + k[3][2]/3.f + k[3][3]/6.f;
616  xOut[4] = xIn[4];
617 
618  fX(mask) = xOut[0];
619  fY(mask) = xOut[1];
620  fTx(mask) = xOut[2];
621  fTy(mask) = xOut[3];
622  fQP(mask) = xOut[4];
623  fZ(mask) += h;
624  // Calculation of the derivatives
625 
626  // derivatives dx/dx and dx/dy
627  // dx/dx
628  Fmat[0] = 1.f;
629  Fmat[5] = 0.f;
630  Fmat[10] = 0.f;
631  Fmat[15] = 0.f;
632  Fmat[20] = 0.f;
633  // dx/dy
634  Fmat[1] = 0.f;
635  Fmat[6] = 1.f;
636  Fmat[11] = 0.f;
637  Fmat[16] = 0.f;
638  Fmat[21] = 0.f;
639  // end of derivatives dx/dx and dx/dy
640 
641  // Derivatives dx/tx
642  x[0] = x0[0] = 0.f;
643  x[1] = x0[1] = 0.f;
644  x[2] = x0[2] = 1.f;
645  x[3] = x0[3] = 0.f;
646  for (unsigned int iStep = 0; iStep < 4; iStep++) { // 2
647  if (iStep > 0) {
648  x[0] = xIn[0] + coef[iStep] * k[0][iStep - 1];
649  x[1] = xIn[1] + coef[iStep] * k[1][iStep - 1];
650  x[3] = xIn[3] + coef[iStep] * k[3][iStep - 1];
651  }
652 
653  k[0][iStep] = x[2] * h;
654  k[1][iStep] = x[3] * h;
655  //k[2][iStep] = (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]) * hCqp;
656  k[3][iStep] = (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]) * hCqp;
657  } // 2
658 
659 
660  Fmat[2] = x0[0] + k[0][0]/6.f + k[0][1]/3.f + k[0][2]/3.f + k[0][3]/6.f;
661 
662  Fmat[7] = x0[1] + k[1][0]/6.f + k[1][1]/3.f + k[1][2]/3.f + k[1][3]/6.f;
663  Fmat[12] = 1.f;
664 
665  Fmat[17] = x0[3] + k[3][0]/6.f + k[3][1]/3.f + k[3][2]/3.f + k[3][3]/6.f;
666  Fmat[22] = 0.f;
667  // end of derivatives dx/dtx
668 
669  // Derivatives dx/ty
670  x[0] = x0[0] = 0.f;
671  x[1] = x0[1] = 0.f;
672  x[2] = x0[2] = 0.f;
673  x[3] = x0[3] = 1.f;
674  for (unsigned int iStep = 0; iStep < 4; iStep++) { // 4
675  if (iStep > 0) {
676  x[0] = xIn[0] + coef[iStep] * k[0][iStep - 1];
677  x[1] = xIn[1] + coef[iStep] * k[1][iStep - 1];
678  x[2] = xIn[2] + coef[iStep] * k[2][iStep - 1];
679  }
680 
681  k[0][iStep] = x[2] * h;
682  k[1][iStep] = x[3] * h;
683  k[2][iStep] = (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]) * hCqp;
684  //k[3][iStep] = (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]) * hCqp;
685  } // 4
686 
687  Fmat[3] = x0[0] + k[0][0]/6.f + k[0][1]/3.f + k[0][2]/3.f + k[0][3]/6.f;
688 
689  Fmat[8] = x0[1] + k[1][0]/6.f + k[1][1]/3.f + k[1][2]/3.f + k[1][3]/6.f;
690 
691  Fmat[13] = x0[2] + k[2][0]/6.f + k[2][1]/3.f + k[2][2]/3.f + k[2][3]/6.f;
692  Fmat[18] = 1.f;
693  Fmat[23] = 0.f;
694  // end of derivatives dx/dty
695 
696  // Derivatives dx/dqp
697  x[0] = x0[0] = 0.f;
698  x[1] = x0[1] = 0.f;
699  x[2] = x0[2] = 0.f;
700  x[3] = x0[3] = 0.f;
701  for (unsigned int iStep = 0; iStep < 4; iStep++) { // 4
702  if (iStep > 0) {
703  x[0] = xIn[0] + coef[iStep] * k[0][iStep - 1];
704  x[1] = xIn[1] + coef[iStep] * k[1][iStep - 1];
705  x[2] = xIn[2] + coef[iStep] * k[2][iStep - 1];
706  x[3] = xIn[3] + coef[iStep] * k[3][iStep - 1];
707  }
708 
709  k[0][iStep] = x[2] * h;
710  k[1][iStep] = x[3] * h;
711  k[2][iStep] = Ax[iStep] * hC +
712  hCqp * (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]);
713  k[3][iStep] = Ay[iStep] * hC +
714  hCqp * (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]);
715  } // 4
716 
717 
718  Fmat[4] = x0[0] + k[0][0]/6.f + k[0][1]/3.f + k[0][2]/3.f + k[0][3]/6.f;
719 
720  Fmat[9] = x0[1] + k[1][0]/6.f + k[1][1]/3.f + k[1][2]/3.f + k[1][3]/6.f;
721 
722  Fmat[14] = x0[2] + k[2][0]/6.f + k[2][1]/3.f + k[2][2]/3.f + k[2][3]/6.f;
723 
724  Fmat[19] = x[3] + k[3][0]/6.f + k[3][1]/3.f + k[3][2]/3.f + k[3][3]/6.f;
725  Fmat[24] = 1.f;
726  // end of derivatives dx/dqp
727 
728 
729  // end calculation of the derivatives
730  // F*C*Ft
731  float_v A = fC20 + Fmat[2] * fC22 + Fmat[3] * fC32 + Fmat[4] * fC42;
732  float_v B = fC30 + Fmat[2] * fC32 + Fmat[3] * fC33 + Fmat[4] * fC43;
733  float_v C = fC40 + Fmat[2] * fC42 + Fmat[3] * fC43 + Fmat[4] * fC44;
734 
735  float_v D = fC21 + Fmat[7] * fC22 + Fmat[8] * fC32 + Fmat[9] * fC42;
736  float_v E = fC31 + Fmat[7] * fC32 + Fmat[8] * fC33 + Fmat[9] * fC43;
737  float_v G = fC41 + Fmat[7] * fC42 + Fmat[8] * fC43 + Fmat[9] * fC44;
738 
739  float_v H = fC22 + Fmat[13] * fC32 + Fmat[14] * fC42;
740  float_v I = fC32 + Fmat[13] * fC33 + Fmat[14] * fC43;
741  float_v J = fC42 + Fmat[13] * fC43 + Fmat[14] * fC44;
742 
743  float_v K = fC43 + Fmat[17] * fC42 + Fmat[19] * fC44;
744 
745  fC00(mask) = fC00 + Fmat[2] * fC20 + Fmat[3] * fC30 + Fmat[4] * fC40 + A * Fmat[2] + B * Fmat[3] + C * Fmat[4];
746  fC10(mask) = fC10 + Fmat[2] * fC21 + Fmat[3] * fC31 + Fmat[4] * fC41 + A * Fmat[7] + B * Fmat[8] + C * Fmat[9];
747  fC20(mask) = A + B * Fmat[13] + C * Fmat[14];
748  fC30(mask) = B + A * Fmat[17] + C * Fmat[19];
749  fC40(mask) = C;
750 
751  fC11(mask) = fC11 + Fmat[7] * fC21 + Fmat[8] * fC31 + Fmat[9] * fC41 + D * Fmat[7] + E * Fmat[8] + G * Fmat[9];
752  fC21(mask) = D + E * Fmat[13] + G * Fmat[14];
753  fC31(mask) = E + D * Fmat[17] + G * Fmat[19];
754  fC41(mask) = G;
755 
756  fC22(mask) = H + I * Fmat[13] + J * Fmat[14];
757  fC32(mask) = I + H * Fmat[17] + J * Fmat[19];
758  fC42(mask) = J;
759 
760  fC33(mask) = fC33 + Fmat[17] * fC32 + Fmat[19] * fC43 + (Fmat[17] * fC22 + fC32 + Fmat[19] * fC42) * Fmat[17] + K * Fmat[19];
761  fC43(mask) = K;
762 
763  fC44(mask) = fC44;
764 
765  return mask;
766 }
767 
768 inline float_m PndFTSCATrackParamVector::TransportToX0( const float_v &x0_out, const L1FieldRegion &F, const float_v &qp0, const float_m &mask )
769 {
770 
771  //cout<<"Extrapolation..."<<endl;
772  //
773  // Part of the analytic extrapolation formula with error (c_light*B*dz)^4/4!
774  //
775  /*cout<<"before extrp fX fY fTx fTy fQP "<<fX[0]<<" "<<fY[0]<<" "<<fTx[0]<<" "<<fTy[0]<<" "<<fQP[0] << " z " << fZ[0] <<endl;
776  cout<<"transport to x0 "<<x0_out[0]<<endl;*/
777  const float_v c_light = 0.000299792458f;
778 // std::cout <<"Extrapolate 0 " <<std::endl;
779 // std::cout << "fX " << fX << " fY " << fY << " fZ " << fZ << std::endl;
780 // std::cout << "fTx " << fTx << " fTy " << fTy << " Qp " << fQP << std::endl;
781 
782  const float_v
783  c1 = 1.f, c2 = 2.f, c3 = 3.f, c4 = 4.f, c6 = 6.f, c9 = 9.f, c15 = 15.f, c18 = 18.f, c45 = 45.f,
784  c2i = 1.f/2.f, c3i = 1.f/3.f, c6i = 1.f/6.f, c12i = 1.f/12.f;
785 
786  float_v dz(Vc::Zero);
787  dz(mask) = (x0_out - fZ);
788  //const float_v dz = (x0_out - fZ);
789  const float_v dz2 = dz*dz;
790  const float_v dz3 = dz2*dz;
791  // construct coefficients
792 
793  const float_v x = fTx;
794  const float_v y = fTy;
795  const float_v xx = x*x;
796  const float_v xy = x*y;
797  const float_v yy = y*y;
798  const float_v y2 = y*c2;
799  const float_v x2 = x*c2;
800  const float_v x4 = x*c4;
801  const float_v xx31 = xx*c3+c1;
802  const float_v xx159 = xx*c15+c9;
803 
804  const float_v Ay = -xx-c1;
805  const float_v Ayy = x*(xx*c3+c3);
806  const float_v Ayz = -c2*xy;
807  const float_v Ayyy = -(c15*xx*xx+c18*xx+c3);
808 
809  const float_v Ayy_fX = c3*xx31;
810  const float_v Ayyy_fX = -x4*xx159;
811 
812  const float_v bx = yy+c1;
813  const float_v Byy = y*xx31;
814  const float_v Byz = c2*xx+c1;
815  const float_v Byyy = -xy*xx159;
816 
817  const float_v Byy_fX = c6*xy;
818  const float_v Byyy_fX = -y*(c45*xx+c9);
819  const float_v Byyy_fY = -x*xx159;
820 
821  // end of coefficients calculation
822 
823  const float_v t2 = c1 + xx + yy;
824  const float_v t = sqrt( t2 );
825  const float_v h = qp0*c_light;
826  const float_v ht = h*t;
827 
828 
829  // get field integrals
830  const float_v ddz = fZ-F.z0;
831  float_v Fx0 = F.cx0 + F.cx1*ddz + F.cx2*ddz*ddz;
832  float_v Fx1 = (F.cx1 + c2*F.cx2*ddz)*dz;
833  float_v Fx2 = F.cx2*dz2;
834  float_v Fy0 = F.cy0 + F.cy1*ddz + F.cy2*ddz*ddz;
835  float_v Fy1 = (F.cy1 + c2*F.cy2*ddz)*dz;
836  float_v Fy2 = F.cy2*dz2;
837  float_v Fz0 = F.cz0 + F.cz1*ddz + F.cz2*ddz*ddz;
838  float_v Fz1 = (F.cz1 + c2*F.cz2*ddz)*dz;
839  float_v Fz2 = F.cz2*dz2;
840 
841  const float_v sx = ( Fx0 + Fx1*c2i + Fx2*c3i );
842  const float_v sy = ( Fy0 + Fy1*c2i + Fy2*c3i );
843  const float_v sz = ( Fz0 + Fz1*c2i + Fz2*c3i );
844 
845  const float_v Sx = ( Fx0*c2i + Fx1*c6i + Fx2*c12i );
846  const float_v Sy = ( Fy0*c2i + Fy1*c6i + Fy2*c12i );
847  const float_v Sz = ( Fz0*c2i + Fz1*c6i + Fz2*c12i );
848 
849  float_v syz;
850  {
851  const float_v
852  d = 1.f/360.f,
853  c00 = 30.f*6.f*d, c01 = 30.f*2.f*d, c02 = 30.f*d,
854  c10 = 3.f*40.f*d, c11 = 3.f*15.f*d, c12 = 3.f*8.f*d,
855  c20 = 2.f*45.f*d, c21 = 2.f*2.f*9.f*d, c22 = 2.f*2.f*5.f*d;
856  syz = Fy0*( c00*Fz0 + c01*Fz1 + c02*Fz2)
857  + Fy1*( c10*Fz0 + c11*Fz1 + c12*Fz2)
858  + Fy2*( c20*Fz0 + c21*Fz1 + c22*Fz2) ;
859  }
860 
861  float_v Syz;
862  {
863  const float_v
864  d = 1.f/2520.f,
865  c00 = 21.f*20.f*d, c01 = 21.f*5.f*d, c02 = 21.f*2.f*d,
866  c10 = 7.f*30.f*d, c11 = 7.f*9.f*d, c12 = 7.f*4.f*d,
867  c20 = 2.f*63.f*d, c21 = 2.f*21.f*d, c22 = 2.f*10.f*d;
868  Syz = Fy0*( c00*Fz0 + c01*Fz1 + c02*Fz2 )
869  + Fy1*( c10*Fz0 + c11*Fz1 + c12*Fz2 )
870  + Fy2*( c20*Fz0 + c21*Fz1 + c22*Fz2 ) ;
871  }
872 
873  const float_v syy = sy*sy*c2i;
874  const float_v syyy = syy*sy*c3i;
875 
876  float_v Syy ;
877  {
878  const float_v
879  d= 1.f/2520.f, c00= 420.f*d, c01= 21.f*15.f*d, c02= 21.f*8.f*d,
880  c03= 63.f*d, c04= 70.f*d, c05= 20.f*d;
881  Syy = Fy0*(c00*Fy0+c01*Fy1+c02*Fy2) + Fy1*(c03*Fy1+c04*Fy2) + c05*Fy2*Fy2 ;
882  }
883 
884  float_v Syyy;
885  {
886  const float_v
887  d = 1.f/181440.f,
888  c000 = 7560*d, c001 = 9*1008*d, c002 = 5*1008*d,
889  c011 = 21*180*d, c012 = 24*180*d, c022 = 7*180*d,
890  c111 = 540*d, c112 = 945*d, c122 = 560*d, c222 = 112*d;
891  const float_v Fy22 = Fy2*Fy2;
892  Syyy = Fy0*( Fy0*(c000*Fy0+c001*Fy1+c002*Fy2)+ Fy1*(c011*Fy1+c012*Fy2)+c022*Fy22 )
893  + Fy1*( Fy1*(c111*Fy1+c112*Fy2)+c122*Fy22) + c222*Fy22*Fy2 ;
894  }
895 
896 
897  const float_v sA1 = sx*xy + sy*Ay + sz*y ;
898  const float_v sA1_fX = sx*y - sy*x2 ;
899  const float_v sA1_fY = sx*x + sz ;
900 
901  const float_v sB1 = sx*bx - sy*xy - sz*x ;
902  const float_v sB1_fX = -sy*y - sz ;
903  const float_v sB1_fY = sx*y2 - sy*x ;
904 
905  const float_v SA1 = Sx*xy + Sy*Ay + Sz*y ;
906  const float_v SA1_fX = Sx*y - Sy*x2 ;
907  const float_v SA1_fY = Sx*x + Sz;
908 
909  const float_v SB1 = Sx*bx - Sy*xy - Sz*x ;
910  const float_v SB1_fX = -Sy*y - Sz;
911  const float_v SB1_fY = Sx*y2 - Sy*x;
912 
913 
914  const float_v sA2 = syy*Ayy + syz*Ayz ;
915  const float_v sA2_fX = syy*Ayy_fX - syz*y2 ;
916  const float_v sA2_fY = -syz*x2 ;
917  const float_v sB2 = syy*Byy + syz*Byz ;
918  const float_v sB2_fX = syy*Byy_fX + syz*x4 ;
919  const float_v sB2_fY = syy*xx31 ;
920 
921  const float_v SA2 = Syy*Ayy + Syz*Ayz ;
922  const float_v SA2_fX = Syy*Ayy_fX - Syz*y2 ;
923  const float_v SA2_fY = -Syz*x2 ;
924  const float_v SB2 = Syy*Byy + Syz*Byz ;
925  const float_v SB2_fX = Syy*Byy_fX + Syz*x4 ;
926  const float_v SB2_fY = Syy*xx31 ;
927 
928  const float_v sA3 = syyy*Ayyy ;
929  const float_v sA3_fX = syyy*Ayyy_fX;
930  const float_v sB3 = syyy*Byyy ;
931  const float_v sB3_fX = syyy*Byyy_fX;
932  const float_v sB3_fY = syyy*Byyy_fY;
933 
934 
935  const float_v SA3 = Syyy*Ayyy ;
936  const float_v SA3_fX = Syyy*Ayyy_fX;
937  const float_v SB3 = Syyy*Byyy ;
938  const float_v SB3_fX = Syyy*Byyy_fX;
939  const float_v SB3_fY = Syyy*Byyy_fY;
940 
941  const float_v ht1 = ht*dz;
942  const float_v ht2 = ht*ht*dz2;
943  const float_v ht3 = ht*ht*ht*dz3;
944  const float_v ht1sA1 = ht1*sA1;
945  const float_v ht1sB1 = ht1*sB1;
946  const float_v ht1SA1 = ht1*SA1;
947  const float_v ht1SB1 = ht1*SB1;
948  const float_v ht2sA2 = ht2*sA2;
949  const float_v ht2SA2 = ht2*SA2;
950  const float_v ht2sB2 = ht2*sB2;
951  const float_v ht2SB2 = ht2*SB2;
952  const float_v ht3sA3 = ht3*sA3;
953  const float_v ht3sB3 = ht3*sB3;
954  const float_v ht3SA3 = ht3*SA3;
955  const float_v ht3SB3 = ht3*SB3;
956 
957  fX (mask) += ((x + ht1SA1 + ht2SA2 + ht3SA3)*dz);
958  fY (mask) += ((y + ht1SB1 + ht2SB2 + ht3SB3)*dz);
959  fTx(mask) += (ht1sA1 + ht2sA2 + ht3sA3);
960  fTy(mask) += (ht1sB1 + ht2sB2 + ht3sB3);
961  fZ (mask) += (dz);
962 
963 // std::cout <<"Extrapolate 1 " <<std::endl;
964 // std::cout << "fX " << fX << " fY " << fY << " fZ " << fZ << std::endl;
965 // std::cout << "fTx " << fTx << " fTy " << fTy << " Qp " << fQP << std::endl;
966 
967  const float_v ctdz = c_light*t*dz;
968  const float_v ctdz2 = c_light*t*dz2;
969 
970  const float_v dqp = fQP - qp0;
971  const float_v t2i = c1*rcp(t2);// /t2;
972  const float_v xt2i = x*t2i;
973  const float_v yt2i = y*t2i;
974  const float_v tmp0 = ht1SA1 + c2*ht2SA2 + c3*ht3SA3;
975  const float_v tmp1 = ht1SB1 + c2*ht2SB2 + c3*ht3SB3;
976  const float_v tmp2 = ht1sA1 + c2*ht2sA2 + c3*ht3sA3;
977  const float_v tmp3 = ht1sB1 + c2*ht2sB2 + c3*ht3sB3;
978 
979  const float_v j02 = dz*(c1 + xt2i*tmp0 + ht1*SA1_fX + ht2*SA2_fX + ht3*SA3_fX);
980  const float_v j12 = dz*( xt2i*tmp1 + ht1*SB1_fX + ht2*SB2_fX + ht3*SB3_fX);
981  const float_v j22 = c1 + xt2i*tmp2 + ht1*sA1_fX + ht2*sA2_fX + ht3*sA3_fX ;
982  const float_v j32 = xt2i*tmp3 + ht1*sB1_fX + ht2*sB2_fX + ht3*sB3_fX ;
983 
984  const float_v j03 = dz*( yt2i*tmp0 + ht1*SA1_fY + ht2*SA2_fY );
985  const float_v j13 = dz*(c1 + yt2i*tmp1 + ht1*SB1_fY + ht2*SB2_fY + ht3*SB3_fY );
986  const float_v j23 = yt2i*tmp2 + ht1*sA1_fY + ht2*sA2_fY ;
987  const float_v j33 = c1 + yt2i*tmp3 + ht1*sB1_fY + ht2*sB2_fY + ht3*sB3_fY ;
988 
989  const float_v j04 = ctdz2*( SA1 + c2*ht1*SA2 + c3*ht2*SA3 );
990  const float_v j14 = ctdz2*( SB1 + c2*ht1*SB2 + c3*ht2*SB3 );
991  const float_v j24 = ctdz *( sA1 + c2*ht1*sA2 + c3*ht2*sA3 );
992  const float_v j34 = ctdz *( sB1 + c2*ht1*sB2 + c3*ht2*sB3 );
993 
994 
995  // extrapolate inverse momentum
996  fX (mask)+=j04*dqp;
997  fY (mask)+=j14*dqp;
998  fTx(mask)+=j24*dqp;
999  fTy(mask)+=j34*dqp;
1000 
1001 // std::cout <<"Extrapolate 2 " <<std::endl;
1002 // std::cout << "fX " << fX << " fY " << fY << " fZ " << fZ << std::endl;
1003 // std::cout << "fTx " << fTx << " fTy " << fTy << " Qp " << fQP << std::endl;
1004  // covariance matrix transport
1005 
1006  const float_v c42 = fC42, c43 = fC43;
1007 
1008  const float_v cj00 = fC00 + fC20*j02 + fC30*j03 + fC40*j04;
1009 // const float_v cj10 = fC10 + fC21*j02 + fC31*j03 + fC41*j04;
1010  const float_v cj20 = fC20 + fC22*j02 + fC32*j03 + c42*j04;
1011  const float_v cj30 = fC30 + fC32*j02 + fC33*j03 + c43*j04;
1012 
1013  const float_v cj01 = fC10 + fC20*j12 + fC30*j13 + fC40*j14;
1014  const float_v cj11 = fC11 + fC21*j12 + fC31*j13 + fC41*j14;
1015  const float_v cj21 = fC21 + fC22*j12 + fC32*j13 + c42*j14;
1016  const float_v cj31 = fC31 + fC32*j12 + fC33*j13 + c43*j14;
1017 
1018 // const float_v cj02 = fC20*j22 + fC30*j23 + fC40*j24;
1019 // const float_v cj12 = fC21*j22 + fC31*j23 + fC41*j24;
1020  const float_v cj22 = fC22*j22 + fC32*j23 + c42*j24;
1021  const float_v cj32 = fC32*j22 + fC33*j23 + c43*j24;
1022 
1023 // const float_v cj03 = fC20*j32 + fC30*j33 + fC40*j34;
1024 // const float_v cj13 = fC21*j32 + fC31*j33 + fC41*j34;
1025  const float_v cj23 = fC22*j32 + fC32*j33 + c42*j34;
1026  const float_v cj33 = fC32*j32 + fC33*j33 + c43*j34;
1027 
1028  fC40(mask)+= (c42*j02 + c43*j03 + fC44*j04); // cj40
1029  fC41(mask)+= (c42*j12 + c43*j13 + fC44*j14); // cj41
1030  fC42(mask) = (c42*j22 + c43*j23 + fC44*j24); // cj42
1031  fC43(mask) = (c42*j32 + c43*j33 + fC44*j34); // cj43
1032 
1033  fC00(mask) = (cj00 + j02*cj20 + j03*cj30 + j04*fC40);
1034  fC10(mask) = (cj01 + j02*cj21 + j03*cj31 + j04*fC41);
1035  fC11(mask) = (cj11 + j12*cj21 + j13*cj31 + j14*fC41);
1036 
1037  //cout<<"transport \n";
1038  //cout<<" C00 " << fC00 << " C10 " << fC10 << " C11 " << fC11 << endl;
1039  fC20(mask) = (j22*cj20 + j23*cj30 + j24*fC40);
1040  fC30(mask) = (j32*cj20 + j33*cj30 + j34*fC40);
1041  fC21(mask) = (j22*cj21 + j23*cj31 + j24*fC41);
1042  fC31(mask) = (j32*cj21 + j33*cj31 + j34*fC41);
1043  fC22(mask) = (j22*cj22 + j23*cj32 + j24*fC42);
1044  fC32(mask) = (j32*cj22 + j33*cj32 + j34*fC42);
1045  fC33(mask) = (j32*cj23 + j33*cj33 + j34*fC43);
1046 
1047  return mask;
1048 }
1049 
1050 inline float_m PndFTSCATrackParamVector::PassMaterial( const L1MaterialInfo &info, const float_v &qp0, const float_m &mask )
1051 {
1052  const float_v mass2 = 0.1395679f*0.1395679f;
1053 
1054  float_v txtx = fTx*fTx;
1055  float_v tyty = fTy*fTy;
1056  float_v txtx1 = txtx + 1.f;
1057  float_v h = txtx + tyty;
1058  float_v t = sqrt(txtx1 + tyty);
1059  float_v h2 = h*h;
1060  float_v qp0t = qp0*t;
1061 
1062  const float_v c1=0.0136f, c2=c1*0.038f, c3=c2*0.5f, c4=-c3/2.0f, c5=c3/3.0f, c6=-c3/4.0f;
1063 
1064  float_v s0 = (c1+c2*info.logRadThick + c3*h + h2*(c4 + c5*h +c6*h2) )*qp0t;
1065  float_v a = ( (t+mass2*qp0*qp0t)*info.RadThick*s0*s0 );
1066  //std::cout<<"info.RadThick info.logRadThick "<<info.RadThick<<" "<<info.logRadThick<<std::endl;
1067 // std::cout <<" a " << a << std::endl;
1068 // a=0.000005;
1069  fC22(mask) += txtx1*a;
1070  fC32(mask) += fTx*fTy*a; fC33(mask) += (1.f + tyty)*a;
1071 
1072  return mask;
1073 }
1074 
1075 inline float_v PndFTSCATrackParamVector::ApproximateBetheBloch( const float_v &bg2 )
1076 {
1077  //
1078  // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
1079  //
1080  // bg2 - (beta*gamma)^2
1081  // kp0 - density [g/cm^3]
1082  // kp1 - density effect first junction point
1083  // kp2 - density effect second junction point
1084  // kp3 - mean excitation energy [GeV]
1085  // kp4 - mean Z/A
1086  //
1087  // The default values for the kp* parameters are for silicon.
1088  // The returned value is in [GeV/(g/cm^2)].
1089  //
1090 
1091  const float_v &kp0 = 2.33f;
1092  const float_v &kp1 = 0.20f;
1093  const float_v &kp2 = 3.00f;
1094  const float_v &kp3 = 173e-9f;
1095  const float_v &kp4 = 0.49848f;
1096 
1097  const float mK = 0.307075e-3f; // [GeV*cm^2/g]
1098  const float _2me = 1.022e-3f; // [GeV/c^2]
1099  const float_v &rho = kp0;
1100  const float_v &x0 = kp1 * 2.303f;
1101  const float_v &x1 = kp2 * 2.303f;
1102  const float_v &mI = kp3;
1103  const float_v &mZA = kp4;
1104  const float_v &maxT = _2me * bg2; // neglecting the electron mass
1105 
1106  //*** Density effect
1107  float_v d2( Vc::Zero );
1108  const float_v x = 0.5f * log( bg2 );
1109  const float_v lhwI = log( 28.816f * 1e-9f * sqrt( rho * mZA ) / mI );
1110 
1111  float_m init = x > x1;
1112 
1113  d2(init) = lhwI + x - 0.5f;
1114  const float_v &r = ( x1 - x ) / ( x1 - x0 );
1115  init = (x > x0) & (x1 > x);
1116  d2(init) = (lhwI + x - 0.5f + ( 0.5f - lhwI - x0 ) * r * r * r);
1117 
1118  return mK*mZA*( 1.f + bg2 ) / bg2*( 0.5f*log( _2me*bg2*maxT/(mI*mI) ) - bg2 / ( 1.f + bg2 ) - d2 );
1119 }
1120 
1121 inline void PndFTSCATrackParamVector::EnergyLossCorrection( const float_v& mass2, const float_v& radThick, float_v &qp0, float_v direction, const float_m &mask)
1122 {
1123  const float_v& p2 = 1.f/(qp0*qp0);
1124  const float_v& E2 = mass2 + p2;
1125 
1126  const float_v& bethe = ApproximateBetheBloch( p2/mass2 );
1127 
1128  float_v tr = sqrt(1.f + fTx*fTx + fTy*fTy) ;
1129 
1130  const float_v& dE = bethe * radThick*tr * 2.33f * 9.34961f;
1131 
1132  const float_v& E2Corrected = (sqrt(E2) + direction*dE) * (sqrt(E2) + direction*dE);
1133  float_v corr = sqrt( p2/( E2Corrected - mass2 ) );
1134  //cout<<"corr "<<corr<<endl;
1135  float_m init = !(corr == corr) || !(mask) ;
1136  corr(init) = 1.f;
1137  qp0(mask) *= corr;
1138  fQP(mask) *= corr;
1139  fC40(mask) *= corr;
1140  fC41(mask) *= corr;
1141  fC42(mask) *= corr;
1142  fC43(mask) *= corr;
1143  fC44(mask) *= corr * corr;
1144 }
1145 
1146 inline float_m PndFTSCATrackParamVector::FilterVtx( const float_v &xV, const float_v& yV, const L1XYMeasurementInfo &info, float_v& extrDx, float_v& extrDy, float_v J[], const float_m &active )
1147 {
1148  float_v zeta0, zeta1, S00, S10, S11, si;
1149  float_v F00, F10, F20, F30, F40, F01, F11, F21, F31, F41 ;
1150  float_v K00, K10, K20, K30, K40, K01, K11, K21, K31, K41;
1151 
1152  zeta0 = fX + extrDx - xV;
1153  zeta1 = fY + extrDy - yV;
1154 
1155  // H = 1 0 J[0] J[1] J[2]
1156  // 0 1 J[3] J[4] J[5]
1157 
1158  // F = CH'
1159  F00 = fC00; F01 = fC10;
1160  F10 = fC10; F11 = fC11;
1161  F20 = J[0]*fC22; F21 = J[3]*fC22;
1162  F30 = J[1]*fC33; F31 = J[4]*fC33;
1163  F40 = J[2]*fC44; F41 = J[5]*fC44;
1164 
1165  S00 = info.C00 + F00 + J[0]*F20 + J[1]*F30 + J[2]*F40;
1166  S10 = info.C10 + F10 + J[3]*F20 + J[4]*F30 + J[5]*F40;
1167  S11 = info.C11 + F11 + J[3]*F21 + J[4]*F31 + J[5]*F41;
1168 
1169  si = 1.f/(S00*S11 - S10*S10);
1170  float_v S00tmp = S00;
1171  S00 = si*S11;
1172  S10 = -si*S10;
1173  S11 = si*S00tmp;
1174 
1175  fChi2(active) += zeta0*zeta0*S00 + 2.f*zeta0*zeta1*S10 + zeta1*zeta1*S11;
1176 
1177  K00 = F00*S00 + F01*S10; K01 = F00*S10 + F01*S11;
1178  K10 = F10*S00 + F11*S10; K11 = F10*S10 + F11*S11;
1179  K20 = F20*S00 + F21*S10; K21 = F20*S10 + F21*S11;
1180  K30 = F30*S00 + F31*S10; K31 = F30*S10 + F31*S11;
1181  K40 = F40*S00 + F41*S10; K41 = F40*S10 + F41*S11;
1182 
1183  fX (active) -= K00*zeta0 + K01*zeta1;
1184  fY (active) -= K10*zeta0 + K11*zeta1;
1185  fTx(active) -= K20*zeta0 + K21*zeta1;
1186  fTy(active) -= K30*zeta0 + K31*zeta1;
1187  fQP(active) -= K40*zeta0 + K41*zeta1;
1188 
1189  fC00(active) -= ( K00*F00 + K01*F01 );
1190  fC10(active) -= ( K10*F00 + K11*F01 );
1191  fC11(active) -= ( K10*F10 + K11*F11 );
1192  fC20(active) = -( K20*F00 + K21*F01 );
1193  fC21(active) = -( K20*F10 + K21*F11 );
1194  fC22(active) -= ( K20*F20 + K21*F21 );
1195  fC30(active) = -( K30*F00 + K31*F01 );
1196  fC31(active) = -( K30*F10 + K31*F11 );
1197  fC32(active) = -( K30*F20 + K31*F21 );
1198  fC33(active) -= ( K30*F30 + K31*F31 );
1199  fC40(active) = -( K40*F00 + K41*F01 );
1200  fC41(active) = -( K40*F10 + K41*F11 );
1201  fC42(active) = -( K40*F20 + K41*F21 );
1202  fC43(active) = -( K40*F30 + K41*F31 );
1203  fC44(active) -= ( K40*F40 + K41*F41 );
1204 
1205  return active;
1206 }
1207 
1208 inline float_m PndFTSCATrackParamVector::TransportJXY0ToX0( const float_v &x0, const L1FieldRegion &F, float_v& extrDx, float_v& extrDy, float_v &J04, float_v &J14, const float_m &active )
1209 {
1210  const float_v c_light = 0.000299792458f, c1 = 1.f, c2i = 0.5f, c6i = 1.f/6.f, c12i = 1.f/12.f;
1211 
1212  const float_v dz = x0 - fZ;
1213  float_v dz2 = dz*dz;
1214  float_v dzc6i = dz*c6i;
1215  float_v dz2c12i = dz2*c12i;
1216 
1217  float_v xx = fTx*fTx;
1218  float_v yy = fTy*fTy;
1219  float_v xy = fTx*fTy;
1220 
1221  float_v Ay = -xx-c1;
1222  float_v bx = yy+c1;
1223 
1224  float_v ctdz2 = c_light*sqrt( c1 + xx + yy )*dz2;
1225 
1226  float_v Sx = F.cx0*c2i + F.cx1*dzc6i + F.cx2*dz2c12i ;
1227  float_v Sy = F.cy0*c2i + F.cy1*dzc6i + F.cy2*dz2c12i ;
1228  float_v Sz = F.cz0*c2i + F.cz1*dzc6i + F.cz2*dz2c12i ;
1229 
1230  extrDx(active) = ( fTx )*dz ;
1231  extrDy(active) = ( fTy )*dz ;
1232  J04(active) = ctdz2 * (Sx*xy + Sy*Ay + Sz*fTy);
1233  J14(active) = ctdz2 * (Sx*bx - Sy*xy - Sz*fTx);
1234 
1235  return active;
1236 }
1237 
1238 #else // PANDA_FTS
1239 
1240 #include "PndFTSCADef.h"
1241 #include "PndFTSVector.h"
1242 #include "PndFTSCAMath.h"
1243 #include "CAX1X2MeasurementInfo.h"
1244 #include "FTSCAStation.h"
1245 
1246 class PndFTSCATrackLinearisationVector;
1247 class PndFTSCAParam;
1248 class FTSCAHit;
1249 class FTSCAHitV;
1250 class FTSCATarget;
1251 
1252 
1253 namespace std
1254 {
1255  template<typename T> struct char_traits;
1256  template<typename _CharT, typename _Traits> class basic_istream;
1257  typedef basic_istream<char, char_traits<char> > istream;
1258  template<typename _CharT, typename _Traits> class basic_ostream;
1259  typedef basic_ostream<char, char_traits<char> > ostream;
1260 } // namespace std
1261 
1269 class PndFTSCATrackParamVector
1270 {
1271  friend std::istream &operator>>( std::istream &, PndFTSCATrackParamVector & );
1272  friend std::ostream &operator<<( std::ostream &, const PndFTSCATrackParamVector & );
1273  public:
1274  PndFTSCATrackParamVector()
1275  : fX( Vc::Zero ),
1276  fSignCosPhi( Vc::Zero ),
1277  fChi2( Vc::Zero ),
1278  fNDF( Vc::Zero )
1279  {
1280  for ( int i = 0; i < 5; ++i ) fP[i].setZero();
1281  for ( int i = 0; i < 15; ++i ) fC[i].setZero();
1282  }
1283 
1284  void ConvertTrackParamToVector( PndFTSCATrackParam t0[float_v::Size], int nTracksV );
1285 
1286  void InitCovMatrix( float_v d2QMom = 0.f );
1287  void InitByTarget( const FTSCATarget& target );
1288  void InitByHit( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_v& dQMom );
1289 
1290  void InitDirection( float_v r0, float_v r1, float_v r2 ) // initialize direction parameters according to a given tangent vector
1291  {
1292  const float_v r = sqrt( r0*r0+r1*r1 );
1293  SetSinPhi( r1/r );
1294  SetSignCosPhi( r0/abs(r0) );
1295  SetDzDs( r2/r );
1296  }
1297 
1298  struct PndFTSCATrackFitParam {
1299  float_v fBethe;
1300  float_v fE;
1301  float_v fTheta2;
1302  float_v fEP2;
1303  float_v fSigmadE2;
1304  float_v fK22;
1305  float_v fK33;
1306  float_v fK43;
1307  float_v fK44;
1308  };
1309 
1310  float_v X() const { return fX; }
1311  float_v Y() const { return fP[0]; }
1312  float_v Z() const { return fP[1]; }
1313  float_v SinPhi() const { return fP[2]; }
1314  float_v DzDs() const { return fP[3]; }
1315  float_v QPt() const { return fP[4]; }
1316 
1317  float_v Angle() const { return fAlpha; }
1318 
1319  float_v X0() const { return X(); }
1320  float_v X1() const { return Y(); }
1321  float_v X2() const { return Z(); }
1322  float_v Tx1() const { return SinPhi()/(SignCosPhi()*sqrt( 1 - SinPhi()*SinPhi() )); } // CHECKME
1323  float_v Tx2() const { return DzDs()/(SignCosPhi()*sqrt( 1 - SinPhi()*SinPhi() )); } // dx2/dx0 = dz/dx
1324  float_v QMomentum() const { return QPt(); } // used for triplets comparison
1325 
1330  float_v SignCosPhi() const { return fSignCosPhi; }
1331  float_v Chi2() const { return fChi2; }
1332  int_v NDF() const { return fNDF; }
1333 
1334  float_v Err2Y() const { return fC[0]; }
1335  float_v Err2Z() const { return fC[2]; }
1336  float_v Err2SinPhi() const { return fC[5]; }
1337  float_v Err2DzDs() const { return fC[9]; }
1338  float_v Err2QPt() const { return fC[14]; }
1339 
1340  float_v GetX() const { return fX; }
1341  float_v GetY() const { return fP[0]; }
1342  float_v GetZ() const { return fP[1]; }
1343  float_v GetSinPhi() const { return fP[2]; }
1344  float_v GetDzDs() const { return fP[3]; }
1345  float_v GetQPt() const { return fP[4]; }
1346  float_v GetSignCosPhi() const { return fSignCosPhi; }
1347  float_v GetChi2() const { return fChi2; }
1348  int_v GetNDF() const { return fNDF; }
1349 
1350  float_v GetKappa( const float_v &Bz ) const { return fP[4]*Bz; }
1351  float_v GetCosPhiPositive() const { return CAMath::Sqrt( float_v( Vc::One ) - SinPhi()*SinPhi() ); }
1352  float_v GetCosPhi() const { return fSignCosPhi*CAMath::Sqrt( float_v( Vc::One ) - SinPhi()*SinPhi() ); }
1353 
1354  float_v Err2X1() const { return fC[0]; }
1355  float_v Err2X2() const { return fC[2]; }
1356  // float_v GetErr2SinPhi() const { return fC[5]; }
1357  // float_v GetErr2DzDs() const { return fC[9]; }
1358  float_v Err2QMomentum() const { return fC[14]; }
1359 
1360  const float_v& Par( int i ) const { return fP[i]; }
1361  const float_v& Cov( int i ) const { return fC[i]; }
1362 
1363  private:
1364  friend class PndFTSCATrackParam;
1365  const float_v *Par() const { return fP; }
1366  const float_v *Cov() const { return fC; }
1367  float_v *Par() { return fP; }
1368  float_v *Cov() { return fC; }
1369  public:
1370 
1371  void SetTrackParam(const PndFTSCATrackParamVector &param, const float_m &m = float_m( true ) )
1372  {
1373  for(int i=0; i<5; i++) fP[i](m) = param.Par()[i];
1374  for(int i=0; i<15; i++) fC[i](m) = param.Cov()[i];
1375  fX(m) = param.X();
1376  fSignCosPhi(m) = param.SignCosPhi();
1377  fChi2(m) = param.GetChi2();
1378  fNDF(static_cast<int_m>(m)) = param.GetNDF();
1379  fAlpha(static_cast<int_m>(m)) = param.Angle();
1380  }
1381 
1382  void SetTrackParamOne(int iV, const PndFTSCATrackParamVector &param, int iVa )
1383  {
1384  for(int i=0; i<5; i++) fP[i][iV] = param.Par()[i][iVa];
1385  for(int i=0; i<15; i++) fC[i][iV] = param.Cov()[i][iVa];
1386  fX[iV] = param.X()[iVa];
1387  fSignCosPhi[iV] = param.SignCosPhi()[iVa];
1388  fChi2[iV] = param.GetChi2()[iVa];
1389  fNDF[iV] = param.GetNDF()[iVa];
1390  fAlpha[iV] = param.Angle()[iVa];
1391  }
1392 
1393  void SetPar( int i, const float_v &v ) { fP[i] = v; }
1394  void SetPar( int i, const float_v &v, const float_m &m ) { fP[i]( m ) = v; }
1395  void SetCov( int i, const float_v &v ) { fC[i] = v; }
1396  void SetCov( int i, const float_v &v, const float_m &m ) { fC[i]( m ) = v; }
1397 
1398  void SetX( const float_v &v ) { fX = v; }
1399  void SetY( const float_v &v ) { fP[0] = v; }
1400  void SetZ( const float_v &v ) { fP[1] = v; }
1401  void SetX( const float_v &v, const float_m &m ) { fX( m ) = v; }
1402  void SetY( const float_v &v, const float_m &m ) { fP[0]( m ) = v; }
1403  void SetZ( const float_v &v, const float_m &m ) { fP[1]( m ) = v; }
1404  void SetSinPhi( const float_v &v ) { fP[2] = v; }
1405  void SetSinPhi( const float_v &v, const float_m &m ) { fP[2]( m ) = v; }
1406  void SetDzDs( const float_v &v ) { fP[3] = v; }
1407  void SetDzDs( const float_v &v, const float_m &m ) { fP[3]( m ) = v; }
1408  void SetQPt( const float_v &v ) { fP[4] = v; }
1409  void SetQMomentum( const float_v &v ) { SetQPt(v); }
1410  void SetQPt( const float_v &v, const float_m &m ) { fP[4]( m ) = v; }
1411  void SetSignCosPhi( const float_v &v ) { fSignCosPhi = v; }
1412  void SetSignCosPhi( const float_v &v, const float_m &m ) { fSignCosPhi(m) = v; }
1413  void SetChi2( const float_v &v ) { fChi2 = v; }
1414  void SetChi2( const float_v &v, const float_m &m ) { fChi2(m) = v; }
1415  void SetNDF( int v ) { fNDF = v; }
1416  void SetNDF( const int_v &v ) { fNDF = v; }
1417  void SetNDF( const int_v &v, const int_m &m ) { fNDF(m) = v; }
1418 
1419  void SetAngle( const float_v &v ) { fAlpha = v; }
1420  void SetAngle( const float_v &v, const float_m &m ) { fAlpha(m) = v; }
1421 
1422  void SetErr2Y( float_v v ) { fC[0] = v; }
1423  void SetErr2Z( float_v v ) { fC[2] = v; }
1424  void SetErr2QPt( float_v v ) { fC[14] = v; }
1425 
1426  float_v GetDist2( const PndFTSCATrackParamVector &t ) const;
1427  float_v GetDistXZ2( const PndFTSCATrackParamVector &t ) const;
1428 
1429 
1430  float_v GetS( const float_v &x, const float_v &y, const float_v &Bz ) const;
1431 
1432  void GetDCAPoint( const float_v &x, const float_v &y, const float_v &z,
1433  float_v *px, float_v *py, float_v *pz, const float_v &Bz ) const;
1434 
1435 
1436  float_m TransportToX0WithMaterial( const float_v &x, const float_v &XOverX0,
1437 const float_v &XThimesRho, const float_v &Bz, const float maxSinPhi = .999f );
1438 
1439  float_m TransportToX0( const float_v &x, const float_v &Bz, const float maxSinPhi = .999f, const float_m &mask = float_m( true ) );
1440 
1441  float_m TransportToX0( const float_v &x, PndFTSCATrackLinearisationVector &t0,
1442  const float_v &Bz, const float maxSinPhi = .999f, float_v *DL = 0, const float_m &mask = float_m( true ) );
1443 
1444  float_m TransportToX0( const float_v &x, const float_v &sinPhi0,
1445  const float_v &Bz, const float_v maxSinPhi = .999f, const float_m &mask = float_m( true ) );
1446 
1447  float_m TransportToX0WithMaterial( const float_v &x, PndFTSCATrackLinearisationVector &t0,
1448  PndFTSCATrackFitParam &par, const float_v &XOverX0,
1449  const float_v &XThimesRho,
1450  const float_v &Bz, const float maxSinPhi = .999f, const float_m &mask = float_m( true ) );
1451 
1452  float_m TransportToX0WithMaterial( const float_v &x,
1453  PndFTSCATrackFitParam &par, const float_v &XOverX0,
1454 const float_v &XThimesRho, const float_v &Bz, const float maxSinPhi = .999f );
1455 
1456  float_m Rotate( const float_v &alpha, PndFTSCATrackLinearisationVector &t0,
1457  const float maxSinPhi = .999f, const float_m &mask = float_m( true ) );
1458  float_m Rotate( const float_v &alpha, const float maxSinPhi = .999f, const float_m &mask = float_m( true ) );
1459  void RotateXY( float_v alpha, float_v &x, float_v &y, float_v &sin, const float_m &mask = float_m( true ) ) const ;
1460 
1461  float_m FilterWithMaterial( const float_v &y, const float_v &z, float_v err2Y, float_v errYZ, float_v err2Z,
1462  float maxSinPhi=0.999f, const float_m &mask = float_m( true ), const int_v& hitNDF = int_v(2) ,const float_v& chi2Cut = 10e10f ); // filters 2-D measurement
1463 
1464  float_m FilterWithMaterial( const float_v &y, const float_v &z, const FTSCAStripInfo &info, float_v err2,
1465  float maxSinPhi=0.999f, const float_m &mask = float_m( true ), const float_v& chi2Cut = 10e10f ); // filters 1-D measurement
1466  float_m FilterWithMaterial( const float_v &y, const float_v &z, const float_v &r, const FTSCAStripInfo &info, float_v err2,
1467  float maxSinPhi=0.999f, const float_m &mask = float_m( true ), const float_v& chi2Cut = 10e10f ); // filters tube measurement
1468 
1469  static float_v ApproximateBetheBloch( const float_v &beta2 );
1470  static float_v BetheBlochGeant( const float_v &bg,
1471  const float_v &kp0 = 2.33f,
1472  const float_v &kp1 = 0.20f,
1473  const float_v &kp2 = 3.00f,
1474  const float_v &kp3 = 173e-9f,
1475  const float_v &kp4 = 0.49848f
1476  );
1477  static float_v BetheBlochSolid( const float_v &bg );
1478  static float_v BetheBlochGas( const float_v &bg );
1479 
1480 
1481  void CalculateFitParameters( PndFTSCATrackFitParam &par, const float_v &mass = 0.13957f );
1482  float_m CorrectForMeanMaterial( const float_v &xOverX0, const float_v &xTimesRho,
1483  const PndFTSCATrackFitParam &par, const float_m &_mask );
1484 
1485  // float_m FilterDelta( const float_m &mask, const float_v &dy, const float_v &dz,
1486  // float_v err2Y, float_v err2Z, const float maxSinPhi = .999f );
1487  // float_m Filter( const float_m &mask, const float_v &y, const float_v &z,
1488  // float_v err2Y, float_v err2Z, const float maxSinPhi = .999f );
1489 
1490 
1491  float_m FilterVtx( const float_v &xV, const float_v& yV, const CAX1X2MeasurementInfo &info, float_v& extrDx, float_v& extrDy, float_v J[], const float_m &active = float_m( true ) );
1492  float_m TransportJ0ToX0( const float_v &x0, const float_v&cBz, float_v& extrDy, float_v& extrDz, float_v J[], const float_m &active = float_m( true ) );
1493 
1494 
1495  float_m Transport( const int_v& ista, const PndFTSCAParam& param, const float_m &mask = float_m( true ) );
1496  float_m Transport( const FTSCAHitV& hit, const PndFTSCAParam& p, const float_m &mask = float_m( true ) );
1497  float_m Filter( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_m &mask = float_m( true ), const float_v& chi2Cut = 10e10f );
1498 
1499  float_m Transport( const FTSCAHit& hit, const PndFTSCAParam& p, const float_m &mask = float_m( true ) );
1500  float_m Filter( const FTSCAHit& hit, const PndFTSCAParam& param, const float_m &mask = float_m( true ), const float_v& chi2Cut = 10e10f );
1501 
1502  float_m AddTarget( const FTSCATarget& target, const float_m &mask = float_m( true ) );
1503 
1504  private:
1505 
1506  float_v fX; // x position
1507  float_v fSignCosPhi; // sign of cosPhi
1508  float_v fP[5]; // 'active' track parameters: Y, Z, SinPhi = dy/sqrt(dx*dx + dy*dy), DzDs = dz/sqrt(dx*dx + dy*dy), q/Pt // if q/pt = 0 then track equation is [dr x e] = 0, where e = { sqrt( ( 1 - SinPhi^2 )/( 1 + DzDs^2 ) ), sqrt( SinPhi^2/( 1 + DzDs^2 ) ), sqrt( 1 /( 1 + DzDs^2 ) ) }
1509  float_v fC[15]; // the covariance matrix for Y,Z,SinPhi,..
1510  float_v fChi2; // the chi^2 value
1511  int_v fNDF; // the Number of Degrees of Freedom
1512 
1513  float_v fAlpha; // coor system
1514 };
1515 
1516 #include "debug.h"
1517 
1518 inline float_m PndFTSCATrackParamVector::TransportToX0( const float_v &x, const float_v &sinPhi0,
1519  const float_v &Bz, const float_v maxSinPhi, const float_m &_mask )
1520 {
1521  //* Transport the track parameters to X=x, using linearization at phi0 with 0 curvature,
1522  //* and the field value Bz
1523  //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
1524  //* linearisation of trajectory t0 is also transported to X=x,
1525  //* returns 1 if OK
1526  //*
1527 
1528  debugKF() << "Start TransportToX0(" << x << ", " << _mask << ")\n" << *this << std::endl;
1529 
1530  const float_v &ey = sinPhi0;
1531  const float_v &dx = x - X();
1532  const float_v &exi = float_v( Vc::One ) * CAMath::RSqrt( float_v( Vc::One ) - ey * ey ); // RSqrt
1533 
1534  const float_v &dxBz = dx * Bz;
1535  const float_v &dS = dx * exi;
1536  const float_v &h2 = dS * exi * exi;
1537  const float_v &h4 = .5f * h2 * dxBz;
1538 //#define LOSE_DEBUG
1539 #ifdef LOSE_DEBUG
1540  std::cout << " TrTo-sinPhi0 = " << sinPhi0 << std::endl;
1541 #endif
1542 // const float_v &sinPhi = SinPhi() * (float_v( Vc::One ) - 0.5f * dxBz * QPt() *dxBz * QPt()/ ( float_v( Vc::One ) - SinPhi()*SinPhi() )) + dxBz * QPt();
1544  const float_v &sinPhi = SinPhi() + dxBz * QPt();
1546 #ifdef LOSE_DEBUG
1547  std::cout << " TrTo-sinPhi = " << sinPhi << std::endl;
1548 #endif
1549  float_m mask = _mask && CAMath::Abs( exi ) <= 1.e4f;
1550  mask &= ( (CAMath::Abs( sinPhi ) <= maxSinPhi) || (maxSinPhi <= 0.f) );
1551 
1552 
1553  fX ( mask ) += dx;
1554  fP[0]( mask ) += dS * ey + h2 * ( SinPhi() - ey ) + h4 * QPt();
1555  fP[1]( mask ) += dS * DzDs();
1556  fP[2]( mask ) = sinPhi;
1557 
1558 
1559  //const float_v c00 = fC[0];
1560  //const float_v c11 = fC[2];
1561  const float_v c20 = fC[3];
1562  //const float_v c21 = fC[4];
1563  const float_v c22 = fC[5];
1564  //const float_v c30 = fC[6];
1565  const float_v c31 = fC[7];
1566  //const float_v c32 = fC[8];
1567  const float_v c33 = fC[9];
1568  const float_v c40 = fC[10];
1569  //const float_v c41 = fC[11];
1570  const float_v c42 = fC[12];
1571  //const float_v c43 = fC[13];
1572  const float_v c44 = fC[14];
1573 
1574  const float_v two( 2.f );
1575 
1576  fC[0] ( mask ) += h2 * h2 * c22 + h4 * h4 * c44
1577  + two * ( h2 * c20 + h4 * ( c40 + h2 * c42 ) );
1578 
1579  //fC[1] ( mask ) += h2 * c21 + h4 * c41 + dS * ( c30 + h2 * c32 + h4 * c43 );
1580  fC[2] ( mask ) += dS * ( two * c31 + dS * c33 );
1581 
1582  fC[3] ( mask ) += h2 * c22 + h4 * c42 + dxBz * ( c40 + h2 * c42 + h4 * c44 );
1583  //fC[4] ( mask ) += dS * c32 + dxBz * ( c41 + dS * c43 );
1584  const float_v &dxBz_c44 = dxBz * c44;
1585  fC[12]( mask ) += dxBz_c44;
1586  fC[5] ( mask ) += dxBz * ( two * c42 + dxBz_c44 );
1587 
1588  //fC[6] ( mask ) += h2 * c32 + h4 * c43;
1589  fC[7] ( mask ) += dS * c33;
1590  //fC[8] ( mask ) += dxBz * c43;
1591  //fC[9] ( mask ) = c33;
1592 
1593  fC[10]( mask ) += h2 * c42 + h4 * c44;
1594  //fC[11]( mask ) += dS * c43;
1595  //fC[13]( mask ) = c43;
1596  //fC[14]( mask ) = c44;
1597 
1598  debugKF() << mask << "\n" << *this << std::endl;
1599  return mask;
1600 }
1601 
1602 #include <assert.h>
1603 
1604 #ifdef NVALGRIND
1605 #define VALGRIND_CHECK_VALUE_IS_DEFINED( mask )
1606 #define VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( v, k )
1607 #define VALGRIND_CHECK_MEM_IS_DEFINED( v, k );
1608 #define VALGRIND_CHECK_MEM_IS_ADDRESSABLE( v, k );
1609 #else
1610 #include <valgrind/memcheck.h>
1611 #define VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( v, k ) \
1612 { \
1613  __typeof__( v + v ) tmp( v ); \
1614  tmp.setZero( !k ); \
1615  VALGRIND_CHECK_VALUE_IS_DEFINED( tmp ); \
1616 }
1617 #endif
1618 
1619 // inline float_m PndFTSCATrackParamVector::FilterDelta( const float_m &mask, const float_v &dy, const float_v &dz,
1620 // float_v err2Y, float_v err2Z, const float maxSinPhi )
1621 // {
1622 // debugKF() << "Kalman filter( " << mask
1623 // << "\n " << dy
1624 // << "\n " << dz
1625 // << "\n " << err2Y
1626 // << "\n " << err2Z
1627 // << "\n):" << std::endl;
1628 
1629 // assert( err2Y > 0.f || !mask );
1630 // assert( err2Z > 0.f || !mask );
1631 // VALGRIND_CHECK_VALUE_IS_DEFINED( mask );
1632 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( dy, mask );
1633 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( dz, mask );
1634 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( err2Y, mask );
1635 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( err2Z, mask );
1636 // #ifndef NVALGRIND
1637 // err2Y.setZero( !mask );
1638 // err2Z.setZero( !mask );
1639 // #endif
1640 // VALGRIND_CHECK_VALUE_IS_DEFINED( maxSinPhi );
1641 
1642 // //* Add the y,z measurement with the Kalman filter
1643 
1644 // const float_v c00 = fC[ 0];
1645 // const float_v c11 = fC[ 2];
1646 // const float_v c20 = fC[ 3];
1647 // const float_v c31 = fC[ 7];
1648 // const float_v c40 = fC[10];
1649 
1650 // VALGRIND_CHECK_VALUE_IS_DEFINED( c00 );
1651 // VALGRIND_CHECK_VALUE_IS_DEFINED( c11 );
1652 // VALGRIND_CHECK_VALUE_IS_DEFINED( c20 );
1653 // VALGRIND_CHECK_VALUE_IS_DEFINED( c40 );
1654 // VALGRIND_CHECK_VALUE_IS_DEFINED( c31 );
1655 
1656 // err2Y += c00;
1657 // err2Z += c11;
1658 // #ifndef NODEBUG
1659 // if ( !( err2Y > 0.f || !mask ).isFull() ) {
1660 // std::cerr << err2Y << mask << ( err2Y > 0.f || !mask ) << c00 << std::endl;
1661 // }
1662 // #endif
1663 // assert( err2Y > 0.f || !mask );
1664 // assert( err2Z > 0.f || !mask );
1665 
1666 // VALGRIND_CHECK_VALUE_IS_DEFINED( fP[0] );
1667 // VALGRIND_CHECK_VALUE_IS_DEFINED( fP[1] );
1668 // VALGRIND_CHECK_VALUE_IS_DEFINED( fP[2] );
1669 
1670 // const float_v &z0 = dy;
1671 // const float_v &z1 = dz;
1672 
1673 // const float_v &mS0 = float_v( Vc::One ) / err2Y;
1674 // const float_v &mS2 = float_v( Vc::One ) / err2Z;
1675 // //const float_v &mS0 = CAMath::Reciprocal( err2Y );
1676 // //const float_v &mS2 = CAMath::Reciprocal( err2Z );
1677 // debugKF() << "delta(mS0): " << CAMath::Abs( float_v( Vc::One ) / err2Y - mS0 ) << std::endl;
1678 // debugKF() << "delta(mS2): " << CAMath::Abs( float_v( Vc::One ) / err2Z - mS2 ) << std::endl;
1679 // assert( mS0 > 0.f || !mask );
1680 // assert( mS2 > 0.f || !mask );
1681 
1682 // // K = CHtS
1683 
1684 // const float_v &k00 = c00 * mS0;
1685 // const float_v &k20 = c20 * mS0;
1686 // const float_v &k40 = c40 * mS0;
1687 
1688 // const float_v &k11 = c11 * mS2;
1689 // const float_v &k31 = c31 * mS2;
1690 
1691 // debugKF() << "delta(k00): " << ( c00 / err2Y - k00 ) << std::endl;
1692 // debugKF() << "delta(k20): " << ( c20 / err2Y - k20 ) << std::endl;
1693 // debugKF() << "delta(k40): " << ( c40 / err2Y - k40 ) << std::endl;
1694 
1695 // debugKF() << "delta(k11): " << ( c11 / err2Z - k11 ) << std::endl;
1696 // debugKF() << "delta(k31): " << ( c31 / err2Z - k31 ) << std::endl;
1697 
1698 // const float_v &sinPhi = fP[2] + k20 * z0 ;
1699 // debugKF() << "delta(sinPhi): " << ( z0 * c20 / err2Y + fP[2] - sinPhi ) << std::endl;
1700 
1701 // assert( maxSinPhi > 0.f );
1702 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( sinPhi, mask );
1703 // const float_m &success = mask && err2Y >= 1.e-8f && err2Z >= 1.e-8f && CAMath::Abs( sinPhi ) < maxSinPhi;
1704 // VALGRIND_CHECK_VALUE_IS_DEFINED( success );
1705 
1706 // fNDF ( static_cast<int_m>( success ) ) += 2;
1707 // fChi2 ( success ) += mS0 * z0 * z0 + mS2 * z1 * z1 ;
1708 
1709 // fP[ 0]( success ) += k00 * z0 ;
1710 // fP[ 1]( success ) += k11 * z1 ;
1711 // fP[ 2]( success ) = sinPhi ;
1712 // fP[ 3]( success ) += k31 * z1 ;
1713 // fP[ 4]( success ) += k40 * z0 ;
1714 
1715 // fC[ 0]( success ) -= k00 * c00 ;
1716 // fC[ 3]( success ) -= k20 * c00 ;
1717 // fC[ 5]( success ) -= k20 * c20 ;
1718 // fC[10]( success ) -= k40 * c00 ;
1719 // fC[12]( success ) -= k40 * c20 ;
1720 // fC[14]( success ) -= k40 * c40 ;
1721 
1722 // fC[ 2]( success ) -= k11 * c11 ;
1723 // fC[ 7]( success ) -= k31 * c11 ;
1724 // fC[ 9]( success ) -= k31 * c31 ;
1725 // #if 1
1726 // const float_m check = ( fC[ 0] >= 0.f ) && ( fC[ 2] >= 0.f ) && ( fC[ 5] >= 0.f ) && ( fC[ 9] >= 0.f ) && ( fC[14] >= 0.f );
1727 // #else
1728 // assert( fC[ 0] >= 0.f );
1729 // assert( fC[ 2] >= 0.f );
1730 // assert( fC[ 5] >= 0.f );
1731 // assert( fC[ 9] >= 0.f );
1732 // assert( fC[14] >= 0.f );
1733 // #endif
1734 // return success && check;
1735 // }
1736 
1737 // inline float_m PndFTSCATrackParamVector::Filter( const float_m &mask, const float_v &y, const float_v &z,
1738 // float_v err2Y, float_v err2Z, const float maxSinPhi )
1739 // {
1740 // debugKF() << "Kalman filter( " << mask
1741 // << "\n " << y
1742 // << "\n " << z
1743 // << "\n " << err2Y
1744 // << "\n " << err2Z
1745 // << "\n):" << std::endl;
1746 
1747 // assert( err2Y > 0.f || !mask );
1748 // assert( err2Z > 0.f || !mask );
1749 // VALGRIND_CHECK_VALUE_IS_DEFINED( mask );
1750 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( y, mask );
1751 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( z, mask );
1752 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( err2Y, mask );
1753 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( err2Z, mask );
1754 // #ifndef NVALGRIND
1755 // err2Y.setZero( !mask );
1756 // err2Z.setZero( !mask );
1757 // #endif
1758 // VALGRIND_CHECK_VALUE_IS_DEFINED( maxSinPhi );
1759 
1760 // //* Add the y,z measurement with the Kalman filter
1761 
1762 // const float_v c00 = fC[ 0];
1763 // const float_v c11 = fC[ 2];
1764 // const float_v c20 = fC[ 3];
1765 // const float_v c31 = fC[ 7];
1766 // const float_v c40 = fC[10];
1767 
1768 // VALGRIND_CHECK_VALUE_IS_DEFINED( c00 );
1769 // VALGRIND_CHECK_VALUE_IS_DEFINED( c11 );
1770 // VALGRIND_CHECK_VALUE_IS_DEFINED( c20 );
1771 // VALGRIND_CHECK_VALUE_IS_DEFINED( c40 );
1772 // VALGRIND_CHECK_VALUE_IS_DEFINED( c31 );
1773 
1774 // err2Y += c00;
1775 // err2Z += c11;
1776 // #ifndef NODEBUG
1777 // if ( !( err2Y > 0.f || !mask ).isFull() ) {
1778 // std::cerr << err2Y << mask << ( err2Y > 0.f || !mask ) << c00 << std::endl;
1779 // }
1780 // #endif
1781 // assert( err2Y > 0.f || !mask );
1782 // assert( err2Z > 0.f || !mask );
1783 
1784 // VALGRIND_CHECK_VALUE_IS_DEFINED( fP[0] );
1785 // VALGRIND_CHECK_VALUE_IS_DEFINED( fP[1] );
1786 // VALGRIND_CHECK_VALUE_IS_DEFINED( fP[2] );
1787 
1788 // const float_v &z0 = y - fP[0];
1789 // const float_v &z1 = z - fP[1];
1790 
1791 // const float_v &mS0 = float_v( Vc::One ) / err2Y;
1792 // const float_v &mS2 = float_v( Vc::One ) / err2Z;
1793 // //const float_v &mS0 = CAMath::Reciprocal( err2Y );
1794 // //const float_v &mS2 = CAMath::Reciprocal( err2Z );
1795 // debugKF() << "delta(mS0): " << CAMath::Abs( float_v( Vc::One ) / err2Y - mS0 ) << std::endl;
1796 // debugKF() << "delta(mS2): " << CAMath::Abs( float_v( Vc::One ) / err2Z - mS2 ) << std::endl;
1797 // assert( mS0 > 0.f || !mask );
1798 // assert( mS2 > 0.f || !mask );
1799 
1800 // // K = CHtS
1801 
1802 // const float_v &k00 = c00 * mS0;
1803 // const float_v &k20 = c20 * mS0;
1804 // const float_v &k40 = c40 * mS0;
1805 
1806 // const float_v &k11 = c11 * mS2;
1807 // const float_v &k31 = c31 * mS2;
1808 
1809 // debugKF() << "delta(k00): " << ( c00 / err2Y - k00 ) << std::endl;
1810 // debugKF() << "delta(k20): " << ( c20 / err2Y - k20 ) << std::endl;
1811 // debugKF() << "delta(k40): " << ( c40 / err2Y - k40 ) << std::endl;
1812 
1813 // debugKF() << "delta(k11): " << ( c11 / err2Z - k11 ) << std::endl;
1814 // debugKF() << "delta(k31): " << ( c31 / err2Z - k31 ) << std::endl;
1815 
1816 // const float_v &sinPhi = fP[2] + k20 * z0 ;
1817 // debugKF() << "delta(sinPhi): " << ( z0 * c20 / err2Y + fP[2] - sinPhi ) << std::endl;
1818 
1819 // assert( maxSinPhi > 0.f );
1820 // VALGRIND_CHECK_MASKED_VECTOR_IS_DEFINED( sinPhi, mask );
1821 // const float_m &success = mask && err2Y >= 1.e-8f && err2Z >= 1.e-8f && CAMath::Abs( sinPhi ) < maxSinPhi;
1822 // VALGRIND_CHECK_VALUE_IS_DEFINED( success );
1823 
1824 // fNDF ( static_cast<int_m>( success ) ) += 2;
1825 // fChi2 ( success ) += mS0 * z0 * z0 + mS2 * z1 * z1 ;
1826 
1827 // fP[ 0]( success ) += k00 * z0 ;
1828 // fP[ 1]( success ) += k11 * z1 ;
1829 // fP[ 2]( success ) = sinPhi ;
1830 // fP[ 3]( success ) += k31 * z1 ;
1831 // fP[ 4]( success ) += k40 * z0 ;
1832 
1833 // fC[ 0]( success ) -= k00 * c00 ;
1834 // fC[ 3]( success ) -= k20 * c00 ;
1835 // fC[ 5]( success ) -= k20 * c20 ;
1836 // fC[10]( success ) -= k40 * c00 ;
1837 // fC[12]( success ) -= k40 * c20 ;
1838 // fC[14]( success ) -= k40 * c40 ;
1839 
1840 // fC[ 2]( success ) -= k11 * c11 ;
1841 // fC[ 7]( success ) -= k31 * c11 ;
1842 // fC[ 9]( success ) -= k31 * c31 ;
1843 // #if 1
1844 // const float_m check = ( fC[ 0] >= 0.f ) && ( fC[ 2] >= 0.f ) && ( fC[ 5] >= 0.f ) && ( fC[ 9] >= 0.f ) && ( fC[14] >= 0.f );
1845 // #else
1846 // assert( fC[ 0] >= 0.f );
1847 // assert( fC[ 2] >= 0.f );
1848 // assert( fC[ 5] >= 0.f );
1849 // assert( fC[ 9] >= 0.f );
1850 // assert( fC[14] >= 0.f );
1851 // #endif
1852 // return success && check;
1853 // }
1854 
1855 inline float_m PndFTSCATrackParamVector::Rotate( const float_v &alpha, const float maxSinPhi, const float_m &mask )
1856 {
1857  //* Rotate the coordinate system in XY on the angle alpha
1858  if ( (abs(alpha) < 1e-6f || !mask).isFull() ) return mask;
1859 
1860  const float_v cA = CAMath::Cos( alpha );
1861  const float_v sA = CAMath::Sin( alpha );
1862  const float_v x = X(), y = Y(), sP = SinPhi(), cP = GetCosPhi();
1863  const float_v cosPhi = cP * cA + sP * sA;
1864  const float_v sinPhi = -cP * sA + sP * cA;
1865 
1866  float_m mReturn = mask && (CAMath::Abs( sinPhi ) < maxSinPhi) && (CAMath::Abs( cosPhi ) > 1.e-2f) && (CAMath::Abs( cP ) > 1.e-2f);
1867 
1868  mReturn &= abs(alpha) < 3.1415f * 0.25f; // allow turn by 45 degree only
1869 
1870  const float_v j0 = cP / cosPhi;
1871  const float_v j2 = cosPhi / cP;
1872 
1873  SetX( x*cA + y*sA, mReturn);
1874  SetY( -x*sA + y*cA, mReturn );
1875  SetSignCosPhi( CAMath::Abs(cosPhi)/cosPhi, mReturn );
1876  SetSinPhi( sinPhi, mReturn );
1877 
1878 
1879  //float J[5][5] = { { j0, 0, 0, 0, 0 }, // Y
1880  // { 0, 1, 0, 0, 0 }, // Z
1881  // { 0, 0, j2, 0, 0 }, // SinPhi
1882  // { 0, 0, 0, 1, 0 }, // DzDs
1883  // { 0, 0, 0, 0, 1 } }; // Kappa
1884  //cout<<"alpha="<<alpha<<" "<<x<<" "<<y<<" "<<sP<<" "<<cP<<" "<<j0<<" "<<j2<<endl;
1885  //cout<<" "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl;
1886  fC[0](mReturn) *= j0 * j0;
1887  fC[1](mReturn) *= j0;
1888  fC[3](mReturn) *= j0;
1889  fC[6](mReturn) *= j0;
1890  fC[10](mReturn) *= j0;
1891 
1892  fC[3](mReturn) *= j2;
1893  fC[4](mReturn) *= j2;
1894  fC[5](mReturn) *= j2 * j2;
1895  fC[8](mReturn) *= j2;
1896  fC[12](mReturn) *= j2;
1897 
1898  fAlpha(mReturn) += alpha;
1899  //cout<<" "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl;
1900  return mReturn;
1901 }
1902 
1903 inline void PndFTSCATrackParamVector::RotateXY( float_v alpha, float_v &x, float_v &y, float_v &sin, const float_m &mask ) const
1904 {
1905  //* Rotate the coordinate system in XY on the angle alpha
1906  if ( (abs(alpha) < 1e-6f || !mask).isFull() ) return;
1907 
1908  const float_v cA = CAMath::Cos( alpha );
1909  const float_v sA = CAMath::Sin( alpha );
1910 
1911  x(mask) = ( X()*cA + Y()*sA );
1912  y(mask) = ( -X()*sA + Y()*cA );
1913  sin(mask) = -GetCosPhi() * sA + SinPhi() * cA;
1914 }
1915 
1916 
1917 inline float_m PndFTSCATrackParamVector::FilterVtx( const float_v &yV, const float_v& zV, const CAX1X2MeasurementInfo &info, float_v& extrDy, float_v& extrDz, float_v J[], const float_m &active )
1918 {
1919  float_v zeta0, zeta1, S00, S10, S11, si;
1920  float_v F00, F10, F20, F30, F40, F01, F11, F21, F31, F41 ;
1921  float_v K00, K10, K20, K30, K40, K01, K11, K21, K31, K41;
1922  float_v& c00 = fC[0];
1923  float_v& c10 = fC[1];
1924  float_v& c11 = fC[2];
1925  float_v& c20 = fC[3];
1926  float_v& c21 = fC[4];
1927  float_v& c22 = fC[5];
1928  float_v& c30 = fC[6];
1929  float_v& c31 = fC[7];
1930  float_v& c32 = fC[8];
1931  float_v& c33 = fC[9];
1932  float_v& c40 = fC[10];
1933  float_v& c41 = fC[11];
1934  float_v& c42 = fC[12];
1935  float_v& c43 = fC[13];
1936  float_v& c44 = fC[14];
1937 
1938  zeta0 = Y() + extrDy - yV;
1939  zeta1 = Z() + extrDz - zV;
1940 
1941  // H = 1 0 J[0] J[1] J[2]
1942  // 0 1 J[3] J[4] J[5]
1943 
1944  // F = CH'
1945  F00 = c00; F01 = c10;
1946  F10 = c10; F11 = c11;
1947  F20 = J[0]*c22; F21 = J[3]*c22;
1948  F30 = J[1]*c33; F31 = J[4]*c33;
1949  F40 = J[2]*c44; F41 = J[5]*c44;
1950 
1951  S00 = info.C00() + F00 + J[0]*F20 + J[1]*F30 + J[2]*F40;
1952  S10 = info.C10() + F10 + J[3]*F20 + J[4]*F30 + J[5]*F40;
1953  S11 = info.C11() + F11 + J[3]*F21 + J[4]*F31 + J[5]*F41;
1954 
1955  si = 1.f/(S00*S11 - S10*S10);
1956  float_v S00tmp = S00;
1957  S00 = si*S11;
1958  S10 = -si*S10;
1959  S11 = si*S00tmp;
1960 
1961  fChi2(active) += zeta0*zeta0*S00 + 2.f*zeta0*zeta1*S10 + zeta1*zeta1*S11;
1962 
1963  K00 = F00*S00 + F01*S10; K01 = F00*S10 + F01*S11;
1964  K10 = F10*S00 + F11*S10; K11 = F10*S10 + F11*S11;
1965  K20 = F20*S00 + F21*S10; K21 = F20*S10 + F21*S11;
1966  K30 = F30*S00 + F31*S10; K31 = F30*S10 + F31*S11;
1967  K40 = F40*S00 + F41*S10; K41 = F40*S10 + F41*S11;
1968 
1969  fP[0](active) -= K00*zeta0 + K01*zeta1;
1970  fP[1](active) -= K10*zeta0 + K11*zeta1;
1971  fP[2](active) -= K20*zeta0 + K21*zeta1;
1972  fP[3](active) -= K30*zeta0 + K31*zeta1;
1973  fP[4](active) -= K40*zeta0 + K41*zeta1;
1974 
1975  c00(active) -= ( K00*F00 + K01*F01 );
1976  c10(active) -= ( K10*F00 + K11*F01 );
1977  c11(active) -= ( K10*F10 + K11*F11 );
1978  c20(active) = -( K20*F00 + K21*F01 );
1979  c21(active) = -( K20*F10 + K21*F11 );
1980  c22(active) -= ( K20*F20 + K21*F21 );
1981  c30(active) = -( K30*F00 + K31*F01 );
1982  c31(active) = -( K30*F10 + K31*F11 );
1983  c32(active) = -( K30*F20 + K31*F21 );
1984  c33(active) -= ( K30*F30 + K31*F31 );
1985  c40(active) = -( K40*F00 + K41*F01 );
1986  c41(active) = -( K40*F10 + K41*F11 );
1987  c42(active) = -( K40*F20 + K41*F21 );
1988  c43(active) = -( K40*F30 + K41*F31 );
1989  c44(active) -= ( K40*F40 + K41*F41 );
1990 
1991  return active;
1992 }
1993 
1994 inline float_m PndFTSCATrackParamVector::TransportJ0ToX0( const float_v &x, const float_v& cBz, float_v& extrDy, float_v& extrDz, float_v J[], const float_m &active )
1995 {
1996  const float_v &ey = SinPhi();
1997  const float_v &dx = x - X();
1998  const float_v &exi = CAMath::RSqrt( float_v( Vc::One ) - ey * ey ); // RSqrt
1999 
2000  const float_v &dxBz = dx * cBz;
2001  const float_v &dS = dx * exi;
2002  const float_v &h2 = dS * exi * exi;
2003  const float_v &h4 = .5f * h2 * dxBz;
2004 
2005  float_m mask = active && CAMath::Abs( exi ) <= 1.e4f;
2006 
2007  extrDy(active) = dS*ey;
2008  extrDz(active) = dS*DzDs();
2009  J[0](active) = dS;
2010  J[1](active) = 0;
2011  J[2](active) = h4;
2012  J[3](active) = dS;
2013  J[4](active) = dS;
2014  J[5](active) = 0;
2015  return active;
2016 }
2017 
2018 std::istream &operator>>( std::istream &in, PndFTSCATrackParamVector &ot );
2019 std::ostream &operator<<( std::ostream &out, const PndFTSCATrackParamVector &ot );
2020 
2021 #endif // !PANDA_FTS
2022 
2023 typedef PndFTSCATrackParamVector TrackParamVector;
2024 
2025 #endif
PndFTSCATrackParamVector::BetheBlochGas
static float_v BetheBlochGas(const float_v &bg)
Definition: PndFTSCATrackParamVector.cxx:998
z0
Double_t z0
Definition: checkhelixhit.C:62
x0
Double_t x0
Definition: checkhelixhit.C:70
c5
c5
Definition: plot_dirc.C:75
p
Double_t p
Definition: anasim.C:58
L1FieldRegion::cy1
float_v cy1
Definition: L1Field.h:133
std::basic_istream
Definition: PndFTSCATrackParamVector.h:1256
mK
double mK
Definition: softtrigger_kin5.C:128
c11
TCanvas * c11
Definition: plot_acceptance_stt.C:188
t2
Int_t t2
Definition: hist-t7.C:106
PndFTSCATrackParamVector::SetCov
void SetCov(int i, const float_v &v)
Definition: PndFTSCATrackParamVector.h:1395
r
double r
Definition: RiemannTest.C:14
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
c4
c4
Definition: plot_dirc.C:71
PndFTSCATrackParamVector::SetAngle
void SetAngle(const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1420
PndFTSCATrackParamVector::InitByHit
void InitByHit(const FTSCAHitV &hit, const PndFTSCAParam &param, const float_v &dQMom)
Definition: PndFTSCATrackParamVector.cxx:652
PndFTSCATrackParamVector::Par
const float_v * Par() const
Definition: PndFTSCATrackParamVector.h:1365
i
Int_t i
Definition: run_full.C:25
m
__m128 m
Definition: P4_F32vec4.h:28
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
corr
PndPidCorrelator * corr
Definition: macro/examples/pid/run_pid_stt.C:49
PndFTSCATrackParamVector::CorrectForMeanMaterial
float_m CorrectForMeanMaterial(const float_v &xOverX0, const float_v &xTimesRho, const PndFTSCATrackFitParam &par, const float_m &_mask)
Definition: PndFTSCATrackParamVector.cxx:1072
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
mrfdata_8b_error::success
static const UInt_t success
Definition: MvdOfflineTBAnalysis_Topix4/mrfdata_8b.h:30
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
L1FieldRegion::cz1
float_v cz1
Definition: L1Field.h:134
PndFTSCATrackParamVector::SetSignCosPhi
void SetSignCosPhi(const float_v &v)
Definition: PndFTSCATrackParamVector.h:1411
L1FieldRegion::cz2
float_v cz2
Definition: L1Field.h:134
h4
TH1F * h4
Definition: anal_point_fwendcap.C:84
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
c10
TCanvas * c10
Definition: plot_invariantmass.C:42
tr
Double_t tr
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:19
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
par
Double_t par[3]
Definition: anaLmdCluster.C:146
log
friend F32vec4 log(const F32vec4 &a)
Definition: P4_F32vec4.h:110
c21
TCanvas * c21
Definition: plot_acceptance_stt.C:247
CAX1X2MeasurementInfo::C11
const float_v & C11() const
Definition: CAX1X2MeasurementInfo.h:18
c2
c2
Definition: plot_dirc.C:39
L1FieldRegion::cx1
float_v cx1
Definition: L1Field.h:132
operator<<
std::ostream & operator<<(std::ostream &o, const PndEventInfo &a)
Definition: PndEventInfo.cxx:68
PndFTSCATrackParamVector::SetTrackParamOne
void SetTrackParamOne(int iV, const PndFTSCATrackParamVector &param, int iVa)
Definition: PndFTSCATrackParamVector.h:1382
PndFTSCATrackParamVector::SetNDF
void SetNDF(const int_v &v, const int_m &m)
Definition: PndFTSCATrackParamVector.h:1417
PndFTSCATrackParamVector::SetSinPhi
void SetSinPhi(const float_v &v)
Definition: PndFTSCATrackParamVector.h:1404
r1
double r1
Definition: RiemannTrackFinderTest.C:14
Zero
static const fvec Zero
Definition: KFParticleBaseSIMD.cxx:23
dz
Double_t dz
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:24
PndFTSCATrackParamVector::GetS
float_v GetS(const float_v &x, const float_v &y, const float_v &Bz) const
Definition: PndFTSCATrackParamVector.cxx:720
PndFTSCATrackParamVector::GetKappa
float_v GetKappa(const float_v &Bz) const
Definition: PndFTSCATrackParamVector.h:1350
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
fZ
Double_t fZ
Definition: PndCaloDraw.cxx:34
PndFTSCATrackParamVector::SetZ
void SetZ(const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1403
v
__m128 v
Definition: P4_F32vec4.h:4
PndFTSCATrackParamVector::Rotate
float_m Rotate(const float_v &alpha, PndFTSCATrackLinearisationVector &t0, const float maxSinPhi=.999f, const float_m &mask=float_m(true))
Definition: PndFTSCATrackParamVector.cxx:1116
PndFTSCATrackParamVector::ApproximateBetheBloch
static float_v ApproximateBetheBloch(const float_v &beta2)
Definition: PndFTSCATrackParamVector.cxx:1019
PndFTSCATrackParamVector::SetSinPhi
void SetSinPhi(const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1405
PndFTSCATrackParamVector::AddTarget
float_m AddTarget(const FTSCATarget &target, const float_m &mask=float_m(true))
Definition: PndFTSCATrackParamVector.cxx:1504
std::basic_ostream
Definition: PndFTSCATrackParamVector.h:1258
h1
TH1F * h1
Definition: anal_hit_digi_cluster_fwendcap.C:121
t1
Int_t t1
Definition: hist-t7.C:106
init
static int init
Definition: ranlxd.cxx:374
dE
Double_t dE
Definition: anasim.C:58
rcp
T rcp(T val)
Definition: PndFTSCADef.h:52
C
int Pic_FED Eff_lEE C()
Definition: Pic_FEP-Eff_lE2.C:1
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
PndFTSCATrackParamVector::FilterVtx
float_m FilterVtx(const float_v &xV, const float_v &yV, const CAX1X2MeasurementInfo &info, float_v &extrDx, float_v &extrDy, float_v J[], const float_m &active=float_m(true))
Definition: PndFTSCATrackParamVector.h:1917
PndFTSCATrackParamVector::GetDistXZ2
float_v GetDistXZ2(const PndFTSCATrackParamVector &t) const
Definition: PndFTSCATrackParamVector.cxx:710
CAMath::RSqrt
static T RSqrt(const T &x)
Definition: PndCAMath.h:38
a
Int_t a
Definition: anaLmdDigi.C:126
PndFTSCATrackParamVector::SetCov
void SetCov(int i, const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1396
PndFTSCATrackParamVector::SetChi2
void SetChi2(const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1414
PndFTSCATrackParamVector::Transport
float_m Transport(const int_v &ista, const PndFTSCAParam &param, const float_m &mask=float_m(true))
Definition: PndFTSCATrackParamVector.cxx:1440
t0
Int_t t0
Definition: hist-t7.C:106
PndFTSCATrackParamVector::InitDirection
void InitDirection(float_v r0, float_v r1, float_v r2)
Definition: PndFTSCATrackParamVector.h:1290
std::ostream
basic_ostream< char, char_traits< char > > ostream
Definition: PndFTSCATrackParamVector.h:1258
c6
TCanvas * c6
Definition: macro/detectors/emc/dedicated/reco_analys.C:140
PndFTSCATrackParamVector::Cov
const float_v * Cov() const
Definition: PndFTSCATrackParamVector.h:1366
PndFTSCATrackParamVector::TransportJ0ToX0
float_m TransportJ0ToX0(const float_v &x0, const float_v &cBz, float_v &extrDy, float_v &extrDz, float_v J[], const float_m &active=float_m(true))
Definition: PndFTSCATrackParamVector.h:1994
F
TFile * F
Definition: anaLmdReco.C:31
y0
Double_t y0
Definition: checkhelixhit.C:71
f
TFile * f
Definition: bump_analys.C:12
L1FieldRegion::cx0
float_v cx0
Definition: L1Field.h:132
PndFTSCATrackParamVector::SetY
void SetY(const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1402
PndFTSCATrackParamVector::SetQPt
void SetQPt(const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1410
PndFTSCATrackParamVector::CalculateFitParameters
void CalculateFitParameters(PndFTSCATrackFitParam &par, const float_v &mass=0.13957f)
Definition: PndFTSCATrackParamVector.cxx:1041
z
Double_t z
Definition: createRootFscGeometryFile.C:146
PndFTSCATrackParamVector::InitCovMatrix
void InitCovMatrix(float_v d2QMom=0.f)
Definition: PndFTSCATrackParamVector.cxx:596
PndFTSCATrackParamVector::operator>>
friend std::istream & operator>>(std::istream &, PndFTSCATrackParamVector &)
Definition: PndFTSCATrackParamVector.cxx:1522
out
TFile * out
Definition: reco_muo.C:20
std::istream
basic_istream< char, char_traits< char > > istream
Definition: PndFTSCATrackParamVector.h:1256
c1
c1
Definition: plot_dirc.C:35
SetField
fRun SetField(fField)
c3
c3
Definition: plot_dirc.C:50
PndFTSCATrackParamVector::BetheBlochSolid
static float_v BetheBlochSolid(const float_v &bg)
Definition: PndFTSCATrackParamVector.cxx:986
PndFTSCATrackParamVector::TransportToX0WithMaterial
float_m TransportToX0WithMaterial(const float_v &x, const float_v &XOverX0, const float_v &XThimesRho, const float_v &Bz, const float maxSinPhi=.999f)
Definition: PndFTSCATrackParamVector.cxx:929
CAX1X2MeasurementInfo::C10
const float_v & C10() const
Definition: CAX1X2MeasurementInfo.h:17
p2
TPad * p2
Definition: hist-t7.C:117
dx
double dx
Definition: runGeaneTestSimBox.cxx:44
PndFTSCATrackParamVector::FilterWithMaterial
float_m FilterWithMaterial(const float_v &y, const float_v &z, float_v err2Y, float_v errYZ, float_v err2Z, float maxSinPhi=0.999f, const float_m &mask=float_m(true), const int_v &hitNDF=int_v(2), const float_v &chi2Cut=10e10f)
Definition: PndFTSCATrackParamVector.cxx:1168
c22
TCanvas * c22
Definition: plot_acceptance_stt.C:253
c20
TCanvas * c20
Definition: plot_acceptance_stt.C:241
L1FieldRegion::cy0
float_v cy0
Definition: L1Field.h:133
fY
Double_t fY
Definition: PndCaloDraw.cxx:34
PndFTSCATrackParamVector::RotateXY
void RotateXY(float_v alpha, float_v &x, float_v &y, float_v &sin, const float_m &mask=float_m(true)) const
Definition: PndFTSCATrackParamVector.h:1903
PndFTSCATrackParamVector::SetQMomentum
void SetQMomentum(const float_v &v)
Definition: PndFTSCATrackParamVector.h:1409
operator>>
TBuffer & operator>>(TBuffer &buf, PndAnaPidSelector *&obj)
Definition: PndAnaPidSelector.cxx:34
h2
TH1F * h2
Definition: anal_hit_digi_cluster_fwendcap.C:122
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndFTSCATrackParamVector::GetDCAPoint
void GetDCAPoint(const float_v &x, const float_v &y, const float_v &z, float_v *px, float_v *py, float_v *pz, const float_v &Bz) const
Definition: PndFTSCATrackParamVector.cxx:737
PndFTSCATrackParamVector::TransportToX0
float_m TransportToX0(const float_v &x, const float_v &Bz, const float maxSinPhi=.999f, const float_m &mask=float_m(true))
Definition: PndFTSCATrackParamVector.cxx:880
L1FieldRegion::cy2
float_v cy2
Definition: L1Field.h:133
L1FieldRegion::cz0
float_v cz0
Definition: L1Field.h:134
PndFTSCATrackParamVector::SetSignCosPhi
void SetSignCosPhi(const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1412
std::char_traits
Definition: PndFTSCATrackParamVector.h:1255
PndFTSCATrackParamVector::BetheBlochGeant
static float_v BetheBlochGeant(const float_v &bg, const float_v &kp0=2.33f, const float_v &kp1=0.20f, const float_v &kp2=3.00f, const float_v &kp3=173e-9f, const float_v &kp4=0.49848f)
Definition: PndFTSCATrackParamVector.cxx:945
PndFTSCATrackParamVector::SetPar
void SetPar(int i, const float_v &v)
Definition: PndFTSCATrackParamVector.h:1393
PndFTSCATrackParamVector::Cov
const float_v & Cov(int i) const
Definition: PndFTSCATrackParamVector.h:1361
t
TTree * t
Definition: bump_analys.C:13
PndFTSCATrackParamVector::Par
const float_v & Par(int i) const
Definition: PndFTSCATrackParamVector.h:1360
PndFTSCATrackParamVector::InitByTarget
void InitByTarget(const FTSCATarget &target)
Definition: PndFTSCATrackParamVector.cxx:619
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
PndFTSCATrackParamVector::SetDzDs
void SetDzDs(const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1407
y
Double_t y
Definition: createRootFscGeometryFile.C:159
alpha
double alpha
Definition: f_Init.h:9
L1FieldRegion::cx2
float_v cx2
Definition: L1Field.h:132
L1FieldRegion::z0
float_v z0
Definition: L1Field.h:135
PndFTSCATrackParamVector::SetTrackParam
void SetTrackParam(const PndFTSCATrackParamVector &param, const float_m &m=float_m(true))
Definition: PndFTSCATrackParamVector.h:1371
PndFTSCATrackParamVector::operator<<
friend std::ostream & operator<<(std::ostream &, const PndFTSCATrackParamVector &)
Definition: PndFTSCATrackParamVector.cxx:1543
TrackParamVector
PndCATrackParamVector TrackParamVector
Definition: PndCATrackParamVector.h:558
PndFTSCATrackParamVector::SetX
void SetX(const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1401
c9
TCanvas * c9
Definition: anal_hit_digi_cluster_fwendcap.C:648
PndFTSCATrackParamVector::GetDist2
float_v GetDist2(const PndFTSCATrackParamVector &t) const
Definition: PndFTSCATrackParamVector.cxx:700
CAX1X2MeasurementInfo::C00
const float_v & C00() const
Definition: CAX1X2MeasurementInfo.h:16
One
static const fvec One
Definition: KFParticleBaseSIMD.cxx:24
PndFTSCATrackParamVector::ConvertTrackParamToVector
void ConvertTrackParamToVector(PndFTSCATrackParam t0[float_v::Size], int nTracksV)
Definition: PndFTSCATrackParamVector.cxx:557
debugKF
PndFTSCANoDebugStream & debugKF()
Definition: debug.h:74
r2
double r2
Definition: RiemannTrackFinderTest.C:19
pz
double pz[39]
Definition: pipisigmas.h:14
PndFTSCATrackParamVector::Filter
float_m Filter(const FTSCAHitV &hit, const PndFTSCAParam &param, const float_m &mask=float_m(true), const float_v &chi2Cut=10e10f)
Definition: PndFTSCATrackParamVector.cxx:1464
PndFTSCATrackParamVector::SetPar
void SetPar(int i, const float_v &v, const float_m &m)
Definition: PndFTSCATrackParamVector.h:1394
Generated on Wed Apr 3 2019 10:02:51 for FairRoot/PandaRoot by   doxygen 1.8.5