1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
|
/* { dg-do compile } */
/* { dg-options "-std=c++11 -w -O2 -fPIC" } */
namespace CLHEP {
static const double meter = 1000.*10;
static const double meter2 = meter*meter;
static const double megaelectronvolt = 1. ;
static const double gigaelectronvolt = 1.e+3;
static const double GeV = gigaelectronvolt;
static const double megavolt = megaelectronvolt;
static const double volt = 1.e-6*megavolt;
static const double tesla = volt*1.e+9/meter2;
}
using CLHEP::GeV;
using CLHEP::tesla;
namespace std {
typedef long int ptrdiff_t;
}
extern "C" {
extern double cos (double __x) throw ();
extern double sin (double __x) throw ();
extern double sqrt (double __x) throw ();
}
namespace std __attribute__ ((__visibility__ ("default"))) {
using ::cos;
using ::sin;
using ::sqrt;
template<class _CharT> struct char_traits;
template<typename _CharT, typename _Traits = char_traits<_CharT> > struct basic_ostream;
typedef basic_ostream<char> ostream;
template<typename _Iterator> struct iterator_traits { };
template<typename _Tp> struct iterator_traits<_Tp*> {
typedef ptrdiff_t difference_type;
typedef _Tp& reference;
};
}
namespace __gnu_cxx __attribute__ ((__visibility__ ("default"))) {
using std::iterator_traits;
template<typename _Iterator, typename _Container> struct __normal_iterator {
_Iterator _M_current;
typedef iterator_traits<_Iterator> __traits_type;
typedef typename __traits_type::difference_type difference_type;
typedef typename __traits_type::reference reference;
explicit __normal_iterator(const _Iterator& __i) : _M_current(__i) { }
reference operator*() const {
return *_M_current;
}
__normal_iterator operator+(difference_type __n) const {
return __normal_iterator(_M_current + __n);
}
};
template<typename _Tp> struct new_allocator {
};
}
namespace std __attribute__ ((__visibility__ ("default"))) {
template<typename _Tp> struct allocator: public __gnu_cxx::new_allocator<_Tp> {
};
struct ios_base { };
template<typename _CharT, typename _Traits> struct basic_ios : public ios_base { };
template<typename _CharT, typename _Traits> struct basic_ostream : virtual public basic_ios<_CharT, _Traits> {
typedef basic_ostream<_CharT, _Traits> __ostream_type;
__ostream_type& operator<<(__ostream_type& (*__pf)(__ostream_type&)) { }
__ostream_type& operator<<(const void* __p) {
return _M_insert(__p);
}
template<typename _ValueT> __ostream_type& _M_insert(_ValueT __v);
};
template<typename _CharT, typename _Traits> inline basic_ostream<_CharT, _Traits>& endl(basic_ostream<_CharT, _Traits>& __os) {
}
}
typedef double G4double;
typedef int G4int;
extern __thread std::ostream *G4cout_p;
struct G4Field;
struct G4FieldManager {
inline G4Field* GetDetectorField() ;
};
namespace CLHEP {
struct Hep3Vector {
Hep3Vector(double x, double y, double z);
inline ~Hep3Vector();
inline double x() const;
inline double y() const;
inline double z() const;
inline double mag() const;
inline Hep3Vector cross(const Hep3Vector &) const;
double dx;
double dy;
double dz;
};
Hep3Vector operator / (const Hep3Vector &, double a);
inline double Hep3Vector::x() const {
return dx;
}
inline double Hep3Vector::y() const {
return dy;
}
inline double Hep3Vector::z() const {
return dz;
}
inline Hep3Vector operator + (const Hep3Vector & a, const Hep3Vector & b) { }
inline Hep3Vector operator * (const Hep3Vector & p, double a) { }
inline double operator * (const Hep3Vector & a, const Hep3Vector & b) { }
inline Hep3Vector::Hep3Vector(double x1, double y1, double z1) : dx(x1), dy(y1), dz(z1) {
}
inline Hep3Vector::~Hep3Vector() { }
inline Hep3Vector Hep3Vector::cross(const Hep3Vector & p) const {
return Hep3Vector(dy*p.dz-p.dy*dz, dz*p.dx-p.dz*dx, dx*p.dy-p.dx*dy);
}
}
typedef CLHEP::Hep3Vector G4ThreeVector;
namespace std __attribute__ ((__visibility__ ("default"))) {
template<typename _Tp, typename _Alloc = std::allocator<_Tp> > struct vector
{
typedef _Tp *pointer;
typedef __gnu_cxx::__normal_iterator<pointer, vector> iterator;
iterator begin() { }
};
}
struct G4TransportationManager {
static G4TransportationManager* GetTransportationManager();
inline G4FieldManager* GetFieldManager() const;
};
struct G4ErrorMatrix {
G4ErrorMatrix(G4int p, G4int q, G4int i);
virtual ~G4ErrorMatrix();
struct G4ErrorMatrix_row {
inline G4ErrorMatrix_row(G4ErrorMatrix&,G4int);
G4double & operator[](G4int);
G4ErrorMatrix& _a;
G4int _r;
};
inline G4ErrorMatrix_row operator[] (G4int);
std::vector<G4double > m;
G4int nrow, ncol;
};
inline G4ErrorMatrix::G4ErrorMatrix_row G4ErrorMatrix::operator[] (G4int r) {
G4ErrorMatrix_row b(*this,r);
return b;
}
inline G4double &G4ErrorMatrix::G4ErrorMatrix_row::operator[](G4int c) {
return *(_a.m.begin()+_r*_a.ncol+c);
}
inline G4ErrorMatrix:: G4ErrorMatrix_row::G4ErrorMatrix_row(G4ErrorMatrix&a, G4int r) : _a(a) {
_r = r;
};
struct G4DynamicParticle {
G4double GetCharge() const;
};
struct G4Step;
struct G4Track {
const G4DynamicParticle* GetDynamicParticle() const;
const G4ThreeVector& GetPosition() const;
G4ThreeVector GetMomentum() const;
const G4Step* GetStep() const;
};
struct G4StepPoint {
const G4ThreeVector& GetPosition() const;
G4ThreeVector GetMomentum() const;
};
struct G4Step {
G4StepPoint* GetPreStepPoint() const;
G4double GetStepLength() const;
};
namespace HepGeom {
template<class T> struct BasicVector3D {
T v_[3];
BasicVector3D(T x1, T y1, T z1) { }
operator T * () {
return v_;
}
T x() const {
return v_[0];
}
T y() const {
return v_[1];
}
T z() const {
return v_[2];
}
T perp2() const { }
T perp() const {
return std::sqrt(perp2());
}
T mag2() const { }
T mag() const {
return std::sqrt(mag2());
}
T theta() const { }
};
inline BasicVector3D<double> operator-(const BasicVector3D<double> & a,const BasicVector3D<double> & b) { }
inline BasicVector3D<double> operator*(const BasicVector3D<double> & v, double a) { }
template<class T> struct Point3D : public BasicVector3D<T> {
explicit Point3D(const double * a) : BasicVector3D<double>(a[0],a[1],a[2]) { }
Point3D(const CLHEP::Hep3Vector & v) : BasicVector3D<double>(v.dx,v.dy,v.dz) { }
};
}
typedef HepGeom::Point3D<G4double> G4Point3D;
namespace HepGeom {
template<class T> struct Vector3D : public BasicVector3D<T> {
Vector3D(const BasicVector3D<double> & v) : BasicVector3D<double>(v) { }
Vector3D(const CLHEP::Hep3Vector & v) : BasicVector3D<double>(v.dx,v.dy,v.dz) { }
operator CLHEP::Hep3Vector () const { }
};
}
typedef HepGeom::Vector3D<G4double> G4Vector3D;
struct G4ErrorFreeTrajState
{
virtual G4int PropagateError( const G4Track* aTrack );
G4int PropagateErrorMSC( const G4Track* aTrack );
};
G4int G4ErrorFreeTrajState::PropagateError( const G4Track* aTrack ) {
G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/10.;
G4Point3D vposPost = aTrack->GetPosition()/10.;
G4Vector3D vpPost = aTrack->GetMomentum()/GeV;
G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition()/10.;
G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV;
G4double pPre = vpPre.mag();
G4double pPost = vpPost.mag();
G4double pInvPre = 1./pPre;
G4double pInvPost = 1./pPost;
G4double deltaPInv = pInvPost - pInvPre;
G4Vector3D vpPreNorm = vpPre * pInvPre;
G4Vector3D vpPostNorm = vpPost * pInvPost;
(*G4cout_p) << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm << std::endl;
G4double sinpPre = std::sin( vpPreNorm.theta() );
G4double sinpPostInv = 1./std::sin( vpPreNorm.theta() );
G4ErrorMatrix transf(5, 5, 0 );
G4double charge = aTrack->GetDynamicParticle()->GetCharge();
G4double h1[3], h2[3];
G4Field* field
= G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField()
;
if( charge != 0. && field )
{
G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.;
G4ThreeVector HPost= G4ThreeVector( h2[0], h2[1], h2[2] ) / tesla *10.;
{
G4double pInvAver = 1./(pInvPre + pInvPost );
G4double CFACT8 = 2.997925E-4;
G4ThreeVector vHAverNorm( (HPre*pInvPre + HPost*pInvPost ) * pInvAver * charge * CFACT8 );
G4double HAver = vHAverNorm.mag();
G4double pAver = (pPre+pPost)*0.5;
G4double QAver = -HAver/pAver;
G4double thetaAver = QAver * stepLengthCm;
G4double sinThetaAver = std::sin(thetaAver);
G4double cosThetaAver = std::cos(thetaAver);
G4double gamma = vHAverNorm * vpPostNorm;
G4ThreeVector AN2 = vHAverNorm.cross( vpPostNorm );
G4double AU = 1./vpPreNorm.perp();
G4ThreeVector vUPre( -AU*vpPreNorm.y(), AU*vpPreNorm.x(), 0. );
G4ThreeVector vVPre( -vpPreNorm.z()*vUPre.y(), vpPreNorm.z()*vUPre.x(), vpPreNorm.x()*vUPre.y() - vpPreNorm.y()*vUPre.x() );
AU = 1./vpPostNorm.perp();
G4ThreeVector vUPost( -AU*vpPostNorm.y(), AU*vpPostNorm.x(), 0. );
G4ThreeVector vVPost( -vpPostNorm.z()*vUPost.y(), vpPostNorm.z()*vUPost.x(), vpPostNorm.x()*vUPost.y() - vpPostNorm.y()*vUPost.x() );
G4Point3D deltaPos( vposPre - vposPost );
G4double QP = QAver * pAver;
G4double ANV = -( vHAverNorm.x()*vUPost.x() + vHAverNorm.y()*vUPost.y() );
G4double ANU = ( vHAverNorm.x()*vVPost.x() + vHAverNorm.y()*vVPost.y() + vHAverNorm.z()*vVPost.z() );
G4double OMcosThetaAver = 1. - cosThetaAver;
G4double TMSINT = thetaAver - sinThetaAver;
G4ThreeVector vHUPre( -vHAverNorm.z() * vUPre.y(), vHAverNorm.z() * vUPre.x(), vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x() );
G4ThreeVector vHVPre( vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(), vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(), vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x() );
transf[0][1] = -deltaPInv/thetaAver* ( TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) + sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + OMcosThetaAver*(vHVPre.x()*vpPostNorm.x()+vHVPre.y()*vpPostNorm.y()+vHVPre.z()*vpPostNorm.z()) );
transf[0][2] = -sinpPre*deltaPInv/thetaAver* ( TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) + sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + OMcosThetaAver*(vHUPre.x()*vpPostNorm.x()+vHUPre.y()*vpPostNorm.y()+vHUPre.z()*vpPostNorm.z()) );
transf[0][3] = -deltaPInv/stepLengthCm*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() );
transf[1][1] = cosThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) + sinThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) + OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())* (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) + ANV*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) - TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) );
transf[1][2] = cosThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) + sinThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) + OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )* (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) + ANV*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) - TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) );
transf[2][0] = -QP*ANU*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())*sinpPostInv *(1.+deltaPInv*pAver);
transf[2][3] = -QAver*ANU*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() )*sinpPostInv;
transf[3][4] = (vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() );
transf[4][0] = pAver*(vVPost.x()*deltaPos.x()+vVPost.y()*deltaPos.y()+vVPost.z()*deltaPos.z()) *(1.+deltaPInv*pAver);
transf[4][1] = ( sinThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) + OMcosThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) + TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())* (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver;
transf[4][2] = ( sinThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) + OMcosThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) + TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())* (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver;
}
}
PropagateErrorMSC( aTrack );
}
|
