FairRoot/PandaRoot
PndFTSCATrackParamVector.cxx
Go to the documentation of this file.
1 // $Id: PndFTSCATrackParamVector.cxx,v 1.1.1.1 2010/07/26 20:55:38 ikulakov Exp $
2 // **************************************************************************
3 // This file is property of and copyright by the ALICE HLT Project *
4 // ALICE Experiment at CERN, All rights reserved. *
5 // *
6 // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
7 // Ivan Kisel <kisel@kip.uni-heidelberg.de> *
8 // for The ALICE HLT Project. *
9 // *
10 // Developed by: Igor Kulakov <I.Kulakov@gsi.de> *
11 // Mykhailo Pugach <M.Pugach@gsi.de> *
12 // Maksym Zyzak <M.Zyzak@gsi.de> *
13 // *
14 // Permission to use, copy, modify and distribute this software and its *
15 // documentation strictly for non-commercial purposes is hereby granted *
16 // without fee, provided that the above copyright notice appears in all *
17 // copies and that both the copyright notice and this permission notice *
18 // appear in the supporting documentation. The authors make no claims *
19 // about the suitability of this software for any purpose. It is *
20 // provided "as is" without express or implied warranty. *
21 // *
22 //***************************************************************************
23 
24 
25 #ifdef PANDA_FTS
26 
27 
28 
29 #include "PndFTSCATrackParamVector.h"
30 #include "PndFTSCAMath.h"
31 #include "PndFTSCAParam.h"
32 #include "FTSCAHitsV.h"
33 #include "FTSCATarget.h"
34 
35 #include <iostream>
36 #include <iomanip>
37 //#ifndef NVALGRIND
38 //#include <valgrind/memcheck.h>
39 //#endif
40 #include <assert.h>
41 #include "debug.h"
42 #ifdef DRAW
43 #include "PndFTSCADisplay.h"
44 #endif
45 void PndFTSCATrackParamVector::ConvertTrackParamToVector( PndFTSCATrackParam t0[float_v::Size],
46  int nTracksV )
47 {
48  float_v tmpVec;
49  int_v tmpVecShort;
50  float_v::Memory tmpFloat;
51  int_v::Memory tmpShort;
52 
53  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Z();
54  tmpVec.load( tmpFloat );
55  SetZ(tmpVec);
56 
57  for(int iP=0; iP<5; iP++)
58  {
59  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Par()[iP];
60  tmpVec.load( tmpFloat );
61  SetPar(iP,tmpVec);
62  }
63  for(int iC=0; iC<15; iC++)
64  {
65  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Cov()[iC];
66  tmpVec.load( tmpFloat );
67  SetCov(iC,tmpVec);
68  }
69  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Chi2();
70  tmpVec.load( tmpFloat );
71  SetChi2(tmpVec);
72  for(int iV=0; iV < nTracksV; iV++) tmpShort[iV] = t0[iV].NDF();
73  tmpVecShort.load( tmpShort );
74  SetNDF(tmpVecShort);
75 
76  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Angle();
77  tmpVec.load( tmpFloat );
78  SetAngle(tmpVec);
79 }
80 
81 
82 void PndFTSCATrackParamVector::InitCovMatrix( float_v /*d2QMom*/ ) //[R.K. 9/2018] unused
83 {
84  SetCov( 0.f, 1.f );
85  SetCov( 1, 0.f );
86  SetCov( 2, 10.f );
87  SetCov( 3, 0.f );
88  SetCov( 4, 0.f );
89  SetCov( 5, 1.f );
90  SetCov( 6, 0.f );
91  SetCov( 7, 0.f );
92  SetCov( 8, 0.f );
93  SetCov( 9, 10.f );
94  SetCov(10, 0.f );
95  SetCov(11, 0.f );
96  SetCov(12, 0.f );
97  SetCov(13, 0.f );
98  SetCov(14, 10.f /*d2QMom*/ );
99 
100  //SetChi2( 0.f );
101 
102  //SetNDF( -5 );
103 }
104 
105 void PndFTSCATrackParamVector::InitByTarget( const FTSCATarget& t )
106 {
107  SetX( t.X1() );
108  SetY( t.X2() );
109  SetZ( t.X0() );
110  SetTx( 0.f );
111  SetTy( 0.f );
112  SetQP( 0.f );
113  SetCov( 0, 1.f/*t.Err2X1()*/ );
114  SetCov( 1, 0.f );
115  SetCov( 2, 10.f /*t.Err2X2()*/ );
116  SetCov( 3, 0.f );
117  SetCov( 4, 0.f );
118  SetCov( 5, 1.f );
119  SetCov( 6, 0.f );
120  SetCov( 7, 0.f );
121  SetCov( 8, 0.f );
122  SetCov( 9, 10.f );
123  SetCov(10, 0.f );
124  SetCov(11, 0.f );
125  SetCov(12, 0.f );
126  SetCov(13, 0.f );
127  SetCov(14, 10.f /*t.Err2QMom()*/ );
128 
129  SetChi2( 0.f );
130 
131  SetNDF( -5 + t.NDF() );
132 
133  SetField( 0, t.B( 0 ), t.ZB( 0 ) );
134  SetField( 1, t.B( 1 ), t.ZB( 1 ) );
135 }
136 
137 void PndFTSCATrackParamVector::InitByHit( const FTSCAHitV& hit, const PndFTSCAParam& /*param*/, const float_v& dQP ) //[R.K. 9/2018] unused
138 {
139  SetX( hit.X1() );
140  SetY( hit.X2() );
141  SetZ( hit.X0() );
142  SetTx( 0.f );
143  SetTy( 0.f );
144  SetQP( 0.f );
145 
146  //const int ista0 = hit.IStation();
147  //const FTSCAStation sta0 = param.Station(ista0);
148  //const L1XYMeasurementInfo &info = sta0.XYInfo;
149 
150  SetCov( 0, hit.Err2X1() );
151  //SetCov( 0, info.C00 );
152 
153  SetCov( 1, 0.f );
154  SetCov( 2, hit.Err2X2() );
155  SetCov( 3, 0.f );
156  SetCov( 4, 0.f );
157  SetCov( 5, 1.f );
158  SetCov( 6, 0.f );
159  SetCov( 7, 0.f );
160  SetCov( 8, 0.f );
161  SetCov( 9, 1.f );
162  SetCov(10, 0.f );
163  SetCov(11, 0.f );
164  SetCov(12, 0.f );
165  SetCov(13, 0.f );
166  SetCov(14, dQP*dQP );
167 
168  SetChi2( 0.f );
169 
170  SetNDF( -4 );
171 }
172 // void PndFTSCATrackParamVector::InitByHit( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_v& dQP )
173 // {
174 // SetX( hit.X1() );
175 // SetY( hit.X2() );
176 // SetZ( hit.X0() );
177 // SetTx( 0.f );
178 // SetTy( 0.f );
179 // SetQP( 0.f );
180 
181 // const int ista0 = hit.IStation();
182 // const FTSCAStation sta0 = param.Station(ista0);
183 // const L1XYMeasurementInfo &info = sta0.XYInfo;
184 
185 // SetCov( 0, info.C00 );
186 // SetCov( 1, info.C10 );
187 // SetCov( 2, info.C11 );
188 // SetCov( 3, 0.f );
189 // SetCov( 4, 0.f );
190 // SetCov( 5, 1.f );
191 // SetCov( 6, 0.f );
192 // SetCov( 7, 0.f );
193 // SetCov( 8, 0.f );
194 // SetCov( 9, 1.f );
195 // SetCov(10, 0.f );
196 // SetCov(11, 0.f );
197 // SetCov(12, 0.f );
198 // SetCov(13, 0.f );
199 // SetCov(14, dQP*dQP );
200 
201 // SetChi2( 0.f );
202 
203 // SetNDF( -3 );
204 
205 // L1FieldValue B;
206 // sta0.fieldSlice.GetFieldValue( hit.X1(), hit.X2(), B );
207 // // SetField( 1, B, sta0.z );
208 
209 // if ( ista0 > 0 ) {
210 // const FTSCAStation staC = param.Station(ista0/2); // The field will be used to propaget to the target. Therefore take station in a middle between target and hit.
211 // const float_v tx1 = hit.X1()/hit.X0(); // suppose target is at 0
212 // const float_v tx2 = hit.X2()/hit.X0(); // suppose target is at 0
213 // staC.fieldSlice.GetFieldValue( tx1*staC.z, tx2*staC.z, B );
214 // // SetField( 0.f, B, staC.z );
215 // }
216 // }
217 /*
218 float_m PndFTSCATrackParamVector::Transport( const FTSCAHitV& hit, const L1FieldRegion& F, const PndFTSCAParam& param, const float_m &mask )
219 {
220  float_m active = mask & hit.IsValid();
221  const char ista = hit.IStation();
222  active &= TransportToX0WithMaterial( hit.X0(), F, param.Station(ista).materialInfo, QP(), active );
223  return active;
224 }
225 
226 float_m PndFTSCATrackParamVector::Transport( const int_v& ista, const int_v& iVSta, const PndFTSCAParam& param, const float_m &mask )
227 {
228  float_m active = mask;
229  const float_v &DZFromStation = ( param.GetNVirtualStations(ista, active) - iVSta ) * param.Station(ista[0]).dZVirtualStation;
230  const float_v &z2 = param.GetX0( ista, active ) - DZFromStation;
231  const float_v &dz = z2 - fZ;
232 
233  L1FieldRegion f;
234  {
235  const CAFieldValue &b2 = param.GetFieldValue( ista, iVSta, fX + fTx*dz, fY + fTy*dz, active );
236  f.Set( b2, z2, fB[1], fZB[1], fB[0], fZB[0] );
237  fB[0].UpdateValue( fB[1], mask );
238  fB[1].UpdateValue( b2, mask );
239  fZB[0](mask) = fZB[1];
240  fZB[1](mask) = z2;
241  }
242  active &= TransportToX0WithMaterial( z2, f, L1MaterialInfo(), QP(), active );
243  return active;
244 }
245 
246 float_m PndFTSCATrackParamVector::Transport( const int_v& ista, const PndFTSCAParam& param, const float_m &mask )
247 {
248  const float_m toTransport = static_cast<float_m>(ista < param.NStations()) & mask;
249  const float_m doNotTransport = mask & !toTransport;
250 
251  float_m active = toTransport;
252 // active &= Rotate( -fAlpha + hit.Angle(), tR, .999f, active );
253  const float_v &zSta = param.GetX0( ista, active );
254  const float_v &dz = zSta - fZ;
255 
256  L1FieldRegion f;
257  {
258  const float_v &z2 = zSta;
259  const CAFieldValue &b2 = param.GetFieldValue( ista, fX + fTx*dz, fY + fTy*dz, active );
260  if ( CAMath::Abs(fZB[0][0] - 147.481f) < 0.01f ) { // temporary, check whether we at first station, then parabola aproximation would be inprecise because dz01 >> dz12
261  f.Set( b2, z2, fB[1], fZB[1] );
262  } else {
263  f.Set( b2, z2, fB[1], fZB[1], fB[0], fZB[0] );
264  }
265  }
266 
267  float_v qp0 = QP();
268  active &= TransportToX0WithMaterial( param.GetX0( ista, active ), f, param.Station(ista[0]).materialInfo, qp0, active );
269  return active || doNotTransport;
270 }*/
271 
272 //THIS INTERFACE WE USE FOR TRANSPORT
273 float_m PndFTSCATrackParamVector::TransportByLine( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_m &mask )
274 {
275  float_m active = mask & hit.IsValid();
276  active &= TransportToX0Line(hit.X0(), mask);
277  active &= PassMaterial( param.Station(0).materialInfo, fQP, active );
278  const float_v mass2 = 0.13957f*0.13957f;
279  float_v direction = -1.f;
280  if(fDirection)
281  direction = 1.f;
282  EnergyLossCorrection(mass2, param.Station(0).materialInfo.RadThick, fQP, direction, active);
283  return active;
284 }
285 
286 float_m PndFTSCATrackParamVector::Transport( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_m &mask )
287 {
288  float_v qp0 = QP();
289  return Transport(hit, param, qp0, mask);
290 }
291 
292 
293 float_m PndFTSCATrackParamVector::Transport( const FTSCAHit& hit, const PndFTSCAParam& param, const float_m &mask )
294 {
295  float_v qp0 = QP();
296  //BUG const float_m bufm = static_cast<float_m> (uint_v::IndexesFromZero() < 1);
297  const float_m bufm1 = static_cast<float_m> ( (fX != 0) && (fZ != 0) );
298  //cout<<"bufm "<<bufm<<endl;
299  //cout<<"bufm1 "<<bufm1<<endl;
300  //BUG return Transport(FTSCAHitV(&hit,bufm), param, qp0, mask);
301  return Transport(FTSCAHitV(&hit,bufm1,true), param, qp0, mask);
302 }
303 
304 
305 //THIS SUB-INTERFACE IS CALLED FROM ALL OTHER TRANSPORT-INTERFACES
306 float_m PndFTSCATrackParamVector::Transport( const FTSCAHitV& hit, const PndFTSCAParam& param, float_v& qp0, const float_m &mask )
307 {
308 
309  float_m active = mask & hit.IsValid();
310 // active &= Rotate( -fAlpha + hit.Angle(), tR, .999f, active );
311  int_v ista(Vc::Zero);
312  ista(active) = hit.IStations();
313 
314  // int_v NVS = param.GetNVirtualStations( ista, active );
315 // NVS(!static_cast<int_m>(active)) = 0;
316 // for( short int iVS = 0; iVS < NVS.max(); iVS++ ) {
317 // const float_m &valid = active && static_cast<float_m>(iVS < NVS);
318 // Transport( ista, iVS, param, valid );
319 // // cout << ista << valid << Z() << endl;
320 // }
321 
322  if(!( ( (ista==32 || ista==16) && int_m(active)).isEmpty() ) && fDirection) //z=606.995; forward direction
323  {
324  int_m maskVS = ((ista == 32) || (ista==16)) && active;
325  float_m activeVS = float_m(maskVS);
326 
327  int_v N(1);
328  N(maskVS) = param.GetNVirtualStations( ista, activeVS);
329 
330  int_v iVS = 0;
331 
332  float_v z0 = fZ;
333  //float_v zVirtualStation(1.f);
334  while( !((iVS<N).isEmpty()) )
335  {
336  float_v zVirtualStation = z0 + (hit.X0() - z0)/float_v(N+1)*float_v(iVS+1);
337 
338  L1FieldRegion virtualF;
339  int_v nVS = N - iVS - int_v(1);
340  UpdateFieldValues(hit, nVS, zVirtualStation, param, virtualF, maskVS);
341  activeVS &= TransportToX0( zVirtualStation, virtualF, qp0, activeVS );
342 // active &= float_m(!maskVS) || activeVS;
343 #ifdef DRAW
344  foreach_bit( unsigned int iV, active )
345 // int iV=0;
346  {
347  if(activeVS[iV])
348  PndFTSCADisplay::Instance().DrawGBPoint( X()[iV], Y()[iV], Z()[iV], kYellow-2, 1.25 );
349  }
350 #endif
351  iVS += int_v(1);
352  }
353  }
354  if(!( ( (ista==31 || ista==15) && int_m(active)).isEmpty() ) && !fDirection) //backward direction
355  {
356  int_m maskVS = ((ista == 31) || (ista==15)) && active;
357  float_m activeVS = float_m(maskVS);
358  int_v N(1);
359  N(maskVS) = param.GetNVirtualStations( ista+1, activeVS);
360 
361  int_v iVS = N - 1;
362 
363  float_v z0 = fZ;
364 
365  float_v zVirtualStation(1.f);
366  while( !((iVS>-1).isEmpty()) )
367  {
368  zVirtualStation(maskVS) = z0 - (z0 -hit.X0())/float_v(N+1)*float_v(N - iVS);
369  L1FieldRegion virtualF;
370  int_v nVS = N - iVS - 1;
371  UpdateFieldValues(hit, nVS, zVirtualStation, param, virtualF, maskVS);
372  activeVS &= TransportToX0( zVirtualStation, virtualF, qp0, activeVS );
373 // active &= float_m(!maskVS) || activeVS;
374 #ifdef DRAW
375  foreach_bit( unsigned int iV, active )
376 // int iV=0;
377  {
378  if(activeVS[iV])
379  PndFTSCADisplay::Instance().DrawGBPoint( X()[iV], Y()[iV], Z()[iV], kMagenta-2, 1.25 );
380  }
381 // std::cout << Z()[0] << " " << zVirtualStation << " " << std::endl;
382 #endif
383  iVS -= int_v(1);
384  }
385  }
386 
387  L1FieldRegion f;
388  UpdateFieldValues(hit, param, f, mask);
389 
390  active &= TransportToX0WithMaterial( hit.X0(), f, param.Station(0).materialInfo, qp0, active ); // suppose material is same for all stations
391 
392  return active;
393 }
394 
395 void PndFTSCATrackParamVector::UpdateFieldValues(const FTSCAHitV& hit, const PndFTSCAParam& param, L1FieldRegion& f, const float_m& mask)
396 {
397  float_m active = mask & hit.IsValid();
398  const float_v &z2 = param.GetX0( hit.IStations(), active );
399  const CAFieldValue &b2 = param.GetFieldValue( hit.IStations(), fX, fY, active );
400  if ( CAMath::Abs(fZB[0][0] - 147.481f) < 0.01f ) {
401  f.Set( b2, z2, fB[1], fZB[1] );
402  } else {
403  f.Set( b2, z2, fB[1], fZB[1], fB[0], fZB[0] );
404  }
405  fB[0].UpdateValue( fB[1], mask );
406  fB[1].UpdateValue( b2, mask );
407  fZB[0](mask) = fZB[1];
408  fZB[1](mask) = z2;
409 }
410 
411 void PndFTSCATrackParamVector::UpdateFieldValues(const FTSCAHitV& hit, int_v& iVrt, float_v& zVirtualStation, const PndFTSCAParam& param, L1FieldRegion& f, const float_m& mask)
412 {
413  float_m active = mask & hit.IsValid();
414  int_v iStation = hit.IStations();
415  if(!fDirection)
416  iStation += 1;
417 
418  const CAFieldValue &b2 = param.GetFieldValue( iStation, iVrt, fX, fY, active );
419 
420  f.Set( b2, zVirtualStation, fB[1], fZB[1], fB[0], fZB[0] );
421 
422  fB[0].UpdateValue( fB[1], mask );
423  fB[1].UpdateValue( b2, mask );
424  fZB[0](mask) = fZB[1];
425  fZB[1](mask) = zVirtualStation;
426 }
427 
428 //THIS INTERFACE WE USE FOR FILTER
429 float_m PndFTSCATrackParamVector::Filter( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_m &mask, const float_v& chi2Cut )
430 {
431  float_m active = mask & hit.IsValid();
432 
433  FTSCAStripInfoVector stripInfo;
434 
435  param.GetStripInfo(stripInfo, hit.IStations(), mask);
436 
437  float_v Chi2Cut;
438  const float_m chi2mask = (fZ<311.f);
439  Chi2Cut(chi2mask) = 10.8f;
440  Chi2Cut(!chi2mask) = 13.8f;
441 
442  //to accept hits on the long distance within the mag-field
443  const float_m accept_longDistance_hit1 = (hit.IStations()>=24 && hit.IStations()<31);
444  Chi2Cut(accept_longDistance_hit1) *= 100.f;
445 
446  const float_m accept_longDistance_hit2 = (hit.IStations()>=31);
447  Chi2Cut(accept_longDistance_hit2) *= 300.f;
448 
449  active &= Filter( hit.XWire1(), hit.XWire2(), hit.RSigned(), stripInfo, hit.Err2R(), active, chi2Cut );
450  return active;
451 }
452 
453 float_m PndFTSCATrackParamVector::Filter( const FTSCAHit& hit, const PndFTSCAParam& param, const float_m &mask, const float_v& chi2Cut )
454 {
455  float_m active = mask;
456 
457  FTSCAStripInfoVector stripInfo;
458  param.GetStripInfo(stripInfo, int(hit.IStation()), active);
459 
460  active &= Filter( hit.XWire1(), hit.XWire2(), hit.RSigned(), stripInfo, hit.Err2R(), active, chi2Cut );
461  return active;
462 }
463 
464 
465 float_m PndFTSCATrackParamVector::AddTarget( const FTSCATarget& target, const float_m &mask )
466 {
467  float_m active = mask;
468  L1FieldRegion fld;
469  if ( fZB[0][0] == 10e10f ) {
470  fld.Set( fB[1], fZB[1], target.B(), target.X0() );
471 // SetField( 0, target.B(), target.X0() );
472  }
473  else {
474  fld.Set( fB[1], fZB[1], fB[0], fZB[0], target.B(), target.X0() );
475  }
476 
477  float_v eX, eY, J04, J14;
478  float_v dz = target.X0() - X0();
479  active &= TransportJXY0ToX0( target.X0(), fld, eX, eY, J04, J14, active );
480  float_v J[6];
481  J[0]= dz; J[1] = 0; J[2]= J04;
482  J[3] = 0; J[4]= dz; J[5]= J14;
483  active &= FilterVtx( target.X1(), target.X2(), L1XYMeasurementInfo( target.Err2X1(), 0.f, target.Err2X2() ), eX, eY, J, active );
484  fNDF (static_cast<int_m>(active)) += target.NDF();
485 
486  return active;
487 }
488 
489 #include <iostream>
490 
491 std::istream &operator>>( std::istream &in, PndFTSCATrackParamVector &t )
492 { // TODO
493  UNUSED_PARAM1(t);
494  // float_v::Memory x, s, p[5], c[15], fChi2;
495  // int_v::Memory fNDF;
496  // for ( int j = 0; j < uint_v::Size; ++j ) {
497  // in >> x[j];
498  // in >> s[j];
499  // for ( int i = 0; i < 5; i++ ) in >> p[i][j];
500  // for ( int i = 0; i < 15; i++ ) in >> c[i][j];
501  // in >> fChi2[j];
502  // in >> fNDF[j];
503  // }
504  // t.fX.load( x );
505  // t.fSignCosPhi.load( s );
506  // for ( int i = 0; i < 5; i++ ) t.fP[i].load( p[i] );
507  // for ( int i = 0; i < 5; i++ ) t.fC[i].load( c[i] );
508  // t.Chi2.load( fChi2 );
509  // t.NDF.load( fNDF );
510  return in;
511 }
512 
513 std::ostream &operator<<( std::ostream &out, const PndFTSCATrackParamVector &t )
514 {
515  out << "PndFTSCATrackParamVector " << std::endl;
516  for ( int i = 0; i < 6; i++ ) out << t.Par(i) << " ";
517  out << std::endl;
518  /*out << "CovMat " << std::endl;
519  //for ( int i = 0; i < 15; i++ )
520  out << t.Cov(0)[0] << std::endl;
521  out << t.Cov(1)[0] << " "<< t.Cov(2)[0] << std::endl;
522  out << t.Cov(3)[0] << " "<< t.Cov(4)[0] <<" "<< t.Cov(5)[0] << std::endl;
523  out << t.Cov(6)[0] << " "<< t.Cov(7)[0] <<" "<< t.Cov(8)[0] <<" "<< t.Cov(9)[0] << std::endl;
524  out << t.Cov(10)[0] <<" "<< t.Cov(11)[0] <<" "<< t.Cov(12)[0] <<" "<< t.Cov(13)[0] <<" "<< t.Cov(14)[0];
525  */
526 
527  return out;
528 }
529 
530 #else // PANDA_FTS
531 
532 #include "PndFTSCATrackParamVector.h"
533 #include "PndFTSCAMath.h"
534 #include "PndFTSCATrackLinearisationVector.h"
535 #include "PndFTSCAParam.h"
536 #include "FTSCAHitsV.h"
537 #include "FTSCATarget.h"
538 
539 #include <iostream>
540 #include <iomanip>
541 #ifndef NVALGRIND
542 #include <valgrind/memcheck.h>
543 #endif
544 #include <assert.h>
545 #include "debug.h"
546 
547 //
548 // Circle in XY:
549 //
550 // kCLight = 0.000299792458;
551 // Kappa = Bz*kCLight*QPt;
552 // R = 1/CAMath::Abs(Kappa);
553 // Xc = X - sin(Phi)/Kappa;
554 // Yc = Y + cos(Phi)/Kappa;
555 //
556 
557 void PndFTSCATrackParamVector::ConvertTrackParamToVector( PndFTSCATrackParam t0[float_v::Size],
558  int nTracksV )
559 {
560  float_v tmpVec;
561  int_v tmpVecShort;
562  float_v::Memory tmpFloat;
563  int_v::Memory tmpShort;
564 
565  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].X();
566  tmpVec.load( tmpFloat );
567  SetX(tmpVec);
568  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].SignCosPhi();
569  tmpVec.load( tmpFloat );
570  SetSignCosPhi(tmpVec);
571 
572  for(int iP=0; iP<5; iP++)
573  {
574  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Par()[iP];
575  tmpVec.load( tmpFloat );
576  SetPar(iP,tmpVec);
577  }
578  for(int iC=0; iC<15; iC++)
579  {
580  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Cov()[iC];
581  tmpVec.load( tmpFloat );
582  SetCov(iC,tmpVec);
583  }
584  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Chi2();
585  tmpVec.load( tmpFloat );
586  SetChi2(tmpVec);
587  for(int iV=0; iV < nTracksV; iV++) tmpShort[iV] = t0[iV].NDF();
588  tmpVecShort.load( tmpShort );
589  SetNDF(tmpVecShort);
590 
591  for(int iV=0; iV < nTracksV; iV++) tmpFloat[iV] = t0[iV].Angle();
592  tmpVec.load( tmpFloat );
593  SetAngle(tmpVec);
594 }
595 
596 void PndFTSCATrackParamVector::InitCovMatrix( float_v d2QMom )
597 {
598  SetCov( 0.f, 100.f );
599  SetCov( 1, 0.f );
600  SetCov( 2, 100.f );
601  SetCov( 3, 0.f );
602  SetCov( 4, 0.f );
603  SetCov( 5, 1.f );
604  SetCov( 6, 0.f );
605  SetCov( 7, 0.f );
606  SetCov( 8, 0.f );
607  SetCov( 9, 1.f );
608  SetCov(10, 0.f );
609  SetCov(11, 0.f );
610  SetCov(12, 0.f );
611  SetCov(13, 0.f );
612  SetCov(14, d2QMom );
613 
614  SetChi2( 0.f );
615 
616  SetNDF( -5 );
617 }
618 
619 void PndFTSCATrackParamVector::InitByTarget( const FTSCATarget& t )
620 {
621  SetX( t.X0() );
622  SetY( t.X1() );
623  SetZ( t.X2() );
624  SetSinPhi( 1.f );
625  SetSignCosPhi( 0.f );
626  SetDzDs( 0.f );
627  SetQPt( 0.f );
628  SetCov( 0.f, t.Err2X1() );
629  SetCov( 1, 0.f );
630  SetCov( 2, t.Err2X2() );
631  SetCov( 3, 0.f );
632  SetCov( 4, 0.f );
633  SetCov( 5, 1.f );
634  SetCov( 6, 0.f );
635  SetCov( 7, 0.f );
636  SetCov( 8, 0.f );
637  SetCov( 9, 1.f );
638  SetCov(10, 0.f );
639  SetCov(11, 0.f );
640  SetCov(12, 0.f );
641  SetCov(13, 0.f );
642  SetCov(14, t.Err2QMom() );
643 
644  SetChi2( 0.f );
645 
646  SetNDF( -5 + t.NDF() );
647 
648  // SetField( 0, t.B(), t.X0() );
649 }
650 
651 
652 void PndFTSCATrackParamVector::InitByHit( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_v& dQMom )
653 {
654  SetX( hit.X0() );
655  SetY( hit.X1() );
656  SetZ( hit.X2() );
657  SetSinPhi( 0.f );
658  SetSignCosPhi( 1.f );
659  SetDzDs( 0.f );
660  SetQPt( 0.f );
661  SetAngle( hit.Angle() );
662 
663  SetCov( 0, hit.Err2X1() );
664  SetCov( 1, 0.f );
665  SetCov( 2, hit.Err2X2() );
666  SetCov( 3, 0.f );
667  SetCov( 4, 0.f );
668  SetCov( 5, 1.f );
669  SetCov( 6, 0.f );
670  SetCov( 7, 0.f );
671  SetCov( 8, 0.f );
672  SetCov( 9, 1.f );
673  SetCov(10, 0.f );
674  SetCov(11, 0.f );
675  SetCov(12, 0.f );
676  SetCov(13, 0.f );
677  SetCov(14, dQMom*dQMom );
678 
679  SetChi2( 0.f );
680 
681  SetNDF( -3 );
682 
683  UNUSED_PARAM1(param);
684  // const int ista0 = hit.IStation();
685  // const FTSCAStation sta0 = param.Station(ista0);
686 
687  // CAFieldValue B;
688  // sta0.fieldSlice.GetFieldValue( hit.X1(), hit.X2(), B );
689  // SetField( 1, B, sta0.z );
690 
691  // if ( ista0 > 0 ) {
692  // const FTSCAStation staC = param.Station(ista0/2); // The field will be used to propaget to the target. Therefore take station in a middle between target and hit.
693  // const float_v tx1 = hit.X1()/hit.X0(); // suppose target is at 0
694  // const float_v tx2 = hit.X2()/hit.X0(); // suppose target is at 0
695  // staC.fieldSlice.GetFieldValue( tx1*staC.z, tx2*staC.z, B );
696  // SetField( 0, B, staC.z );
697  // }
698 }
699 
700 float_v PndFTSCATrackParamVector::GetDist2( const PndFTSCATrackParamVector &t ) const
701 {
702  // get squared distance between tracks
703 
704  const float_v &dx = GetX() - t.GetX();
705  const float_v &dy = GetY() - t.GetY();
706  const float_v &dz = GetZ() - t.GetZ();
707  return dx*dx + dy*dy + dz*dz;
708 }
709 
710 float_v PndFTSCATrackParamVector::GetDistXZ2( const PndFTSCATrackParamVector &t ) const
711 {
712  // get squared distance between tracks in X&Z
713 
714  const float_v &dx = GetX() - t.GetX();
715  const float_v &dz = GetZ() - t.GetZ();
716  return dx*dx + dz*dz;
717 }
718 
719 
720 float_v PndFTSCATrackParamVector::GetS( const float_v &_x, const float_v &_y, const float_v &Bz ) const
721 {
722  //* Get XY path length to the given point
723 
724  const float_v &k = GetKappa( Bz );
725  const float_v &ex = GetCosPhi();
726  const float_v &ey = GetSinPhi();
727  const float_v &x = _x - GetX();
728  const float_v &y = _y - GetY();
729  float_v dS = x * ex + y * ey;
730  const float_m &mask = CAMath::Abs( k ) > 1.e-4f;
731  if ( !mask.isEmpty() ) {
732  dS( mask ) = CAMath::ATan2( k * dS, 1 + k * ( x * ey - y * ex ) ) / k;
733  }
734  return dS;
735 }
736 
737 void PndFTSCATrackParamVector::GetDCAPoint( const float_v &x, const float_v &y, const float_v &z,
738  float_v *xp, float_v *yp, float_v *zp, const float_v &Bz ) const
739 {
740  //* Get the track point closest to the (x,y,z)
741 
742  const float_v &x0 = GetX();
743  const float_v &y0 = GetY();
744  const float_v &k = GetKappa( Bz );
745  const float_v &ex = GetCosPhi();
746  const float_v &ey = GetSinPhi();
747  const float_v &dx = x - x0;
748  const float_v &dy = y - y0;
749  const float_v &ax = dx * k + ey;
750  const float_v &ay = dy * k - ex;
751  const float_v &a = sqrt( ax * ax + ay * ay );
752  *xp = x0 + ( dx - ey * ( ( dx * dx + dy * dy ) * k - 2.f * ( -dx * ey + dy * ex ) ) / ( a + 1.f ) ) / a;
753  *yp = y0 + ( dy + ex * ( ( dx * dx + dy * dy ) * k - 2.f * ( -dx * ey + dy * ex ) ) / ( a + 1.f ) ) / a;
754  const float_v s = GetS( x, y, Bz );
755  *zp = GetZ() + GetDzDs() * s;
756  const float_v dZ = CAMath::Abs( GetDzDs() * CAMath::TwoPi() / k );
757  const float_m mask = CAMath::Abs( k ) > 1.e-2f && dZ > .1f;
758  ( *zp )( mask ) += CAMath::Round( ( z - *zp ) / dZ ) * dZ;
759 }
760 
761 
762 //*
763 //* Transport routines
764 //*
765 
766 
767 float_m PndFTSCATrackParamVector::TransportToX0( const float_v &x, PndFTSCATrackLinearisationVector &t0, const float_v &Bz, const float maxSinPhi, float_v *DL, const float_m &_mask )
768 {
769  //* Transport the track parameters to X=x, using linearization at t0, and the field value Bz
770  //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
771  //* linearisation of trajectory t0 is also transported to X=x,
772  //* returns 1 if OK
773  //*
774 // std::cout << t0.QPt()[0] << " " << fP[4][0] << std::endl;
775  const float_v ex = t0.CosPhi();
776  const float_v ey = t0.SinPhi();
777  const float_v k = t0.QPt() * Bz;
778  const float_v dx = x - X();
779 
780  float_v ey1 = k * dx + ey;
781 
782  // check for intersection with X=x
783 
784  float_v ex1 = CAMath::Sqrt( 1.f - ey1 * ey1 );
785  ex1( ex < Vc::Zero ) = -ex1;
786 
787  const float_v dx2 = dx * dx;
788  const float_v ss = ey + ey1;
789  const float_v cc = ex + ex1;
790 
791  const float_v cci = 1.f / cc;
792  const float_v exi = 1.f / ex;
793  const float_v ex1i = 1.f / ex1;
794 
795  float_m mask = _mask && CAMath::Abs( ey1 ) <= maxSinPhi && CAMath::Abs( cc ) >= 1.e-4f && CAMath::Abs( ex ) >= 1.e-4f && CAMath::Abs( ex1 ) >= 1.e-4f;
796 
797  const float_v tg = ss * cci; // tan((phi1+phi)/2)
798 
799  const float_v dy = dx * tg;
800  float_v dl = dx * CAMath::Sqrt( 1.f + tg * tg );
801 
802  dl( cc < Vc::Zero ) = -dl;
803  const float_v dSin = CAMath::Max( float_v( -1.f ), CAMath::Min( float_v( Vc::One ), dl * k * 0.5f ) );
804  const float_v dS = ( CAMath::Abs( k ) > 1.e-4f ) ? ( 2 * CAMath::ASin( dSin ) / k ) : dl;
805  const float_v dz = dS * t0.DzDs();
806 
807  if ( DL ) {
808  ( *DL )( mask ) = -dS * CAMath::Sqrt( 1.f + t0.DzDs() * t0.DzDs() );
809  }
810 
811  const float_v d[3] = { fP[2] - t0.SinPhi(), fP[3] - t0.DzDs(), fP[4] - t0.QPt() };
812 
813  //float H0[5] = { 1,0, h2, 0, h4 };
814  //float H1[5] = { 0, 1, 0, dS, 0 };
815  //float H2[5] = { 0, 0, 1, 0, dxBz };
816  //float H3[5] = { 0, 0, 0, 1, 0 };
817  //float H4[5] = { 0, 0, 0, 0, 1 };
818 
819  const float_v h2 = dx * ( 1.f + ey * ey1 + ex * ex1 ) * exi * ex1i * cci;
820  const float_v h4 = dx2 * ( cc + ss * ey1 * ex1i ) * cci * cci * Bz;
821  const float_v dxBz = dx * Bz;
822 
823  ex1( !mask ) = ex;
824  ey1( !mask ) = ey;
825 
826 // std::cout << "QPt " << t0.QPt()[0] << " " << fP[4][0] <<
827 // " CosPhi " << t0.CosPhi()[0] << " " << ex1[0] << std::endl;
828 
829  t0.SetCosPhi( ex1 );
830  t0.SetSinPhi( ey1 );
831 
832  fX( mask ) = X() + dx;
833  fP[0]( mask ) = Y() + dy + h2 * d[0] + h4 * d[2];
834  fP[1]( mask ) = Z() + dz + dS * d[1];
835  fP[2]( mask ) = t0.SinPhi() + d[0] + dxBz * d[2];
836  fP[2]( CAMath::Abs(fP[2]) > maxSinPhi ) = t0.SinPhi();
837 
838  const float_v c00 = fC[0];
839  const float_v c10 = fC[1];
840  const float_v c11 = fC[2];
841  const float_v c20 = fC[3];
842  const float_v c21 = fC[4];
843  const float_v c22 = fC[5];
844  const float_v c30 = fC[6];
845  const float_v c31 = fC[7];
846  const float_v c32 = fC[8];
847  const float_v c33 = fC[9];
848  const float_v c40 = fC[10];
849  const float_v c41 = fC[11];
850  const float_v c42 = fC[12];
851  const float_v c43 = fC[13];
852  const float_v c44 = fC[14];
853 
854  fC[0]( mask ) = ( c00 + h2 * h2 * c22 + h4 * h4 * c44
855  + 2.f * ( h2 * c20 + h4 * c40 + h2 * h4 * c42 ) );
856 
857  fC[1]( mask ) = c10 + h2 * c21 + h4 * c41 + dS * ( c30 + h2 * c32 + h4 * c43 );
858  fC[2]( mask ) = c11 + 2.f * dS * c31 + dS * dS * c33;
859 
860  fC[3]( mask ) = c20 + h2 * c22 + h4 * c42 + dxBz * ( c40 + h2 * c42 + h4 * c44 );
861  fC[4]( mask ) = c21 + dS * c32 + dxBz * ( c41 + dS * c43 );
862  fC[5]( mask ) = c22 + 2.f * dxBz * c42 + dxBz * dxBz * c44;
863 
864  fC[6]( mask ) = c30 + h2 * c32 + h4 * c43;
865  fC[7]( mask ) = c31 + dS * c33;
866  fC[8]( mask ) = c32 + dxBz * c43;
867  fC[9]( mask ) = c33;
868 
869  fC[10]( mask ) = c40 + h2 * c42 + h4 * c44;
870  fC[11]( mask ) = c41 + dS * c43;
871  fC[12]( mask ) = c42 + dxBz * c44;
872  fC[13]( mask ) = c43;
873  fC[14]( mask ) = c44;
874 //std::cout << fC[10] << " " << fC[11] <<" " << h2 << " " << h4 << " " << c44 << " " << mask << " " << Bz << std::endl;
875  debugKF() << mask << "\n" << *this << std::endl;
876  return mask;
877 }
878 
879 
880 float_m PndFTSCATrackParamVector::TransportToX0( const float_v &x, const float_v &Bz,
881  const float maxSinPhi, const float_m &mask )
882 {
883  //* Transport the track parameters to X=x
884 
885 // assert( ( x == 0 && mask ).isEmpty() );
886 
887  PndFTSCATrackLinearisationVector t0( *this );
888 
889  return TransportToX0( x, t0, Bz, maxSinPhi, 0, mask );
890 }
891 
892 
893 
894 float_m PndFTSCATrackParamVector::TransportToX0WithMaterial( const float_v &x, PndFTSCATrackLinearisationVector &t0, PndFTSCATrackFitParam &par, const float_v &XOverX0,
895 const float_v &XThimesRho, const float_v &Bz, const float maxSinPhi, const float_m &mask_ )
896 {
897  //* Transport the track parameters to X=x taking into account material budget
898 // UNUSED_PARAM1(par);
899  // const float kRho = 1.025e-3f;//0.9e-3;
900 // const float kRadLen = 29.532f;//28.94;
901 // const float kRhoOverRadLen = kRho / kRadLen;
902  // const float kRhoOverRadLen = 7.68e-5;
903  float_v dl;
904 
905  float_m mask = mask_ && TransportToX0( x, t0, Bz, maxSinPhi, &dl, mask_ );
906 // if ( !mask.isEmpty() ) {
907 // CorrectForMeanMaterial( dl * kRhoOverRadLen, dl * kRho, par, mask );
908 // }
909 
910  const float_v dzds2 = t0.DzDs()*t0.DzDs();
911  const float_v sphi2 = t0.SinPhi()*t0.SinPhi();
912  float_v koeff = sqrt( (1 + dzds2)*( 1 + sphi2 ) ); // = dr/dx
913 // mask &= CorrectForMeanMaterial( XOverX0*koeff, XThimesRho*koeff, par, mask );
914  CorrectForMeanMaterial( XOverX0*koeff, XThimesRho*koeff, par, mask );
915 
916  return mask;
917 }
918 
919 
920 float_m PndFTSCATrackParamVector::TransportToX0WithMaterial( const float_v &x, PndFTSCATrackFitParam &par, const float_v &XOverX0,
921 const float_v &XThimesRho, const float_v &Bz, const float maxSinPhi )
922 {
923  //* Transport the track parameters to X=x taking into account material budget
924 
925  PndFTSCATrackLinearisationVector t0( *this );
926  return TransportToX0WithMaterial( x, t0, par, XOverX0, XThimesRho, Bz, maxSinPhi );
927 }
928 
929 float_m PndFTSCATrackParamVector::TransportToX0WithMaterial( const float_v &x, const float_v &XOverX0,
930 const float_v &XThimesRho, const float_v &Bz, const float maxSinPhi )
931 {
932  //* Transport the track parameters to X=x taking into account material budget
933 
934  PndFTSCATrackFitParam par;
935  CalculateFitParameters( par );
936  return TransportToX0WithMaterial( x, par, XOverX0, XThimesRho, Bz, maxSinPhi );
937 }
938 
939 
940 //*
941 //* Multiple scattering and energy losses
942 //*
943 
944 
945 float_v PndFTSCATrackParamVector::BetheBlochGeant( const float_v &bg2,
946  const float_v &kp0,
947  const float_v &kp1,
948  const float_v &kp2,
949  const float_v &kp3,
950  const float_v &kp4 )
951 {
952  //
953  // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
954  //
955  // bg2 - (beta*gamma)^2
956  // kp0 - density [g/cm^3]
957  // kp1 - density effect first junction point
958  // kp2 - density effect second junction point
959  // kp3 - mean excitation energy [GeV]
960  // kp4 - mean Z/A
961  //
962  // The default values for the kp* parameters are for silicon.
963  // The returned value is in [GeV/(g/cm^2)].
964  //
965 
966  const float mK = 0.307075e-3f; // [GeV*cm^2/g]
967  const float _2me = 1.022e-3f; // [GeV/c^2]
968  const float_v &rho = kp0;
969  const float_v &x0 = kp1 * 2.303f;
970  const float_v &x1 = kp2 * 2.303f;
971  const float_v &mI = kp3;
972  const float_v &mZA = kp4;
973  const float_v &maxT = _2me * bg2; // neglecting the electron mass
974 
975  //*** Density effect
976  float_v d2( Vc::Zero );
977  const float_v x = 0.5f * CAMath::Log( bg2 );
978  const float_v lhwI = CAMath::Log( 28.816f * 1e-9f * CAMath::Sqrt( rho * mZA ) / mI );
979  d2( x > x1 ) = lhwI + x - 0.5f;
980  const float_v &r = ( x1 - x ) / ( x1 - x0 );
981  d2( x > x0 && x <= x1 ) = lhwI + x - 0.5f + ( 0.5f - lhwI - x0 ) * r * r * r;
982 
983  return mK*mZA*( float_v( Vc::One ) + bg2 ) / bg2*( 0.5f*CAMath::Log( maxT * maxT / ( mI*mI ) ) - bg2 / ( float_v( Vc::One ) + bg2 ) - d2 );
984 }
985 
986 float_v PndFTSCATrackParamVector::BetheBlochSolid( const float_v &bg )
987 {
988  //------------------------------------------------------------------
989  // This is an approximation of the Bethe-Bloch formula,
990  // reasonable for solid materials.
991  // All the parameters are, in fact, for Si.
992  // The returned value is in [GeV]
993  //------------------------------------------------------------------
994 
995  return BetheBlochGeant( bg );
996 }
997 
998 float_v PndFTSCATrackParamVector::BetheBlochGas( const float_v &bg )
999 {
1000  //------------------------------------------------------------------
1001  // This is an approximation of the Bethe-Bloch formula,
1002  // reasonable for gas materials.
1003  // All the parameters are, in fact, for Ne.
1004  // The returned value is in [GeV]
1005  //------------------------------------------------------------------
1006 
1007  const float_v rho = 0.9e-3f;
1008  const float_v x0 = 2.f;
1009  const float_v x1 = 4.f;
1010  const float_v mI = 140.e-9f;
1011  const float_v mZA = 0.49555f;
1012 
1013  return BetheBlochGeant( bg, rho, x0, x1, mI, mZA );
1014 }
1015 
1016 
1017 
1018 
1019 float_v PndFTSCATrackParamVector::ApproximateBetheBloch( const float_v &beta2 )
1020 {
1021  //------------------------------------------------------------------
1022  // This is an approximation of the Bethe-Bloch formula with
1023  // the density effect taken into account at beta*gamma > 3.5
1024  // (the approximation is reasonable only for solid materials)
1025  //------------------------------------------------------------------
1026 
1027  const float_v &beta2_1subBeta2 = beta2 / ( float_v( Vc::One ) - beta2 ); // beta2 * CAMath::Reciprocal( float_v( Vc::One ) - beta2 );
1028  const float_v &_0p000153_beta2 = 0.153e-3f / beta2;
1029  const float log_3p5mul5940 = 9.942227380852058f; // log( 3.5 * 5940 )
1030  const float log_5940 = 8.68946441235669f; // log( 5940 )
1031  const float_v log_beta2_1subBeta2 = CAMath::Log( beta2_1subBeta2 );
1032 
1033  float_v ret = _0p000153_beta2 * ( log_5940 + log_beta2_1subBeta2 - beta2 );
1034  ret( beta2_1subBeta2 > 3.5f*3.5f ) =
1035  _0p000153_beta2 * ( log_3p5mul5940 + 0.5f * log_beta2_1subBeta2 - beta2 );
1036  ret.setZero( beta2 >= float_v( Vc::One ) );
1037  return ret;
1038 }
1039 
1040 
1041 void PndFTSCATrackParamVector::CalculateFitParameters( PndFTSCATrackFitParam &par, const float_v &mass )
1042 {
1043  //*!
1044 
1045  const float_v p2 = ( float_v( Vc::One ) + fP[3] * fP[3] );
1046  const float_v k2 = fP[4] * fP[4];
1047  const float_v mass2 = mass * mass;
1048  const float_v beta2 = p2 / ( p2 + mass2 * k2 );
1049 
1050  float_v pp2 = 10000.f; pp2( k2 > 1.e-8f ) = p2 / k2; // impuls 2
1051 
1052  //par.fBethe = BetheBlochGas( pp2/mass2);
1053  // par.fBethe = ApproximateBetheBloch( pp2 / mass2 ); Why pp2 / mass2?, should be beta2.
1054  par.fBethe = ApproximateBetheBloch( beta2 );
1055  par.fE = CAMath::Sqrt( pp2 + mass2 );
1056  par.fTheta2 = 14.1f * 14.1f / ( beta2 * pp2 * 1.e6f );
1057  par.fEP2 = par.fE / pp2;
1058 
1059  // Approximate energy loss fluctuation (M.Ivanov)
1060 
1061  const float knst = 0.07f; // To be tuned.
1062  par.fSigmadE2 = knst * par.fEP2 * fP[4];
1063  par.fSigmadE2 = par.fSigmadE2 * par.fSigmadE2;
1064 
1065  par.fK22 = ( float_v( Vc::One ) + fP[3] * fP[3] );
1066  par.fK33 = par.fK22 * par.fK22;
1067  par.fK43 = fP[3] * fP[4] * par.fK22;
1068  par.fK44 = fP[3] * fP[3] * fP[4] * fP[4];
1069 }
1070 
1071 
1072 float_m PndFTSCATrackParamVector::CorrectForMeanMaterial( const float_v &xOverX0, const float_v &xTimesRho, const PndFTSCATrackFitParam &par, const float_m &_mask )
1073 {
1074  //------------------------------------------------------------------
1075  // This function corrects the track parameters for the crossed material.
1076  // "xOverX0" - X/X0, the thickness in units of the radiation length.
1077  // "xTimesRho" - is the product length*density (g/cm^2).
1078  //------------------------------------------------------------------
1079 
1080  float_v &fC22 = fC[5];
1081  float_v &fC33 = fC[9];
1082  float_v &fC40 = fC[10];
1083  float_v &fC41 = fC[11];
1084  float_v &fC42 = fC[12];
1085  float_v &fC43 = fC[13];
1086  float_v &fC44 = fC[14];
1087 
1088  //Energy losses************************
1089 
1090  const float_v &dE = par.fBethe * xTimesRho;
1091 // assert( (dE >= float_v(Vc::Zero) || !_mask ).isFull() );
1092  float_m mask = _mask && dE <= 0.3f * par.fE; //30% energy loss is too much!
1093  const float_v &corr = ( float_v( Vc::One ) - par.fEP2 * dE );
1094  mask &= corr >= 0.3f && corr <= 1.3f;
1095 
1096  fP[4]( mask ) *= corr;
1097  fC40 ( mask ) *= corr;
1098  fC41 ( mask ) *= corr;
1099  fC42 ( mask ) *= corr;
1100  fC43 ( mask ) *= corr;
1101  fC44 ( mask ) *= corr * corr;
1102  fC44 ( mask ) += par.fSigmadE2 * CAMath::Abs( dE );
1103 
1104  //Multiple scattering******************
1105 
1106  assert( (xOverX0 >= float_v(Vc::Zero) || !mask ).isFull() );
1107  const float_v &theta2 = par.fTheta2 * xOverX0;
1108  fC22( mask ) += theta2 * par.fK22 * ( float_v( Vc::One ) - fP[2] * fP[2] );
1109  fC33( mask ) += theta2 * par.fK33;
1110  fC43( mask ) += theta2 * par.fK43;
1111  fC44( mask ) += theta2 * par.fK44;
1112 
1113  return mask;
1114 }
1115 
1116 float_m PndFTSCATrackParamVector::Rotate( const float_v &alpha, PndFTSCATrackLinearisationVector &t0,
1117  const float maxSinPhi, const float_m &mask )
1118 {
1119  //* Rotate the coordinate system in XY on the angle alpha
1120  if ( (CAMath::Abs(alpha) < 1e-6f || !mask).isFull() ) return mask;
1121 
1122  const float_v cA = CAMath::Cos( alpha );
1123  const float_v sA = CAMath::Sin( alpha );
1124  const float_v x0 = X(), y0 = Y(), sP = t0.SinPhi(), cP = t0.CosPhi();
1125  const float_v cosPhi = cP * cA + sP * sA;
1126  const float_v sinPhi = -cP * sA + sP * cA;
1127 
1128  float_m ReturnMask(mask);
1129  ReturnMask &= (!( CAMath::Abs( sinPhi ) > maxSinPhi || CAMath::Abs( cosPhi ) < 1.e-4f || CAMath::Abs( cP ) < 1.e-4f ));
1130 
1131  float_v tmp = alpha*0.15915494f;// 1/(2.f*3.1415f);
1132  ReturnMask &= abs(tmp - round(tmp)) < 0.167f; // allow turn by 60 degree only TODO scalar
1133 
1134  //float J[5][5] = { { j0, 0, 0, 0, 0 }, // Y
1135  // { 0, 1, 0, 0, 0 }, // Z
1136  // { 0, 0, j2, 0, 0 }, // SinPhi
1137  // { 0, 0, 0, 1, 0 }, // DzDs
1138  // { 0, 0, 0, 0, 1 } }; // Kappa
1139 
1140  const float_v j0 = cP / cosPhi;
1141  const float_v j2 = cosPhi / cP;
1142  const float_v d = SinPhi() - sP;
1143 
1144  SetX( x0*cA + y0*sA, ReturnMask );
1145  SetY(-x0*sA + y0*cA, ReturnMask );
1146  t0.SetCosPhi( cosPhi );
1147  t0.SetSinPhi( sinPhi );
1148 
1149  SetSinPhi( sinPhi + j2*d, ReturnMask );
1150 
1151  fC[0](ReturnMask) *= j0 * j0;
1152  fC[1](ReturnMask) *= j0;
1153  fC[3](ReturnMask) *= j0;
1154  fC[6](ReturnMask) *= j0;
1155  fC[10](ReturnMask) *= j0;
1156 
1157  fC[3](ReturnMask) *= j2;
1158  fC[4](ReturnMask) *= j2;
1159  fC[5](ReturnMask) *= j2 * j2;
1160  fC[8](ReturnMask) *= j2;
1161  fC[12](ReturnMask) *= j2;
1162 
1163  fAlpha(ReturnMask) += alpha;
1164 
1165  return ReturnMask;
1166 }
1167 
1168 float_m PndFTSCATrackParamVector::FilterWithMaterial( const float_v &y, const float_v &z, float_v err2Y, float_v errYZ, float_v err2Z, float maxSinPhi, const float_m &mask, const int_v& hitNDF, const float_v& chi2Cut )
1169 {
1170  assert( maxSinPhi > 0.f );
1171  //* Add the y,z measurement with the Kalman filter
1172 
1173  const float_v c00 = fC[0];
1174  const float_v c10 = fC[1];
1175  const float_v c11 = fC[2];
1176  const float_v c20 = fC[3];
1177  const float_v c21 = fC[4];
1178 // float c22 = fC[5];
1179  const float_v c30 = fC[6];
1180  const float_v c31 = fC[7];
1181 // float c32 = fC[8];
1182 // float c33 = fC[9];
1183  const float_v c40 = fC[10];
1184  const float_v c41 = fC[11];
1185 // float c42 = fC[12];
1186 // float c43 = fC[13];
1187 // float c44 = fC[14];
1188 
1189 
1190  const float_v
1191  z0 = y - fP[0],
1192  z1 = z - fP[1];
1193 
1194  // foreach_bit( int i, mask ) {
1195  // cout << i << " e = " << err2Y[i] << " " << errYZ[i] << " " << err2Z[i] << " c = " << c00[i] << " " << c10[i] << " " << c11[i] << " z = " << z0[i] << " " << z1[i] << " chi2 = " << fChi2[i]<< endl;
1196  // }
1197 
1198  err2Y += c00;
1199  err2Z += c11;
1200  errYZ += c10;
1201  float_v d = float_v( Vc::One ) / ( err2Y*err2Z - errYZ*errYZ );
1202 
1203 
1204  float_m success = mask;// if ( ISUNLIKELY( err2Y < 1.e-8f ) || ISUNLIKELY( err2Z < 1.e-8f ) ) return 0;
1205  success &= (err2Y > 1.e-8f ) && ( err2Z > 1.e-8f );
1206 
1207  const float_v mS0 = err2Z*d;
1208  const float_v mS1 = -errYZ*d;
1209  const float_v mS2 = err2Y*d;
1210 
1211  // K = CHtS
1212 
1213  const float_v
1214  k00 = c00 * mS0 + c10*mS1, k01 = c00 * mS1 + c10*mS2,
1215  k10 = c10 * mS0 + c11*mS1, k11 = c10 * mS1 + c11*mS2,
1216  k20 = c20 * mS0 + c21*mS1, k21 = c20 * mS1 + c21*mS2,
1217  k30 = c30 * mS0 + c31*mS1, k31 = c30 * mS1 + c31*mS2,
1218  k40 = c40 * mS0 + c41*mS1, k41 = c40 * mS1 + c41*mS2;
1219 
1220  const float_v sinPhi = fP[2] + k20 * z0 + k21 * z1;
1221 
1222  success &= CAMath::Abs( sinPhi ) < maxSinPhi;
1223 
1224  fNDF( static_cast<int_m>(success) ) += hitNDF;
1225  fChi2(success) += mS0 * z0 * z0 + 2 * z0 * z1 * mS1 + mS2 * z1 * z1;
1226  success &= fChi2 < chi2Cut;
1227  if ( success.isEmpty() ) return success; // TODO move upper
1228 
1229  // foreach_bit( int i, success ) {
1230  // cout << i << " " << mS0[i] << " " << mS1[i] << " " << mS2[i] << " " << (mS0 * z0 * z0)[i] << " " << (2 * z0 * z1 * mS1)[i] << " " << (mS2 * z1 * z1)[i] << " " << fChi2[i] << endl;
1231  // }
1232  // std::cout << success << std::endl;
1233  fP[ 0](success) += k00 * z0 + k01 * z1;
1234  fP[ 1](success) += k10 * z0 + k11 * z1;
1235  fP[ 2](success) = sinPhi ;
1236  fP[ 3](success) += k30 * z0 + k31 * z1;
1237  fP[ 4](success) += k40 * z0 + k41 * z1;
1238  fC[ 0](success) -= (k00 * c00 + k01 * c10); //c00
1239 
1240  fC[ 1](success) -= (k10 * c00 + k11 * c10); //c10
1241  fC[ 2](success) -= (k10 * c10 + k11 * c11); //c11
1242 
1243  fC[ 3](success) -= (k20 * c00 + k21 * c10); //c20
1244  fC[ 4](success) -= (k20 * c10 + k21 * c11); //c21
1245  fC[ 5](success) -= (k20 * c20 + k21 * c21); //c22
1246 
1247  fC[ 6](success) -= (k30 * c00 + k31 * c10); //c30
1248  fC[ 7](success) -= (k30 * c10 + k31 * c11); //c31
1249  fC[ 8](success) -= (k30 * c20 + k31 * c21); //c32
1250  fC[ 9](success) -= (k30 * c30 + k31 * c31); //c33
1251 
1252  fC[10](success) -= (k40 * c00 + k41 * c10); //c40
1253  fC[11](success) -= (k40 * c10 + k41 * c11); //c41
1254 // std::cout << fC[11] << " " << fC[10] << std::endl;
1255 
1256  fC[12](success) -= (k40 * c20 + k41 * c21); //c42
1257  fC[13](success) -= (k40 * c30 + k41 * c31); //c43
1258  fC[14](success) -= (k40 * c40 + k41 * c41); //c44
1259 
1260  return success;
1261 }
1262 
1263 float_m PndFTSCATrackParamVector::FilterWithMaterial( const float_v &y, const float_v &z, const FTSCAStripInfo &info, float_v err2, float maxSinPhi, const float_m &mask, const float_v& chi2Cut )
1264  //* Adds the 1-D measurement with the Kalman filter
1265  //* @beta is angle between the strip and z-axis, clockwise. The strip equation {y',z'} is: (y'-y)*cosB + (z'-z)*sinB = 0
1266 {
1267  assert( maxSinPhi > 0.f );
1268  float_m success = mask;
1269  success &= ( err2 > 1.e-8f );
1270 
1271  const float_v& c00 = fC[0];
1272  const float_v& c10 = fC[1];
1273  const float_v& c11 = fC[2];
1274 
1275  const float_v& sb = info.sin;
1276  const float_v& cb = info.cos;
1277 
1278  const float_v u = cb*y + sb*z;
1279  const float_v zeta = cb*fP[0] + sb*fP[1] - u;
1280 
1281  // F = CH'
1282  const float_v F0 = cb*c00 + sb*c10;
1283  const float_v F1 = cb*c10 + sb*c11;
1284 
1285  const float_v HCH = ( F0*cb + F1*sb );
1286 
1287  const float_v wi = 1/( err2 + HCH );
1288  const float_v zetawi = zeta *wi;
1289 
1290  fChi2(success) += zeta * zetawi;
1291  success &= fChi2 < chi2Cut;
1292  if ( success.isEmpty() ) return success;
1293 
1294  const float_v& c20 = fC[3];
1295  const float_v& c21 = fC[4];
1296  const float_v& c30 = fC[6];
1297  const float_v& c31 = fC[7];
1298  const float_v& c40 = fC[10];
1299  const float_v& c41 = fC[11];
1300 
1301  const float_v F2 = cb*c20 + sb*c21;
1302  const float_v F3 = cb*c30 + sb*c31;
1303  const float_v F4 = cb*c40 + sb*c41;
1304 
1305 
1306  const float_v K1 = F1*wi;
1307  const float_v K2 = F2*wi;
1308  const float_v K3 = F3*wi;
1309  const float_v K4 = F4*wi;
1310 
1311 
1312  const float_v sinPhi = fP[2] - F2*zetawi;
1313 
1314  success &= CAMath::Abs( sinPhi ) < maxSinPhi;
1315 
1316  fNDF( static_cast<int_m>(success) ) += 1; // TODO
1317 
1318  fP[ 0](success) -= F0*zetawi;
1319  fP[ 1](success) -= F1*zetawi;
1320  fP[ 2](success) = sinPhi ;
1321  fP[ 3](success) -= F3*zetawi;
1322  fP[ 4](success) -= F4*zetawi;
1323  fC[ 0](success) -= F0*F0*wi;
1324 
1325  fC[ 1](success) -= K1*F0;
1326  fC[ 2](success) -= K1*F1;
1327 
1328  fC[ 3](success) -= K2*F0;
1329  fC[ 4](success) -= K2*F1;
1330  fC[ 5](success) -= K2*F2;
1331 
1332  fC[ 6](success) -= K3*F0;
1333  fC[ 7](success) -= K3*F1;
1334  fC[ 8](success) -= K3*F2;
1335  fC[ 9](success) -= K3*F3;
1336 
1337  fC[10](success) -= K4*F0;
1338  fC[11](success) -= K4*F1;
1339 
1340  fC[12](success) -= K4*F2;
1341  fC[13](success) -= K4*F3;
1342  fC[14](success) -= K4*F4;
1343 
1344  return success;
1345 }
1346 
1347 float_m PndFTSCATrackParamVector::FilterWithMaterial( const float_v &y0, const float_v &z0, const float_v &r, const FTSCAStripInfo &info, float_v err2, float maxSinPhi, const float_m &mask, const float_v& chi2Cut )
1348  //* Adds the tube measurement with the Kalman filter
1349  //* @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 }
1350 {
1351  // linearize track in current point, which must be x == x0
1352  // 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]|
1353  const float_v& etx = sqrt( 1 - SinPhi()*SinPhi() );
1354  const float_v& ety = SinPhi();
1355  const float_v& etz = 1.f;
1356  const float_v& eox = - ety * info.cos - etz * info.sin; // e_o[x] = - e_t[y] * e_s[z] + e_t[z] * e_s[y]
1357  const float_v& eoy = etx * info.cos; // e_o[y] = - 0 + e_t[x] * e_s[z]
1358  const float_v& eoz = etx * info.sin; // e_o[z] = - e_t[x] * e_s[y] + 0
1359 
1360  const float_v& iEo = rsqrt( eox*eox + eoy*eoy + eoz*eoz );
1361  const float_v& h0 = eoy*iEo;
1362  const float_v& h1 = eoz*iEo;
1363 
1364  assert( maxSinPhi > 0.f );
1365  float_m success = mask;
1366  success &= ( err2 > 1.e-8f );
1367 
1368  const float_v& c00 = fC[0];
1369  const float_v& c10 = fC[1];
1370  const float_v& c11 = fC[2];
1371 
1372  const float_v zeta = h0*(fP[0] - y0) + h1*(fP[1] - z0) - r;
1373 
1374  // F = CH'
1375  const float_v F0 = h0*c00 + h1*c10;
1376  const float_v F1 = h0*c10 + h1*c11;
1377 
1378  const float_v HCH = ( F0*h0 + F1*h1 );
1379 
1380  const float_v wi = 1/( err2 + HCH );
1381  const float_v zetawi = zeta *wi;
1382 
1383  fChi2(success) += zeta * zetawi;
1384  success &= fChi2 < chi2Cut;
1385  if ( success.isEmpty() ) return success;
1386 
1387  const float_v& c20 = fC[3];
1388  const float_v& c21 = fC[4];
1389  const float_v& c30 = fC[6];
1390  const float_v& c31 = fC[7];
1391  const float_v& c40 = fC[10];
1392  const float_v& c41 = fC[11];
1393 
1394  const float_v F2 = h0*c20 + h1*c21;
1395  const float_v F3 = h0*c30 + h1*c31;
1396  const float_v F4 = h0*c40 + h1*c41;
1397 
1398 
1399  const float_v K1 = F1*wi;
1400  const float_v K2 = F2*wi;
1401  const float_v K3 = F3*wi;
1402  const float_v K4 = F4*wi;
1403 
1404 
1405  const float_v sinPhi = fP[2] - F2*zetawi;
1406 
1407  success &= CAMath::Abs( sinPhi ) < maxSinPhi;
1408 
1409  fNDF( static_cast<int_m>(success) ) += 1; // TODO
1410 
1411  fP[ 0](success) -= F0*zetawi;
1412  fP[ 1](success) -= F1*zetawi;
1413  fP[ 2](success) = sinPhi ;
1414  fP[ 3](success) -= F3*zetawi;
1415  fP[ 4](success) -= F4*zetawi;
1416  fC[ 0](success) -= F0*F0*wi;
1417 
1418  fC[ 1](success) -= K1*F0;
1419  fC[ 2](success) -= K1*F1;
1420 
1421  fC[ 3](success) -= K2*F0;
1422  fC[ 4](success) -= K2*F1;
1423  fC[ 5](success) -= K2*F2;
1424 
1425  fC[ 6](success) -= K3*F0;
1426  fC[ 7](success) -= K3*F1;
1427  fC[ 8](success) -= K3*F2;
1428  fC[ 9](success) -= K3*F3;
1429 
1430  fC[10](success) -= K4*F0;
1431  fC[11](success) -= K4*F1;
1432 
1433  fC[12](success) -= K4*F2;
1434  fC[13](success) -= K4*F3;
1435  fC[14](success) -= K4*F4;
1436 
1437  return success;
1438 }
1439 
1440 float_m PndFTSCATrackParamVector::Transport( const int_v& ista, const PndFTSCAParam& param, const float_m &mask )
1441 {
1442  float_m active = mask;
1443 
1444  PndFTSCATrackFitParam fitPar;
1445  CalculateFitParameters( fitPar );
1446  PndFTSCATrackLinearisationVector tE( *this );
1447  active &= TransportToX0WithMaterial( param.GetX0( ista, active ), tE, fitPar, param.GetXOverX0(int_v(ista),active), param.GetXTimesRho(int_v(ista),active), param.cBz(), 0.999f, active );
1448  return active;
1449 }
1450 
1451 float_m PndFTSCATrackParamVector::Transport( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_m &mask )
1452 {
1453  if ( ((CAMath::Abs(-fAlpha + hit.Angle()) < 1e-6f && CAMath::Abs(-fX + hit.X0()) < 2e-4f) || !mask).isFull() ) return mask;
1454  float_m active = mask & hit.IsValid();
1455  PndFTSCATrackLinearisationVector tR( *this );
1456  active &= Rotate( -fAlpha + hit.Angle(), tR, .999f, active );
1457 
1458  PndFTSCATrackFitParam fitPar;
1459  CalculateFitParameters( fitPar );
1460  active &= TransportToX0WithMaterial( hit.X0(), tR, fitPar, param.GetXOverX0(hit.IStations(),active), param.GetXTimesRho(hit.IStations(),active), param.cBz(), 0.999f, active );
1461  return active;
1462 }
1463 
1464 float_m PndFTSCATrackParamVector::Filter( const FTSCAHitV& hit, const PndFTSCAParam& param, const float_m &mask, const float_v& chi2Cut )
1465 {
1466  const int iSta = hit.IStations()[hit.IsValid().firstOne()];
1467  if ( param.Station( iSta ).NDF == 1 ) {
1468 #ifdef DRIFT_TUBES
1469  return FilterWithMaterial( hit.XWire1(), hit.XWire2(), hit.RSigned(), param.Station( iSta ).f, hit.Err2R(), 0.999f, mask, chi2Cut );
1470 #else
1471  assert(0);
1472 #endif
1473  }
1474  return FilterWithMaterial( hit.X1(), hit.X2(), hit.Err2X1(), hit.ErrX12(), hit.Err2X2(), 0.999f, mask, int_v(param.Station( hit.IStations()[0] ).NDF), chi2Cut );
1475 }
1476 
1477 float_m PndFTSCATrackParamVector::Transport( const FTSCAHit& hit, const PndFTSCAParam& param, const float_m &mask )
1478 {
1479  if ( ((CAMath::Abs(-fAlpha + hit.Angle()) < 1e-6f && CAMath::Abs(-fX + hit.X0()) < 2e-4f) || !mask).isFull() ) return mask;
1480 
1481  float_m active = mask;
1482 
1483  PndFTSCATrackLinearisationVector tR( *this );
1484  active &= Rotate( -fAlpha + hit.Angle(), tR, .999f, active );
1485 
1486  PndFTSCATrackFitParam fitPar;
1487  CalculateFitParameters( fitPar );
1488  active &= TransportToX0WithMaterial( hit.X0(), tR, fitPar, param.GetXOverX0(hit.IStation()), param.GetXTimesRho(hit.IStation()), param.cBz(), 0.999f, active );
1489  return active;
1490 }
1491 
1492 float_m PndFTSCATrackParamVector::Filter( const FTSCAHit& hit, const PndFTSCAParam& param, const float_m &mask, const float_v& chi2Cut )
1493 {
1494  if ( param.Station( hit.IStation() ).NDF == 1 ) {
1495 #ifdef DRIFT_TUBES
1496  return FilterWithMaterial( hit.XWire1(), hit.XWire2(), hit.RSigned(), param.Station( hit.IStation() ).f, hit.Err2R(), 0.999f, mask, chi2Cut );
1497 #else
1498  assert(0);
1499 #endif
1500  }
1501  return FilterWithMaterial( hit.X1(), hit.X2(), hit.Err2X1(), hit.ErrX12(), hit.Err2X2(), 0.999f, mask, int_v(param.Station( hit.IStation() ).NDF), chi2Cut );
1502 }
1503 
1504 float_m PndFTSCATrackParamVector::AddTarget( const FTSCATarget& target, const float_m &mask )
1505 {
1506  float_m active = mask;
1507 
1508  float_v eY, eZ;
1509  float_v dz = target.X0() - X0();
1510  float_v J[6];
1511  // H = 1 0 J[0] J[1] J[2]
1512  // 0 1 J[3] J[4] J[5]
1513  active &= TransportJ0ToX0( target.X0(), static_cast<float_v>(target.B()), eY, eZ, J, active );
1514  active &= FilterVtx( target.X1(), target.X2(), CAX1X2MeasurementInfo( target.Err2X1(), 0.f, target.Err2X2() ), eY, eZ, J, active );
1515  fNDF (static_cast<int_m>(active)) += target.NDF();
1516 
1517  return active;
1518 }
1519 
1520 #include <iostream>
1521 
1522 std::istream &operator>>( std::istream &in, PndFTSCATrackParamVector &t )
1523 {
1524  float_v::Memory x, s, p[5], c[15], chi2;
1525  int_v::Memory ndf;
1526  for ( int j = 0; j < uint_v::Size; ++j ) {
1527  in >> x[j];
1528  in >> s[j];
1529  for ( int i = 0; i < 5; i++ ) in >> p[i][j];
1530  for ( int i = 0; i < 15; i++ ) in >> c[i][j];
1531  in >> chi2[j];
1532  in >> ndf[j];
1533  }
1534  t.fX.load( x );
1535  t.fSignCosPhi.load( s );
1536  for ( int i = 0; i < 5; i++ ) t.fP[i].load( p[i] );
1537  for ( int i = 0; i < 5; i++ ) t.fC[i].load( c[i] );
1538  t.fChi2.load( chi2 );
1539  t.fNDF.load( ndf );
1540  return in;
1541 }
1542 
1543 std::ostream &operator<<( std::ostream &out, const PndFTSCATrackParamVector &t )
1544 {
1545  if ( out == std::cerr ) {
1546  out << "------------------------------ Track Param ------------------------------"
1547  << "\n X: " << t.X()
1548  << "\n SignCosPhi: " << t.SignCosPhi()
1549  << "\n Chi2: " << t.Chi2()
1550  << "\n NDF: " << t.NDF()
1551  << "\n Y: " << t.Par()[0]
1552  << "\n Z: " << t.Par()[1]
1553  << "\n SinPhi: " << t.Par()[2]
1554  << "\n DzDs: " << t.Par()[3]
1555  << "\n q/Pt: " << t.Par()[4]
1556  << "\nCovariance Matrix\n";
1557  int i = 0;
1558  out << std::setprecision( 2 );
1559  for ( int step = 1; step <= 5; ++step ) {
1560  int end = i + step;
1561  for ( ; i < end; ++i ) {
1562  out << t.Cov()[i] << '\t';
1563  }
1564  out << "\n";
1565  }
1566  out << std::setprecision( 6 );
1567  return out << std::endl;
1568  }
1569  for ( int j = 0; j < uint_v::Size; ++j ) {
1570  out << t.X()[j] << " "
1571  << t.SignCosPhi()[j] << " "
1572  << t.Chi2()[j] << " "
1573  << t.NDF()[j]
1574  << std::endl;
1575  for ( int i = 0; i < 5; i++ ) out << t.Par()[i][j] << " ";
1576  out << std::endl;
1577  for ( int i = 0; i < 15; i++ ) out << t.Cov()[i][j] << " ";
1578  out << std::endl;
1579  }
1580  return out;
1581 }
1582 
1583 #endif // !PANDA_FTS
CAMath::ASin
static T ASin(const T &x)
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
FTSCAHit::X2
float X2() const
Definition: FTSCAHits.h:43
FTSCAHitV::X2
float_v X2() const
Definition: FTSCAHitsV.h:51
p
Double_t p
Definition: anasim.C:58
std::basic_istream
Definition: PndFTSCATrackParamVector.h:1256
mK
double mK
Definition: softtrigger_kin5.C:128
dy
double dy
Definition: runGeaneTestSimBox.cxx:45
FTSCAHitV::Err2X1
float_v Err2X1() const
Definition: FTSCAHitsV.h:56
c11
TCanvas * c11
Definition: plot_acceptance_stt.C:188
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
PndFTSCAParam::Station
const FTSCAStation & Station(short i) const
Definition: PndFTSCAParam.h:45
PndFTSCAParam::GetStripInfo
void GetStripInfo(FTSCAStripInfoVector &stripInfo, const int_v iStation, const float_m &mask) const
Definition: PndFTSCAParam.cxx:148
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
FTSCAStation::f
FTSCAStripInfo f
Definition: FTSCAStation.h:32
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
FTSCAHit::IStation
char IStation() const
Definition: FTSCAHits.h:33
FTSCATarget::X1
float X1() const
Definition: FTSCATarget.h:39
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
FTSCATarget::NDF
int NDF() const
Definition: FTSCATarget.h:48
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
FTSCAHit::Angle
float Angle() const
Definition: FTSCAHits.h:71
PndFTSCATrackParamVector::SetSignCosPhi
void SetSignCosPhi(const float_v &v)
Definition: PndFTSCATrackParamVector.h:1411
PndFTSCADisplay::DrawGBPoint
void DrawGBPoint(float x, float y, float z, int Start=1, Size_t width=1)
Definition: PndFTSCADisplay.cxx:532
h4
TH1F * h4
Definition: anal_point_fwendcap.C:84
s
TLorentzVector s
Definition: Pnd2DStar.C:50
c10
TCanvas * c10
Definition: plot_invariantmass.C:42
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
par
Double_t par[3]
Definition: anaLmdCluster.C:146
c21
TCanvas * c21
Definition: plot_acceptance_stt.C:247
PndFTSCAParam::GetX0
float GetX0(short iSt) const
Definition: PndFTSCAParam.h:61
FTSCAHitV::IsValid
float_m IsValid() const
Definition: FTSCAHitsV.h:38
operator<<
std::ostream & operator<<(std::ostream &o, const PndEventInfo &a)
Definition: PndEventInfo.cxx:68
CAMath::Round
static T Round(const T &x)
PndFTSCADisplay::Instance
static PndFTSCADisplay & Instance()
Definition: PndFTSCADisplay.cxx:98
PndFTSCATrackParamVector::SetSinPhi
void SetSinPhi(const float_v &v)
Definition: PndFTSCATrackParamVector.h:1404
PndFTSCAParam::GetXOverX0
float GetXOverX0(short iSt) const
Definition: PndFTSCAParam.h:65
Zero
static const fvec Zero
Definition: KFParticleBaseSIMD.cxx:23
dz
Double_t dz
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:24
FTSCAHitV::X1
float_v X1() const
Definition: FTSCAHitsV.h:50
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
FTSCAHitV::Angle
float_v Angle() const
Definition: FTSCAHitsV.h:73
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
FTSCAHit::X0
float X0() const
Definition: FTSCAHits.h:41
PndFTSCATrackParamVector::AddTarget
float_m AddTarget(const FTSCATarget &target, const float_m &mask=float_m(true))
Definition: PndFTSCATrackParamVector.cxx:1504
FTSCAHit::X1
float X1() const
Definition: FTSCAHits.h:42
std::basic_ostream
Definition: PndFTSCATrackParamVector.h:1258
h1
TH1F * h1
Definition: anal_hit_digi_cluster_fwendcap.C:121
dE
Double_t dE
Definition: anasim.C:58
FTSCATarget::Err2X1
float Err2X1() const
Definition: FTSCATarget.h:42
FTSCAHitV::IStations
int_v IStations() const
Definition: FTSCAHitsV.h:45
FTSCATarget::B
const CAFieldValue & B(int i=0) const
Definition: FTSCATarget.h:46
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
FTSCAHit::Err2X2
float Err2X2() const
Definition: FTSCAHits.h:52
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
b2
Double_t b2
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:22
PndFTSCATrackParamVector::GetDistXZ2
float_v GetDistXZ2(const PndFTSCATrackParamVector &t) const
Definition: PndFTSCATrackParamVector.cxx:710
a
Int_t a
Definition: anaLmdDigi.C:126
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
PndFTSCAParam::cBz
float cBz() const
Definition: PndFTSCAParam.h:49
CAMath::Min
static T Min(const T &x, const T &y)
Definition: PndCAMath.h:35
FTSCATarget::X2
float X2() const
Definition: FTSCATarget.h:40
FTSCATarget::ZB
const float_v & ZB(int i) const
Definition: FTSCATarget.h:47
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
y0
Double_t y0
Definition: checkhelixhit.C:71
rsqrt
friend F32vec4 rsqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:32
CAMath::ATan2
static T ATan2(const T &y, const T &x)
f
TFile * f
Definition: bump_analys.C:12
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
out
TFile * out
Definition: reco_muo.C:20
FTSCATarget::X0
float X0() const
Definition: FTSCATarget.h:38
SetField
fRun SetField(fField)
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
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
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
PndFTSCAParam::GetXTimesRho
float GetXTimesRho(short iSt) const
Definition: PndFTSCAParam.h:66
c20
TCanvas * c20
Definition: plot_acceptance_stt.C:241
CAMath::Max
static T Max(const T &x, const T &y)
Definition: PndCAMath.h:36
fY
Double_t fY
Definition: PndCaloDraw.cxx:34
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
CAMath::Log
static T Log(const T &x)
Definition: PndCAMath.h:40
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
CAMath::TwoPi
static float TwoPi()
Definition: PndCAMath.h:61
L1FieldRegion::Set
void Set(const CAFieldValue &B0, const float_v B0z, const CAFieldValue &B1, const float_v B1z, const CAFieldValue &B2, const float_v B2z)
Definition: L1Field.h:147
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
FTSCATarget::Err2X2
float Err2X2() const
Definition: FTSCATarget.h:43
FTSCAHitV::X0
float_v X0() const
Definition: FTSCAHitsV.h:49
PndFTSCATrackParamVector::InitByTarget
void InitByTarget(const FTSCATarget &target)
Definition: PndFTSCATrackParamVector.cxx:619
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
y
Double_t y
Definition: createRootFscGeometryFile.C:159
alpha
double alpha
Definition: f_Init.h:9
FTSCAHitV::Err2X2
float_v Err2X2() const
Definition: FTSCAHitsV.h:58
FTSCAHitV::ErrX12
float_v ErrX12() const
Definition: FTSCAHitsV.h:57
FTSCAHit::Err2X1
float Err2X1() const
Definition: FTSCAHits.h:50
PndFTSCATrackParamVector::GetDist2
float_v GetDist2(const PndFTSCATrackParamVector &t) const
Definition: PndFTSCATrackParamVector.cxx:700
One
static const fvec One
Definition: KFParticleBaseSIMD.cxx:24
FTSCATarget::Err2QMom
float Err2QMom() const
Definition: FTSCATarget.h:44
PndFTSCATrackParamVector::ConvertTrackParamToVector
void ConvertTrackParamToVector(PndFTSCATrackParam t0[float_v::Size], int nTracksV)
Definition: PndFTSCATrackParamVector.cxx:557
debugKF
PndFTSCANoDebugStream & debugKF()
Definition: debug.h:74
FTSCAHit::ErrX12
float ErrX12() const
Definition: FTSCAHits.h:51
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
Generated on Wed Apr 3 2019 08:43:56 for FairRoot/PandaRoot by   doxygen 1.8.5