1 /* { dg-do compile } */ 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 } 278