BOSS 8.0.0
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 "Ext_xp_err.h"
17#include <math.h>
18#include <stdlib.h>
19
20// Constants.
21
22static const double Cinv( 3335.640952 ); // 1/c in cm*GeV/KGauss unit.
23// static const double C( 0.0002997924582 );// Light speed in KGauss/(cm*GeV) unit.
24static const double C( 0.002997924582 ); // Light speed in tesla/(cm*Gev) unit. (by Wangll)
25static const double Small( 0.01 ); // small number.
26static const double Eps( 1.0e-12 ); // Infinitesimal number.
27static const double Infinite( 1.0e10 ); // Infinite number.
28
29static const int Ndim_herr( 5 ); // Dimension of helix error matrix.
30
31// Default constructor
32Ext_xp_err::Ext_xp_err()
33 : m_xv( 3 )
34 , m_pv( 3 )
35 , m_xp_jcb( Ndim_err, Ndim_err, 0 )
36 , m_h2xp_jcb( Ndim_err, Ndim_herr, 0 )
37 , m_q( 0 )
38 , m_mass2( 0 ) {}
39
40// Copy constructor
41Ext_xp_err::Ext_xp_err( const Ext_xp_err& err )
42 : Ext_errmx( (Ext_errmx)err )
43 , m_xv( err.m_xv )
44 , m_pv( err.m_pv )
45 , m_xp_jcb( err.m_xp_jcb )
46 , m_h2xp_jcb( err.m_h2xp_jcb )
47 , m_q( err.m_q )
48 , m_mass2( err.m_mass2 ) {}
49
50/*
51 Set error matrix( symmetry matrix ). This initializes the error matrix.
52*/
53void Ext_xp_err::set_err( const HepSymMatrix& err, const Hep3Vector& xv, const Hep3Vector& pv,
54 const double& q, const double& mass )
55// (Inputs)
56// err -- Track error (symmetry) matrix for (x,y,z,px,py,pz) form.
57// xv -- coordinate(x,y,z) at the initialization.
58// pv -- momentum(px,py,pz) at the initialization.
59// q -- charge of the particle, either 1 or -1.
60// mass - Mass of the particle in GeV/c**2 unit.
61//
62{
63 Ext_errmx* errmx_ptr = (Ext_errmx*)this;
64 errmx_ptr->put_err( err );
65
66 m_xv = xv;
67 m_pv = pv;
68 m_q = q;
69 m_mass2 = mass * mass;
70
71 // We do NOT check the validity of the error matrix here.
72}
73
74/*
75 Transform the error matrix for the transformation: x -> x'
76*/
77bool Ext_xp_err::move( const Hep3Vector& xv1, const Hep3Vector& pv1, const Hep3Vector& B,
78 const int ms_on, const double chi_cc ) {
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 ) { p_inv = 1.0 / p_abs; }
97 else
98 {
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. )
144 { // B != 0 and q !=0 case
145 m_xp_jcb( 4, 4 ) += -cx * ( pxy * Bz - pzx * By ); // dpx'/dpx
146 m_xp_jcb( 4, 5 ) = cx * ( ( pz2 + px2 ) * Bz + pyz * By ); // dpx'/dpy
147 m_xp_jcb( 4, 6 ) = -cx * ( ( px2 + py2 ) * By + pyz * Bz ); // dpx'/dpz
148
149 m_xp_jcb( 5, 4 ) = -cx * ( ( py2 + pz2 ) * Bz + pzx * Bx ); // dpy'/dpx
150 m_xp_jcb( 5, 5 ) += -cx * ( pyz * Bx - pxy * Bz ); // dpy'/dpy
151 m_xp_jcb( 5, 6 ) = cx * ( ( px2 + py2 ) * Bx + pzx * Bz ); // dpy'/dpz
152
153 m_xp_jcb( 6, 4 ) = cx * ( ( py2 + pz2 ) * By + pxy * Bx ); // dpz'/dpx
154 m_xp_jcb( 6, 5 ) = -cx * ( ( pz2 + px2 ) * Bx + pxy * By ); // dpz'/dpy
155 m_xp_jcb( 6, 6 ) += -cx * ( pzx * By - pyz * Bx ); // dpz'/dpz
156 }
157
158 // error transformation.
159
160 m_err = m_err.similarity( m_xp_jcb );
161
162 // Multiple scattering.
163
164 if ( ms_on )
165 {
166 double beta_p_inv( ( m_mass2 + p2 ) / ( p2 * p2 ) ); // 1/(p*beta)**2
167 double th2( 100000.0 * ms_coeff * ms_coeff * beta_p_inv * dx ); // mutiply 100000.0 by
168 // L.L.Wang due to units
169 // problem
170 double m11( th2 * dx * dx / 3.0 );
171 double m12( 0.5 * th2 * dx * p_abs );
172 double m22( th2 * p2 );
173
174 if ( p_abs > Eps )
175 {
176 double c1( px * p_inv );
177 double c2( py * p_inv );
178 double c3( pz * p_inv );
179 double c12( -c1 * c2 );
180 double c13( -c1 * c3 );
181 double c23( -c2 * c3 );
182 double s1s( 1 - c1 * c1 );
183 double s2s( 1 - c2 * c2 );
184 double s3s( 1 - c3 * c3 );
185
186 m_err.fast( 1, 1 ) += m11 * s1s; // error(x,x)
187
188 m_err.fast( 2, 1 ) += m11 * c12; // error(y,x)
189 m_err.fast( 2, 2 ) += m11 * s2s; // error(y,y)
190
191 m_err.fast( 3, 1 ) += m11 * c13; // error(z,x)
192 m_err.fast( 3, 2 ) += m11 * c23; // error(z,y)
193 m_err.fast( 3, 3 ) += m11 * s3s; // error(z,z)
194
195 m_err.fast( 4, 1 ) += m12 * s1s; // error(px,x)
196 m_err.fast( 4, 2 ) += m12 * c12; // error(px,y)
197 m_err.fast( 4, 3 ) += m12 * c13; // error(px,z)
198 m_err.fast( 4, 4 ) += m22 * s1s; // error(px,px)
199
200 m_err.fast( 5, 1 ) += m12 * c12; // error(py,x)
201 m_err.fast( 5, 2 ) += m12 * s2s; // error(py,y)
202 m_err.fast( 5, 3 ) += m12 * c23; // error(py,z)
203 m_err.fast( 5, 4 ) += m22 * c12; // error(py,px)
204 m_err.fast( 5, 5 ) += m22 * s2s; // error(py,py)
205
206 m_err.fast( 6, 1 ) += m12 * c13; // error(pz,x)
207 m_err.fast( 6, 2 ) += m12 * c23; // error(pz,y)
208 m_err.fast( 6, 3 ) += m12 * s3s; // error(pz,z)
209 m_err.fast( 6, 4 ) += m22 * c13; // error(pz,px)
210 m_err.fast( 6, 5 ) += m22 * c23; // error(pz,py)
211 m_err.fast( 6, 6 ) += m22 * s3s; // error(pz,pz)
212 }
213 }
214
215 // Now store the new xv and pv.
216 m_xv = xv1;
217 m_pv = pv1;
218 return 1;
219}
220
221// Assign "=" operator
222Ext_xp_err& Ext_xp_err::operator=( const Ext_xp_err& err ) {
223 if ( this != &err )
224 {
225 *( (Ext_errmx*)this ) = err;
226 m_xv = err.m_xv;
227 m_pv = err.m_pv;
228 m_xp_jcb = err.m_xp_jcb;
229 m_h2xp_jcb = err.m_h2xp_jcb;
230 m_q = err.m_q;
231 m_mass2 = err.m_mass2;
232 }
233 return *this;
234}
235
236/*
237 ostream "<<" operator
238*/
239std::ostream& operator<<( std::ostream& s, const Ext_xp_err& err ) {
240 s << "m_err: " << (Ext_errmx)err << std::endl
241 << "m_xv: " << err.m_xv << std::endl
242 << " abs(xv)= " << err.m_xv.mag() << std::endl
243 << "m_pv: " << err.m_pv << std::endl
244 << " abs(pv)= " << err.m_pv.mag() << std::endl
245 << "m_xp_jcb: " << err.m_xp_jcb << std::endl
246 << "m_h2xp_jcb: " << err.m_h2xp_jcb << std::endl
247 << "m_q: " << err.m_q << std::endl
248 << "m_mass: " << sqrt( err.m_mass2 ) << std::endl;
249 return s;
250}
p2
double p2[4]
Definition BabayagaNLO.cxx:41
mass
double mass
Definition CosmicGenerator.cxx:128
operator<<
std::ostream & operator<<(std::ostream &s, const Ext_xp_err &err)
Definition Ext_xp_err.cxx:239
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
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::m_err
HepSymMatrix m_err
Definition Ext_errmx.h:128
Ext_errmx::put_err
void put_err(const double error[])
Definition Ext_errmx.cxx:47
Ext_errmx::Ext_errmx
Ext_errmx()
Definition Ext_errmx.cxx:28
Ext_xp_err::operator=
Ext_xp_err & operator=(const Ext_xp_err &xp_err)
Definition Ext_xp_err.cxx:222
Ext_xp_err::Ext_xp_err
Ext_xp_err()
Definition Ext_xp_err.cxx:32
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:53
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:77