FairRoot/PandaRoot
PndFTSCATrackParam.h
Go to the documentation of this file.
1 
2  //-*- Mode: C++ -*-
3  // $Id: PndFTSCATrackParam.h,v 1.5 2011/05/20 16:28:05 fisyak Exp $
4  // ************************************************************************
5  // This file is property of and copyright by the ALICE HLT Project *
6  // ALICE Experiment at CERN, All rights reserved. *
7  // See cxx source for full Copyright notice *
8  // *
9  //*************************************************************************
10 
11 
12 #ifndef PNDFTSCATRACKPARAM_H
13 #define PNDFTSCATRACKPARAM_H
14 #include <string.h>
15 #include "PndFTSCAMath.h"
17 #include "PndFTSCADef.h"
18 #include "FTSCAHits.h"
19 
20 #ifdef PANDA_FTS
21 
30 {
33  public:
34 
36  PndFTSCATrackParam( const TrackParamVector &v, int i )
37  : x( v.X()[i] ), y( v.Y()[i] ), tx( v.Tx()[i] ), ty( v.Ty()[i] ), qp( v.QP()[i] ), z( v.Z()[i] ),
38  C00( v.Cov(0)[i] ),
39  C10( v.Cov(1)[i] ),
40  C11( v.Cov(2)[i] ),
41  C20( v.Cov(3)[i] ),
42  C21( v.Cov(4)[i] ),
43  C22( v.Cov(5)[i] ),
44  C30( v.Cov(6)[i] ),
45  C31( v.Cov(7)[i] ),
46  C32( v.Cov(8)[i] ),
47  C33( v.Cov(9)[i] ),
48  C40( v.Cov(10)[i] ),
49  C41( v.Cov(11)[i] ),
50  C42( v.Cov(12)[i] ),
51  C43( v.Cov(13)[i] ),
52  C44( v.Cov(14)[i] ),
53  chi2( v.Chi2()[i] ),
54  ndf( v.NDF()[i] ),
55  fAlpha( v.Angle()[i] )
56  {}
57 
58  float X0() const { return z; }
59  float X1() const { return x; }
60  float X2() const { return y; }
61  float Tx1() const { return tx; }
62  float Tx2() const { return ty; }
63 
64  float X() const { return x; }
65  float Y() const { return y; }
66  float Z() const { return z; }
67  float Tx() const { return tx; }
68  float Ty() const { return ty; }
69  float QP() const { return qp; }
70 
71  float Chi2() const { return chi2; }
72  int NDF() const { return ndf; }
73 
74  float QMomentum() const { return QP(); } // used for triplets comparison
75  float Angle() const { return fAlpha; }
76 
77  float Err2X1() const { return Err2X(); }
78  float Err2X2() const { return Err2Y(); }
79  float Err2Tx1() const { return Err2Tx(); }
80  float Err2Tx2() const { return Err2Ty(); }
81 
82  float Err2X() const { return C00; }
83  float Err2Y() const { return C11; }
84  float Err2Tx() const { return C22; }
85  float Err2Ty() const { return C33; }
86  float Err2QP() const { return C44; }
87  float Err2QMomentum() const { return Err2QP(); }
88 
89  float* Par() { return &x; }
90  float* Cov() { return &C00; }
91 
92  float Par( int i ) const {
93  switch ( i ) {
94  case 0: return x;
95  case 1: return y;
96  case 2: return tx;
97  case 3: return ty;
98  case 4: return qp;
99  }
100  assert(0); return -1;
101  }
102 
103  float Cov( int i ) const {
104  switch ( i ) {
105  case 0: return C00;
106  case 1: return C10;
107  case 2: return C11;
108  case 3: return C20;
109  case 4: return C21;
110  case 5: return C22;
111  case 6: return C30;
112  case 7: return C31;
113  case 8: return C32;
114  case 9: return C33;
115  case 10: return C40;
116  case 11: return C41;
117  case 12: return C42;
118  case 13: return C43;
119  case 14: return C44;
120  }
121  assert(0); return -1;
122  }
123 
124  void SetCov( int i, float v ) {
125  switch ( i ) {
126  case 0: C00 = v; break;
127  case 1: C10 = v; break;
128  case 2: C11 = v; break;
129  case 3: C20 = v; break;
130  case 4: C21 = v; break;
131  case 5: C22 = v; break;
132  case 6: C30 = v; break;
133  case 7: C31 = v; break;
134  case 8: C32 = v; break;
135  case 9: C33 = v; break;
136  case 10: C40 = v; break;
137  case 11: C41 = v; break;
138  case 12: C42 = v; break;
139  case 13: C43 = v; break;
140  case 14: C44 = v; break;
141  }
142  }
143 
144  void SetX( float v ) { x = v; }
145  void SetY( float v ) { y = v; }
146  void SetZ( float v ) { z = v; }
147 
148  void SetTX( float v ) { tx = v; }
149  void SetTY( float v ) { ty = v; }
150  void SetQP( float v ) { qp = v; }
151 
152  void SetChi2( float v ) { chi2 = v; }
153  void SetNDF( int v ) { ndf = v; }
154 
155 
156  void Print() const;
157 
158  bool TransportToX0( const float &x_, const float &Bz ) {
159  UNUSED_PARAM1(Bz);
160  return TransportToX0Line( x_ );
161  };
162 
163  bool Transport( const FTSCAHit& hit, const PndFTSCAParam& param )
164  {
166  v.ConvertTrackParamToVector( this, 1 );
167  const float_m &mask = v.Transport( hit, param );
168  *this = PndFTSCATrackParam( v, 0 );
169  return mask[0];
170  };
171 
172  bool Filter( const FTSCAHit& hit, const PndFTSCAParam& param )
173  {
175  v.ConvertTrackParamToVector( this, 1 );
176  const float_m &mask = v.Filter( hit, param );
177  *this = PndFTSCATrackParam( v, 0 );
178  return mask[0];
179  };
180 
181  void Reset() { x=y=tx=ty=qp=z=
182  C00=
183  C10= C11=
184  C20= C21= C22=
185  C30= C31= C32= C33=
186  C40= C41= C42= C43= C44=0.f;
187  ndf = 0;
188  chi2 = -1;}
189 
190  bool TransportToX0Line( float x0 );
191 
192  bool IsValid() const { return chi2 != -1; }
193  void SetAsInvalid(){ chi2 = -1; }
194 
195  bool Rotate( float ){ return 1; };// don't need rotation in CBM
196 
197 
198  private:
199 
200  float x,y,tx,ty,qp,z,
201  C00,
202  C10, C11,
203  C20, C21, C22,
204  C30, C31, C32, C33,
205  C40, C41, C42, C43, C44;
206 
207  float chi2; // the chi^2 value
208  short ndf; // the Number of Degrees of Freedom
209 
210  float fAlpha; // coor system
211 };
212 
213 inline bool PndFTSCATrackParam::TransportToX0Line( float x0_out )
214 {
215  float dz = (x0_out - z);
216 
217  x += tx*dz;
218  y += ty*dz;
219  z += dz;
220 
221  const float dzC32_in = dz * C32;
222 
223  C21 += dzC32_in;
224  C10 += dz * ( C21 + C30 );
225 
226  const float C20_in = C20;
227 
228  C20 += dz * C22;
229  C00 += dz * ( C20 + C20_in );
230 
231  const float C31_in = C31;
232 
233  C31 += dz * C33;
234  C11 += dz * ( C31 + C31_in );
235  C30 += dzC32_in;
236 
237  C40 += dz * C42;
238  C41 += dz * C43;
239 
240  return 1;
241 }
242 
243 #else // PANDA_FTS
244 
246 {
249  public:
250 
253  : fX( v.X()[i] ),
254  fSignCosPhi( v.SignCosPhi()[i] ),
255  fChi2( v.Chi2()[i] ),
256  fNDF( v.NDF()[i] ),
257  fAlpha( v.Angle()[i] )
258  {
259  for ( int j = 0; j < 5; ++j ) fP[j] = v.Par()[j][i];
260  for ( int j = 0; j < 15; ++j ) fC[j] = v.Cov()[j][i];
261  }
262 
263  float X0() const { return fX; }
264  float X1() const { return Y(); }
265  float X2() const { return Z(); }
266 
267  float X() const { return fX; }
268  float Y() const { return fP[0]; }
269  float Z() const { return fP[1]; }
270  float SinPhi() const { return fP[2]; }
271  float DzDs() const { return fP[3]; }
272  float QPt() const { return fP[4]; }
273  float QP() const { return QPt() / sqrt(DzDs()*DzDs() + 1); }
274  float QMomentum() const { return QPt(); }
275 
276  float SignCosPhi() const { return fSignCosPhi; }
277  float Chi2() const { return fChi2; }
278  int NDF() const { return fNDF; }
279 
280  float Err2X1() const { return Err2Y(); }
281  float Err2X2() const { return Err2Z(); }
282 
283  float Err2Y() const { return fC[0]; }
284  float Err2Z() const { return fC[2]; }
285  float Err2SinPhi() const { return fC[5]; }
286  float Err2DzDs() const { return fC[9]; }
287  float Err2QPt() const { return fC[14]; }
288 
289  float Angle() const { return fAlpha; }
290 
291  float Kappa( float Bz ) const { return fP[4]*Bz; }
292  float CosPhi() const { return fSignCosPhi*CAMath::Sqrt( 1 - SinPhi()*SinPhi() ); }
293 
294  float Err2QMomentum() const { return fC[14]; }
295 
296  const float *Par() const { return fP; }
297  const float *Cov() const { return fC; }
298 
299  bool GetXYZPxPyPzQ( float& x, float& y, float& z, float& px, float& py, float& pz, int& q, float cov[21] ) const;
300 
301  void SetSinPhi( float v ) { fP[2] = v; }
302 
303  void GetDCAPoint( float x, float y, float z,
304  float &px, float &py, float &pz, float Bz ) const;
305 
306 
307  bool TransportToX0( float x, float Bz, float maxSinPhi = .999 );
308 
309  bool Rotate( float alpha, float maxSinPhi = .999 );
310 
311  bool Transport( const FTSCAHit& hit, const PndFTSCAParam& param );
312 
313  bool IsValid() const { return fChi2 != -1; }
314  void SetAsInvalid(){ fChi2 = -1; }
315 
316  private:
317  void Reset() {
318  fX = 0;
319  fSignCosPhi = 0;
320  for(int i=0; i < 5; i++) fP[i] = 0;
321  for(int i=0; i < 15; i++) fC[i] = 0;
322  fChi2 = 0;
323  fNDF = 0;
324  }
325  float S( float x, float y, float Bz ) const;
326 
327  float fX; // x position
328  float fSignCosPhi; // sign of cosPhi // phi = arctg (Dy/Dx)
329  float fP[5]; // 'active' track parameters: Y, Z, SinPhi = dy/sqrt(dx*dx + dy*dy), DzDs = dz/sqrt(dx*dx + dy*dy), q/Pt
330  float fC[15]; // the covariance matrix for Y,Z,SinPhi,..
331  float fChi2; // the chi^2 value
332  int fNDF; // the Number of Degrees of Freedom
333 
334  float fAlpha; // coor system
335 };
336 
337 inline bool PndFTSCATrackParam::TransportToX0( float x, float Bz, float maxSinPhi )
338 {
339  //* Transport the track parameters to X=x, using linearization at t0, and the field value Bz
340  //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
341  //* linearisation of trajectory t0 is also transported to X=x,
342  //* returns 1 if OK
343  //*
344 
345  const float ex = CosPhi();
346  const float ey = SinPhi();
347  const float k = QPt() * Bz;
348  const float dx = x - X();
349 
350  const float ey1 = k * dx + ey;
351 
352  // check for intersection with X=x
353 
354  if ( CAMath::Abs( ey1 ) > maxSinPhi ) return 0;
355 
356  float ex1 = CAMath::Sqrt( 1.f - ey1 * ey1 );
357  if ( ex < 0 ) ex1 = -ex1;
358 
359  const float dx2 = dx * dx;
360  const float ss = ey + ey1;
361  const float cc = ex + ex1;
362 
363  if ( ( CAMath::Abs( cc ) < 1.e-4 || CAMath::Abs( ex ) < 1.e-4 || CAMath::Abs( ex1 ) < 1.e-4 ) ) return 0;
364 
365  const float cci = 1.f / cc;
366  const float exi = 1.f / ex;
367  const float ex1i = 1.f / ex1;
368 
369  const float tg = ss * cci; // tan((phi1+phi)/2)
370 
371  const float dy = dx * tg;
372  float dl = dx * CAMath::Sqrt( 1.f + tg * tg );
373 
374  if ( cc < 0 ) dl = -dl;
375  float dSin = dl * k * 0.5;
376  if ( dSin > 1.f ) dSin = 1.f;
377  if ( dSin < -1.f ) dSin = -1.f;
378  const float dS = ( CAMath::Abs( k ) > 1.e-4 ) ? ( 2 * CAMath::ASin( dSin ) / k ) : dl;
379  const float dz = dS * DzDs();
380 
381 
382  //float H0[5] = { 1,0, h2, 0, h4 };
383  //float H1[5] = { 0, 1, 0, dS, 0 };
384  //float H2[5] = { 0, 0, 1, 0, dxBz };
385  //float H3[5] = { 0, 0, 0, 1, 0 };
386  //float H4[5] = { 0, 0, 0, 0, 1 };
387 
388  const float h2 = dx * ( 1.f + ey * ey1 + ex * ex1 ) * exi * ex1i * cci;
389  const float h4 = dx2 * ( cc + ss * ey1 * ex1i ) * cci * cci * Bz;
390  const float dxBz = dx * Bz;
391 
392  fX = X() + dx;
393  fP[0] = Y() + dy;
394  fP[1] = Z() + dz;
395  fP[2] = ey1;
396 
397 #if 1
398  const float c00 = fC[0];
399  const float c10 = fC[1];
400  const float c11 = fC[2];
401  const float c20 = fC[3];
402  const float c21 = fC[4];
403  const float c22 = fC[5];
404  const float c30 = fC[6];
405  const float c31 = fC[7];
406  const float c32 = fC[8];
407  const float c33 = fC[9];
408  const float c40 = fC[10];
409  const float c41 = fC[11];
410  const float c42 = fC[12];
411  const float c43 = fC[13];
412  const float c44 = fC[14];
413 
414  fC[0] = ( c00 + h2 * h2 * c22 + h4 * h4 * c44
415  + 2.f * ( h2 * c20 + h4 * c40 + h2 * h4 * c42 ) );
416 
417  fC[1] = c10 + h2 * c21 + h4 * c41 + dS * ( c30 + h2 * c32 + h4 * c43 );
418  fC[2] = c11 + 2.f * dS * c31 + dS * dS * c33;
419 
420  fC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * ( c40 + h2 * c42 + h4 * c44 );
421  fC[4] = c21 + dS * c32 + dxBz * ( c41 + dS * c43 );
422  fC[5] = c22 + 2.f * dxBz * c42 + dxBz * dxBz * c44;
423 
424  fC[6] = c30 + h2 * c32 + h4 * c43;
425  fC[7] = c31 + dS * c33;
426  fC[8] = c32 + dxBz * c43;
427  fC[9] = c33;
428 
429  fC[10] = c40 + h2 * c42 + h4 * c44;
430  fC[11] = c41 + dS * c43;
431  fC[12] = c42 + dxBz * c44;
432  fC[13] = c43;
433  fC[14] = c44;
434 #else
435  fC[0] = ( fC[0] + h2 * h2 * fC[5] + h4 * h4 * fC[14]
436  + 2 * ( h2 * fC[3] + h4 * fC[10] + h2 * h4 * fC[12] ) );
437 
438  fC[1] = fC[1] + h2 * fC[4] + h4 * fC[11] + dS * ( fC[6] + h2 * fC[8] + h4 * fC[13] );
439  fC[2] = fC[2] + 2 * dS * fC[7] + dS * dS * fC[9];
440 
441  fC[3] = fC[3] + h2 * fC[5] + h4 * fC[12] + dxBz * ( fC[10] + h2 * fC[12] + h4 * fC[14] );
442  fC[4] = fC[4] + dS * fC[8] + dxBz * ( fC[11] + dS * fC[13] );
443  fC[5] = fC[5] + 2 * dxBz * fC[12] + dxBz * dxBz * fC[14];
444 
445  fC[6] = fC[6] + h2 * fC[8] + h4 * fC[13];
446  fC[7] = fC[7] + dS * fC[9];
447  fC[8] = fC[8] + dxBz * fC[13];
448  fC[9] = fC[9];
449 
450  fC[10] = fC[10] + h2 * fC[12] + h4 * fC[14];
451  fC[11] = fC[11] + dS * fC[13];
452  fC[12] = fC[12] + dxBz * fC[14];
453  fC[13] = fC[13];
454  fC[14] = fC[14];
455 #endif
456 
457  return 1;
458 }
459 
460 #endif // PANDA_FTS
461 
463 
466 
467 #endif
static T ASin(const T &x)
MechFsc Print()
Double_t x0
Definition: checkhelixhit.C:70
double dy
TCanvas * c11
const float_v & Cov(int i) const
bool Transport(const FTSCAHit &hit, const PndFTSCAParam &param)
Int_t i
Definition: run_full.C:25
void GetDCAPoint(float x, float y, float z, float &px, float &py, float &pz, float Bz) const
bool Rotate(float alpha, float maxSinPhi=.999)
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
float SignCosPhi() const
float Err2QPt() const
TH1F * h4
TCanvas * c10
PndCATrackParam TrackParam
const float * Par() const
TCanvas * c21
std::ostream & operator<<(std::ostream &o, const PndEventInfo &a)
friend std::istream & operator>>(std::istream &, PndFTSCATrackParam &)
friend std::ostream & operator<<(std::ostream &, const PndFTSCATrackParam &)
PndFTSCATrackParam(const TrackParamVector &v, int i)
__m128 v
Definition: P4_F32vec4.h:4
static T Abs(const T &x)
Definition: PndCAMath.h:39
const float * Cov() const
float_m Transport(const int_v &ista, const PndFTSCAParam &param, const float_m &mask=float_m(true))
bool GetXYZPxPyPzQ(float &x, float &y, float &z, float &px, float &py, float &pz, int &q, float cov[21]) const
float QMomentum() const
TFile * f
Definition: bump_analys.C:12
Double_t z
bool TransportToX0(float x, float Bz, float maxSinPhi=.999)
TFile * out
Definition: reco_muo.C:20
float Err2QMomentum() const
double dx
TCanvas * c22
TCanvas * c20
float Err2DzDs() const
TBuffer & operator>>(TBuffer &buf, PndAnaPidSelector *&obj)
Double_t x
const float_v & Par(int i) const
float S(float x, float y, float Bz) const
PndSdsMCPoint * hit
Definition: anasim.C:70
Double_t y
double alpha
Definition: f_Init.h:9
void ConvertTrackParamToVector(PndFTSCATrackParam t0[float_v::Size], int nTracksV)
float Err2SinPhi() const
double pz[39]
Definition: pipisigmas.h:14
float Kappa(float Bz) const
float_m Filter(const FTSCAHitV &hit, const PndFTSCAParam &param, const float_m &mask=float_m(true), const float_v &chi2Cut=10e10f)