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