BOSS 7.0.6
BESIII Offline Software System
Loading...
Searching...
No Matches
Ext_xp_err.cxx
Go to the documentation of this file.
1// File: Ext_xp_err.cc
2//
3// Descrition: Handle extrapolation error matrix( x, y, z, px, py, pz ).
4// The used coordinate system is the cartesian BESIII coordinate system.
5//
6// This takes a base class, Ext_errmx.
7//
8// The initialization has to be done by Ext_xp_err::set_err.
9//
10//
11// Modified from BELLE by Wang Liangliang
12//
13// Data: 2005.3.30
14//
15
16#include <math.h>
17#include <stdlib.h>
18#include "TrkExtAlg/Ext_xp_err.h"
19
20
21
22// Constants.
23
24static const double Cinv( 3335.640952 );// 1/c in cm*GeV/KGauss unit.
25//static const double C( 0.0002997924582 );// Light speed in KGauss/(cm*GeV) unit.
26static const double C( 0.002997924582 );// Light speed in tesla/(cm*Gev) unit. (by Wangll)
27static const double Small( 0.01 );// small number.
28static const double Eps(1.0e-12);// Infinitesimal number.
29static const double Infinite( 1.0e10 );// Infinite number.
30
31static const int Ndim_herr( 5 );// Dimension of helix error matrix.
32
33// Default constructor
34Ext_xp_err::Ext_xp_err() :
35 m_xv(3), m_pv(3), m_xp_jcb(Ndim_err,Ndim_err,0),
36 m_h2xp_jcb(Ndim_err,Ndim_herr,0), m_q(0), m_mass2(0){}
37
38// Copy constructor
39Ext_xp_err::Ext_xp_err( const Ext_xp_err &err ) :
40 Ext_errmx( (Ext_errmx)err ), m_xv(err.m_xv), m_pv(err.m_pv),
41 m_xp_jcb(err.m_xp_jcb), m_h2xp_jcb(err.m_h2xp_jcb), m_q( err.m_q ),
42 m_mass2( err.m_mass2 ){}
43
44/*
45 Set error matrix( symmetry matrix ). This initializes the error matrix.
46*/
47void Ext_xp_err::set_err( const HepSymMatrix &err, const Hep3Vector &xv,
48 const Hep3Vector &pv,
49 const double &q, const double &mass )
50// (Inputs)
51// err -- Track error (symmetry) matrix for (x,y,z,px,py,pz) form.
52// xv -- coordinate(x,y,z) at the initialization.
53// pv -- momentum(px,py,pz) at the initialization.
54// q -- charge of the particle, either 1 or -1.
55// mass - Mass of the particle in GeV/c**2 unit.
56//
57{
58 Ext_errmx *errmx_ptr = (Ext_errmx *)this;
59 errmx_ptr->put_err( err );
60
61 m_xv = xv;
62 m_pv = pv;
63 m_q = q;
64 m_mass2 = mass * mass;
65
66// We do NOT check the validity of the error matrix here.
67
68}
69
70
71
72/*
73 Transform the error matrix for the transformation: x -> x'
74*/
75bool Ext_xp_err::move( const Hep3Vector &xv1, const Hep3Vector &pv1,
76 const Hep3Vector &B, const int ms_on,
77 const double chi_cc)
78{
79// (Inputs)
80// xv1 -- Coordinate(x',y',z') after transformation.
81// pv1 -- Momentum(px',py',pz') after transformation.
82// B -- B-field(Bx,By,Bz) at the transformation.
83// ms_on -- Flag to switch on/off the multiple scattering effect.
84// ms_on = 0 for switch off. ms_on = 1 for switch on.
85// chi_cc -- Constant of Moliere theory.
86//
87// (Outputs)
88// return - = 1 for success, = 0 for failed.
89//
90 double dx( ( xv1 - m_xv ).mag() );
91 double dp( ( pv1 - m_pv ).mag() );
92 double p2( m_pv.mag2() );
93 double p_abs( sqrt( p2 ) );
94
95 double p_inv;
96 if( p_abs > Small && pv1.mag() > Small ){
97 p_inv = 1.0 / p_abs;
98 } else {
99 p_inv = Infinite;
100 return 0;
101 }
102
103 double ms_coeff( 2.557 * chi_cc );
104 bool with_B( ( B.mag() > Small ) ? 1 : 0 );
105 double p2inv( p_inv * p_inv );
106 double p3inv( p2inv * p_inv );
107 double fdx( dx * p3inv );
108 double cx( 100.*m_q * C * fdx );//*100 due to units problem by L.L.Wang
109 double fdp( dp * p_inv );
110
111 double px( m_pv.x() );
112 double py( m_pv.y() );
113 double pz( m_pv.z() );
114 double px2( px * px );
115 double py2( py * py );
116 double pz2( pz * pz );
117 double pxy( px * py );
118 double pyz( py * pz );
119 double pzx( pz * px );
120 double Bx( B.x() );
121 double By( B.y() );
122 double Bz( B.z() );
123
124 m_xp_jcb( 1, 1 ) = 1.0; // dx'/dx
125 m_xp_jcb( 1, 4 ) = fdx * ( py2 + pz2 ); // dx'/dpx
126 m_xp_jcb( 1, 5 ) = - fdx * pxy; // dx'/dpy
127 m_xp_jcb( 1, 6 ) = - fdx * pzx; // dx'/dpz
128
129 m_xp_jcb( 2, 2 ) = 1.0; // dy'/dy
130 m_xp_jcb( 2, 4 ) = - fdx * pxy; // dy'/dpx
131 m_xp_jcb( 2, 5 ) = fdx * ( pz2 + px2 ); // dy'/dpy
132 m_xp_jcb( 2, 6 ) = - fdx * pyz; // dy'/dpz
133
134 m_xp_jcb( 3, 3 ) = 1.0; // dz'/dz
135 m_xp_jcb( 3, 4 ) = - fdx * pzx; // dz'/dpx
136 m_xp_jcb( 3, 5 ) = - fdx * pyz; // dz'/dpy
137 m_xp_jcb( 3, 6 ) = fdx * ( px2 + py2 ); // dz'/dpz
138
139 m_xp_jcb( 4, 4 ) = 1.0 - fdp; // dx'/dx energy loss
140 m_xp_jcb( 5, 5 ) = 1.0 - fdp; // dy'/dy energy loss
141 m_xp_jcb( 6, 6 ) = 1.0 - fdp; // dz'/dz energy loss
142
143 if( with_B && m_q!=0. ){ // B != 0 and q !=0 case
144 m_xp_jcb( 4, 4 ) += - cx * ( pxy * Bz - pzx * By ); // dpx'/dpx
145 m_xp_jcb( 4, 5 ) = cx * ( ( pz2 + px2 ) * Bz + pyz * By ); // dpx'/dpy
146 m_xp_jcb( 4, 6 ) = - cx * ( ( px2 + py2 ) * By + pyz * Bz ); // dpx'/dpz
147
148 m_xp_jcb( 5, 4 ) = - cx * ( ( py2 + pz2 ) * Bz + pzx * Bx ); // dpy'/dpx
149 m_xp_jcb( 5, 5 ) += - cx * ( pyz * Bx - pxy * Bz ); // dpy'/dpy
150 m_xp_jcb( 5, 6 ) = cx * ( ( px2 + py2 ) * Bx + pzx * Bz ); // dpy'/dpz
151
152 m_xp_jcb( 6, 4 ) = cx * ( ( py2 + pz2 ) * By + pxy * Bx ); // dpz'/dpx
153 m_xp_jcb( 6, 5 ) = - cx * ( ( pz2 + px2 ) * Bx + pxy * By ); // dpz'/dpy
154 m_xp_jcb( 6, 6 ) += - cx * ( pzx * By - pyz * Bx ); // dpz'/dpz
155 }
156
157// error transformation.
158
159 m_err = m_err.similarity( m_xp_jcb );
160
161// Multiple scattering.
162
163 if( ms_on ){
164 double beta_p_inv( ( m_mass2 + p2 ) / ( p2 * p2 ) ); // 1/(p*beta)**2
165 double th2( 100000.0 * ms_coeff * ms_coeff * beta_p_inv * dx );//mutiply 100000.0 by L.L.Wang
166 //due to units problem
167 double m11( th2 * dx * dx /3.0);
168 double m12( 0.5 * th2 * dx * p_abs );
169 double m22( th2 * p2 );
170
171 if( p_abs > Eps ){
172 double c1( px*p_inv );
173 double c2( py*p_inv );
174 double c3( pz*p_inv );
175 double c12( - c1*c2 );
176 double c13( - c1*c3 );
177 double c23( - c2*c3 );
178 double s1s( 1 - c1*c1 );
179 double s2s( 1 - c2*c2 );
180 double s3s( 1 - c3*c3 );
181
182 m_err.fast( 1, 1 ) += m11*s1s; // error(x,x)
183
184 m_err.fast( 2, 1 ) += m11*c12; // error(y,x)
185 m_err.fast( 2, 2 ) += m11*s2s; // error(y,y)
186
187 m_err.fast( 3, 1 ) += m11*c13; // error(z,x)
188 m_err.fast( 3, 2 ) += m11*c23; // error(z,y)
189 m_err.fast( 3, 3 ) += m11*s3s; // error(z,z)
190
191 m_err.fast( 4, 1 ) += m12*s1s; // error(px,x)
192 m_err.fast( 4, 2 ) += m12*c12; // error(px,y)
193 m_err.fast( 4, 3 ) += m12*c13; // error(px,z)
194 m_err.fast( 4, 4 ) += m22*s1s; // error(px,px)
195
196 m_err.fast( 5, 1 ) += m12*c12; // error(py,x)
197 m_err.fast( 5, 2 ) += m12*s2s; // error(py,y)
198 m_err.fast( 5, 3 ) += m12*c23; // error(py,z)
199 m_err.fast( 5, 4 ) += m22*c12; // error(py,px)
200 m_err.fast( 5, 5 ) += m22*s2s; // error(py,py)
201
202 m_err.fast( 6, 1 ) += m12*c13; // error(pz,x)
203 m_err.fast( 6, 2 ) += m12*c23; // error(pz,y)
204 m_err.fast( 6, 3 ) += m12*s3s; // error(pz,z)
205 m_err.fast( 6, 4 ) += m22*c13; // error(pz,px)
206 m_err.fast( 6, 5 ) += m22*c23; // error(pz,py)
207 m_err.fast( 6, 6 ) += m22*s3s; // error(pz,pz)
208 }
209
210 }
211
212// Now store the new xv and pv.
213 m_xv = xv1;
214 m_pv = pv1;
215 return 1;
216}
217
218//Assign "=" operator
219Ext_xp_err &Ext_xp_err::operator=( const Ext_xp_err &err )
220{
221 if( this != &err ){
222 *((Ext_errmx *)this) = err;
223 m_xv = err.m_xv;
224 m_pv = err.m_pv;
225 m_xp_jcb = err.m_xp_jcb;
226 m_h2xp_jcb = err.m_h2xp_jcb;
227 m_q = err.m_q;
228 m_mass2 = err.m_mass2;
229 }
230 return *this;
231}
232
233/*
234 ostream "<<" operator
235*/
236std::ostream &operator<<( std::ostream &s, const Ext_xp_err &err )
237{
238 s << "m_err: " << (Ext_errmx)err
239 << std::endl
240 << "m_xv: " << err.m_xv << std::endl << " abs(xv)= " << err.m_xv.mag()
241 << std::endl
242 << "m_pv: " << err.m_pv << std::endl << " abs(pv)= " << err.m_pv.mag()
243 << std::endl
244 << "m_xp_jcb: " << err.m_xp_jcb
245 << std::endl
246 << "m_h2xp_jcb: " << err.m_h2xp_jcb
247 << std::endl
248 << "m_q: " << err.m_q
249 << std::endl
250 << "m_mass: " << sqrt(err.m_mass2)
251 << std::endl;
252 return s;
253}
mass
double mass
Definition: CosmicGenerator.cxx:138
operator<<
std::ostream & operator<<(std::ostream &s, const Ext_xp_err &err)
Definition: Ext_xp_err.cxx:236
Ext_xp_err.h
s
XmlRpcServer s
Definition: HelloServer.cpp:11
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition: KKsem.h:33
m_q
NTuple::Item< double > m_q
Definition: MdcHistItem.h:84
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition: RRes.h:29
Ext_errmx
Definition: Ext_errmx.h:36
Ext_errmx::m_err
HepSymMatrix m_err
Definition: Ext_errmx.h:133
Ext_errmx::put_err
void put_err(const double error[])
Definition: Ext_errmx.cxx:45
Ext_xp_err
Definition: Ext_xp_err.h:29
Ext_xp_err::operator=
Ext_xp_err & operator=(const Ext_xp_err &xp_err)
Definition: Ext_xp_err.cxx:219
Ext_xp_err::Ext_xp_err
Ext_xp_err()
Definition: Ext_xp_err.cxx:34
Ext_xp_err::set_err
void set_err(const HepSymMatrix &err, const Hep3Vector &xv, const Hep3Vector &pv, const double &q, const double &mass)
Definition: Ext_xp_err.cxx:47
Ext_xp_err::move
bool move(const Hep3Vector &xv1, const Hep3Vector &pv1, const Hep3Vector &B, const int ms_on, const double chi_cc)
Definition: Ext_xp_err.cxx:75