1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // $Id: G4ErrorFreeTrajState.cc,v 1.7 2007/09/24 16:25:57 arce Exp $ |
---|
27 | // GEANT4 tag $Name: $ |
---|
28 | // |
---|
29 | // ------------------------------------------------------------ |
---|
30 | // GEANT 4 class implementation file |
---|
31 | // ------------------------------------------------------------ |
---|
32 | // |
---|
33 | #include "G4ErrorFreeTrajState.hh" |
---|
34 | #include "G4ErrorFreeTrajParam.hh" |
---|
35 | #include "G4ErrorSurfaceTrajState.hh" |
---|
36 | |
---|
37 | #include "G4ErrorMatrix.hh" |
---|
38 | #include <iomanip> |
---|
39 | |
---|
40 | #include "G4Field.hh" |
---|
41 | #include "G4FieldManager.hh" |
---|
42 | #include "G4TransportationManager.hh" |
---|
43 | #include "G4GeometryTolerance.hh" |
---|
44 | #include "G4Material.hh" |
---|
45 | #include "G4ErrorPropagatorData.hh" |
---|
46 | |
---|
47 | //------------------------------------------------------------------------ |
---|
48 | G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4String& partType, const G4Point3D& pos, const G4Vector3D& mom, const G4ErrorTrajErr& errmat) : G4ErrorTrajState( partType, pos, mom, errmat ) |
---|
49 | { |
---|
50 | fTrajParam = G4ErrorFreeTrajParam( pos, mom ); |
---|
51 | Init(); |
---|
52 | } |
---|
53 | |
---|
54 | |
---|
55 | //------------------------------------------------------------------------ |
---|
56 | G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4ErrorSurfaceTrajState& tpSD ) : G4ErrorTrajState( tpSD.GetParticleType(), tpSD.GetPosition(), tpSD.GetMomentum() ) |
---|
57 | { |
---|
58 | // G4ThreeVector planeNormal = tpSD.GetPlaneNormal(); |
---|
59 | // G4double fPt = tpSD.GetMomentum()*planeNormal;//mom projected on normal to plane |
---|
60 | // G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters(); |
---|
61 | // G4ThreeVector Psc = fPt * planeNormal + tpSDparam.GetPU()*tpSDparam.GetVectorU() + tpSD.GetPV()*tpSD.GetVectorW(); |
---|
62 | |
---|
63 | fTrajParam = G4ErrorFreeTrajParam( fPosition, fMomentum ); |
---|
64 | Init(); |
---|
65 | |
---|
66 | //----- Get the error matrix in SC coordinates |
---|
67 | G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters(); |
---|
68 | G4double mom = fMomentum.mag(); |
---|
69 | G4double mom2 = fMomentum.mag2(); |
---|
70 | G4double TVW1 = std::sqrt( mom2 / ( mom2 + tpSDparam.GetPV()*tpSDparam.GetPV() + tpSDparam.GetPV()*tpSDparam.GetPV()) ); |
---|
71 | G4ThreeVector vTVW( TVW1, tpSDparam.GetPV()/mom * TVW1, tpSDparam.GetPW()/mom * TVW1 ); |
---|
72 | G4Vector3D vectorU = tpSDparam.GetVectorV().cross( tpSDparam.GetVectorW() ); |
---|
73 | G4Vector3D vTN = vTVW.x()*vectorU + vTVW.y()*tpSDparam.GetVectorV() + vTVW.z()*tpSDparam.GetVectorW(); |
---|
74 | |
---|
75 | #ifdef G4EVERBOSE |
---|
76 | if( iverbose >= 5){ |
---|
77 | G4double pc2 = std::asin( vTN.z() ); |
---|
78 | G4double pc3 = std::atan (vTN.y()/vTN.x()); |
---|
79 | |
---|
80 | G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda() << " diff " << pc2-GetParameters().GetLambda() << G4endl; |
---|
81 | G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi() << " diff " << pc3-GetParameters().GetPhi() << G4endl; |
---|
82 | } |
---|
83 | #endif |
---|
84 | |
---|
85 | //--- Get the unit vectors perp to P |
---|
86 | G4double cosl = std::cos( GetParameters().GetLambda() ); |
---|
87 | if (cosl < 1.E-30) cosl = 1.E-30; |
---|
88 | G4double cosl1 = 1./cosl; |
---|
89 | G4Vector3D vUN(-vTN.y()*cosl1, vTN.x()*cosl1, 0. ); |
---|
90 | G4Vector3D vVN(-vTN.z()*vUN.y(), vTN.z()*vUN.x(), cosl ); |
---|
91 | |
---|
92 | G4Vector3D vUperp = G4Vector3D( -fMomentum.y(), fMomentum.x(), 0.); |
---|
93 | G4Vector3D vVperp = vUperp.cross( fMomentum ); |
---|
94 | vUperp *= 1./vUperp.mag(); |
---|
95 | vVperp *= 1./vVperp.mag(); |
---|
96 | |
---|
97 | #ifdef G4EVERBOSE |
---|
98 | if( iverbose >= 5){ |
---|
99 | G4cout << " CHECK: vUN " << vUN << " = " << vUperp << " diff " << (vUN-vUperp).mag() << G4endl; |
---|
100 | G4cout << " CHECK: vVN " << vVN << " = " << vVperp << " diff " << (vVN-vVperp).mag() << G4endl; |
---|
101 | } |
---|
102 | #endif |
---|
103 | |
---|
104 | //get the dot products of vectors perpendicular to direction and vector defining SD plane |
---|
105 | G4double dUU = vUperp * tpSD.GetVectorV(); |
---|
106 | G4double dUV = vUperp * tpSD.GetVectorW(); |
---|
107 | G4double dVU = vVperp * tpSD.GetVectorV(); |
---|
108 | G4double dVV = vVperp * tpSD.GetVectorW(); |
---|
109 | |
---|
110 | |
---|
111 | //--- Get transformation first |
---|
112 | G4ErrorMatrix transfM(5, 5, 1 ); |
---|
113 | //--- Get magnetic field |
---|
114 | const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField(); |
---|
115 | G4ThreeVector dir = fTrajParam.GetDirection(); |
---|
116 | G4double invCosTheta = 1./std::cos( dir.theta() ); |
---|
117 | |
---|
118 | if( fCharge != 0 |
---|
119 | && field ) { |
---|
120 | G4double pos1[3]; pos1[0] = fPosition.x()*cm; pos1[1] = fPosition.y()*cm; pos1[2] = fPosition.z()*cm; |
---|
121 | G4double h1[3]; |
---|
122 | field->GetFieldValue( pos1, h1 ); |
---|
123 | G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.; |
---|
124 | G4double magHPre = HPre.mag(); |
---|
125 | G4double invP = 1./fMomentum.mag(); |
---|
126 | G4double magHPreM = magHPre * invP; |
---|
127 | if( magHPre != 0. ) { |
---|
128 | G4double magHPreM2 = fCharge / magHPre; |
---|
129 | |
---|
130 | G4double Q = -magHPreM * c_light; |
---|
131 | G4double sinz = -HPre*vUperp * magHPreM2; |
---|
132 | G4double cosz = HPre*vVperp * magHPreM2; |
---|
133 | |
---|
134 | transfM[1][3] = -Q*dir.y()*sinz; |
---|
135 | transfM[1][4] = -Q*dir.z()*sinz; |
---|
136 | transfM[2][3] = -Q*dir.y()*cosz*invCosTheta; |
---|
137 | transfM[2][4] = -Q*dir.z()*cosz*invCosTheta; |
---|
138 | } |
---|
139 | } |
---|
140 | |
---|
141 | transfM[0][0] = 1.; |
---|
142 | transfM[1][1] = dir.x()*dVU; |
---|
143 | transfM[1][2] = dir.x()*dVV; |
---|
144 | transfM[2][1] = dir.x()*dUU*invCosTheta; |
---|
145 | transfM[2][2] = dir.x()*dUV*invCosTheta; |
---|
146 | transfM[3][3] = dUU; |
---|
147 | transfM[3][4] = dUV; |
---|
148 | transfM[4][3] = dVU; |
---|
149 | transfM[4][4] = dVV; |
---|
150 | |
---|
151 | fError = G4ErrorTrajErr( tpSD.GetError().similarity( transfM ) ); |
---|
152 | |
---|
153 | #ifdef G4EVERBOSE |
---|
154 | if( iverbose >= 1) G4cout << "error matrix SD2SC " << fError << G4endl; |
---|
155 | if( iverbose >= 4) G4cout << "G4ErrorFreeTrajState from SD " << *this << G4endl; |
---|
156 | #endif |
---|
157 | } |
---|
158 | |
---|
159 | |
---|
160 | //------------------------------------------------------------------------ |
---|
161 | void G4ErrorFreeTrajState::Init() |
---|
162 | { |
---|
163 | theTSType = G4eTS_FREE; |
---|
164 | BuildCharge(); |
---|
165 | theTransfMat = G4ErrorMatrix(5,5,0); |
---|
166 | //- theFirstStep = true; |
---|
167 | } |
---|
168 | |
---|
169 | //------------------------------------------------------------------------ |
---|
170 | void G4ErrorFreeTrajState::Dump( std::ostream& out ) const |
---|
171 | { |
---|
172 | out << *this; |
---|
173 | } |
---|
174 | |
---|
175 | //------------------------------------------------------------------------ |
---|
176 | G4int G4ErrorFreeTrajState::Update( const G4Track* aTrack ) |
---|
177 | { |
---|
178 | G4int ierr = 0; |
---|
179 | fTrajParam.Update( aTrack ); |
---|
180 | UpdatePosMom( aTrack->GetPosition(), aTrack->GetMomentum() ); |
---|
181 | return ierr; |
---|
182 | |
---|
183 | } |
---|
184 | |
---|
185 | |
---|
186 | //------------------------------------------------------------------------ |
---|
187 | std::ostream& operator<<(std::ostream& out, const G4ErrorFreeTrajState& ts) |
---|
188 | { |
---|
189 | out.setf(std::ios::fixed,std::ios::floatfield); |
---|
190 | |
---|
191 | |
---|
192 | ts.DumpPosMomError( out ); |
---|
193 | |
---|
194 | out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl; |
---|
195 | |
---|
196 | return out; |
---|
197 | |
---|
198 | } |
---|
199 | |
---|
200 | |
---|
201 | //------------------------------------------------------------------------ |
---|
202 | G4int G4ErrorFreeTrajState::PropagateError( const G4Track* aTrack ) |
---|
203 | { |
---|
204 | G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm; |
---|
205 | G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); |
---|
206 | |
---|
207 | if( std::fabs(stepLengthCm) <= kCarTolerance/cm ) return 0; |
---|
208 | |
---|
209 | #ifdef G4EVERBOSE |
---|
210 | if( iverbose >= 2 )G4cout << " G4ErrorFreeTrajState::PropagateError " << G4endl; |
---|
211 | #endif |
---|
212 | |
---|
213 | // * *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES |
---|
214 | G4Point3D vposPost = aTrack->GetPosition()/cm; |
---|
215 | G4Vector3D vpPost = aTrack->GetMomentum()/GeV; |
---|
216 | // G4Point3D vposPre = fPosition/cm; |
---|
217 | // G4Vector3D vpPre = fMomentum/GeV; |
---|
218 | G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition()/cm; |
---|
219 | G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV; |
---|
220 | //correct to avoid propagation along Z |
---|
221 | if( vpPre.mag() == vpPre.z() ) vpPre.setX( 1.E-6*MeV ); |
---|
222 | if( vpPost.mag() == vpPost.z() ) vpPost.setX( 1.E-6*MeV ); |
---|
223 | |
---|
224 | G4double pPre = vpPre.mag(); |
---|
225 | G4double pPost = vpPost.mag(); |
---|
226 | #ifdef G4EVERBOSE |
---|
227 | if( iverbose >= 2 ) { |
---|
228 | G4cout << "G4EP: vposPre " << vposPre << G4endl |
---|
229 | << "G4EP: vposPost " << vposPost << G4endl; |
---|
230 | G4cout << "G4EP: vpPre " << vpPre << G4endl |
---|
231 | << "G4EP: vpPost " << vpPost << G4endl; |
---|
232 | G4cout << " err start step " << fError << G4endl; |
---|
233 | G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl; |
---|
234 | } |
---|
235 | #endif |
---|
236 | |
---|
237 | if( pPre == 0. || pPost == 0 ) return 2; |
---|
238 | G4double pInvPre = 1./pPre; |
---|
239 | G4double pInvPost = 1./pPost; |
---|
240 | G4double deltaPInv = pInvPost - pInvPre; |
---|
241 | |
---|
242 | G4Vector3D vpPreNorm = vpPre * pInvPre; |
---|
243 | G4Vector3D vpPostNorm = vpPost * pInvPost; |
---|
244 | // if( iverbose >= 2 ) G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm << G4endl; |
---|
245 | //return if propagation along Z?? |
---|
246 | if( 1. - std::fabs(vpPostNorm.z()) < kCarTolerance ) return 4; |
---|
247 | G4double sinpPre = std::sin( vpPreNorm.theta() ); //cosine perpendicular to pPre = sine pPre |
---|
248 | G4double sinpPost = std::sin( vpPostNorm.theta() ); //cosine perpendicular to pPost = sine pPost |
---|
249 | G4double sinpPostInv = 1./std::sin( vpPreNorm.theta() ); |
---|
250 | |
---|
251 | #ifdef G4EVERBOSE |
---|
252 | if( iverbose >= 2 ) G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl; |
---|
253 | #endif |
---|
254 | //* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR |
---|
255 | //* *** NEUTRAL PARTICLE OR FIELDFREE REGION |
---|
256 | G4ErrorMatrix transf(5, 5, 0 ); |
---|
257 | |
---|
258 | transf[3][2] = stepLengthCm * sinpPost; |
---|
259 | transf[4][1] = stepLengthCm; |
---|
260 | for( size_t ii=0;ii < 5; ii++ ){ |
---|
261 | transf[ii][ii] = 1.; |
---|
262 | } |
---|
263 | #ifdef G4EVERBOSE |
---|
264 | if( iverbose >= 2 ) { |
---|
265 | G4cout << "G4EP: transf matrix neutral " << transf; |
---|
266 | } |
---|
267 | #endif |
---|
268 | |
---|
269 | // charge X propagation direction |
---|
270 | G4double charge = aTrack->GetDynamicParticle()->GetCharge(); |
---|
271 | if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ) { |
---|
272 | charge *= -1.; |
---|
273 | } |
---|
274 | // G4cout << " charge " << charge << G4endl; |
---|
275 | //t check if particle has charge |
---|
276 | //t if( charge == 0 ) goto 45; |
---|
277 | // check if the magnetic field is = 0. |
---|
278 | |
---|
279 | //position is from geant4, it is assumed to be in mm (for debugging, eventually it will not be transformed) |
---|
280 | G4double pos1[3]; pos1[0] = vposPre.x()*cm; pos1[1] = vposPre.y()*cm; pos1[2] = vposPre.z()*cm; |
---|
281 | G4double pos2[3]; pos2[0] = vposPost.x()*cm; pos2[1] = vposPost.y()*cm; pos2[2] = vposPost.z()*cm; |
---|
282 | G4double h1[3], h2[3]; |
---|
283 | |
---|
284 | const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField(); |
---|
285 | if( !field ) return 0; //goto 45 |
---|
286 | |
---|
287 | // calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION |
---|
288 | if( charge != 0. && field ) { |
---|
289 | |
---|
290 | field->GetFieldValue( pos1, h1 ); |
---|
291 | field->GetFieldValue( pos2, h2 ); |
---|
292 | G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.; //10. is to get same dimensions as GEANT3 (kilogauss) |
---|
293 | G4ThreeVector HPost= G4ThreeVector( h2[0], h2[1], h2[2] ) / tesla *10.; |
---|
294 | G4double magHPre = HPre.mag(); |
---|
295 | G4double magHPost = HPost.mag(); |
---|
296 | #ifdef G4EVERBOSE |
---|
297 | if( iverbose >= 2 ) G4cout << "G4EP: HPre " << HPre << G4endl |
---|
298 | << "G4EP: HPost " << HPost << G4endl; |
---|
299 | #endif |
---|
300 | |
---|
301 | if( magHPre + magHPost != 0. ) { |
---|
302 | |
---|
303 | //* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2 |
---|
304 | G4double gam; |
---|
305 | if( magHPost != 0. ){ |
---|
306 | gam = HPost * vpPostNorm / magHPost; |
---|
307 | }else { |
---|
308 | gam = HPre * vpPreNorm / magHPre; |
---|
309 | } |
---|
310 | |
---|
311 | // G4eMagneticLimitsProcess will limit the step, but based on an straight line trajectory |
---|
312 | G4double alphaSqr = 1. - gam * gam; |
---|
313 | G4double diffHSqr = ( HPre * pInvPre - HPost * pInvPost ).mag2(); |
---|
314 | G4double delhp6Sqr = 300.*300.; |
---|
315 | #ifdef G4EVERBOSE |
---|
316 | if( iverbose >= 2 ) G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr << " diffHSqr " << diffHSqr << G4endl; |
---|
317 | #endif |
---|
318 | if( diffHSqr * alphaSqr > delhp6Sqr ) return 3; |
---|
319 | |
---|
320 | |
---|
321 | //* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT |
---|
322 | G4double pInvAver = 1./(pInvPre + pInvPost ); |
---|
323 | G4double CFACT8 = 2.997925E-4; |
---|
324 | //G4double HAver |
---|
325 | G4ThreeVector vHAverNorm( (HPre*pInvPre + HPost*pInvPost ) * pInvAver * charge * CFACT8 ); |
---|
326 | G4double HAver = vHAverNorm.mag(); |
---|
327 | G4double invHAver = 1./HAver; |
---|
328 | vHAverNorm *= invHAver; |
---|
329 | #ifdef G4EVERBOSE |
---|
330 | if( iverbose >= 2 ) G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver << " charge " << charge<< G4endl; |
---|
331 | #endif |
---|
332 | |
---|
333 | G4double pAver = (pPre+pPost)*0.5; |
---|
334 | G4double QAver = -HAver/pAver; |
---|
335 | G4double thetaAver = QAver * stepLengthCm; |
---|
336 | G4double sinThetaAver = std::sin(thetaAver); |
---|
337 | G4double cosThetaAver = std::cos(thetaAver); |
---|
338 | G4double gamma = vHAverNorm * vpPostNorm; |
---|
339 | G4ThreeVector AN2 = vHAverNorm.cross( vpPostNorm ); |
---|
340 | |
---|
341 | #ifdef G4EVERBOSE |
---|
342 | if( iverbose >= 2 ) G4cout << " G4EP: AN2 " << AN2 << G4endl; |
---|
343 | #endif |
---|
344 | G4double AU = 1./vpPreNorm.perp(); |
---|
345 | //t G4ThreeVector vU( vpPreNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU ); |
---|
346 | G4ThreeVector vUPre( -AU*vpPreNorm.y(), |
---|
347 | AU*vpPreNorm.x(), |
---|
348 | 0. ); |
---|
349 | G4ThreeVector vVPre( -vpPreNorm.z()*vUPre.y(), |
---|
350 | vpPreNorm.z()*vUPre.x(), |
---|
351 | vpPreNorm.x()*vUPre.y() - vpPreNorm.y()*vUPre.x() ); |
---|
352 | |
---|
353 | // |
---|
354 | AU = 1./vpPostNorm.perp(); |
---|
355 | //t G4ThreeVector vU( vpPostNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU ); |
---|
356 | G4ThreeVector vUPost( -AU*vpPostNorm.y(), |
---|
357 | AU*vpPostNorm.x(), |
---|
358 | 0. ); |
---|
359 | G4ThreeVector vVPost( -vpPostNorm.z()*vUPost.y(), |
---|
360 | vpPostNorm.z()*vUPost.x(), |
---|
361 | vpPostNorm.x()*vUPost.y() - vpPostNorm.y()*vUPost.x() ); |
---|
362 | #ifdef G4EVERBOSE |
---|
363 | //- G4cout << " vpPostNorm " << vpPostNorm << G4endl; |
---|
364 | if( iverbose >= 2 ) G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre << " vUPost " << vUPost << " vVPost " << vVPost << G4endl; |
---|
365 | #endif |
---|
366 | G4Point3D deltaPos( vposPre - vposPost ); |
---|
367 | |
---|
368 | // * *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2 |
---|
369 | // * *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT |
---|
370 | // * *** TAKEN INTO ACCOUNT |
---|
371 | |
---|
372 | G4double QP = QAver * pAver; // = -HAver |
---|
373 | #ifdef G4EVERBOSE |
---|
374 | if( iverbose >= 2) G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver << G4endl; |
---|
375 | #endif |
---|
376 | G4double ANV = -( vHAverNorm.x()*vUPost.x() + vHAverNorm.y()*vUPost.y() ); |
---|
377 | G4double ANU = ( vHAverNorm.x()*vVPost.x() + vHAverNorm.y()*vVPost.y() + vHAverNorm.z()*vVPost.z() ); |
---|
378 | G4double OMcosThetaAver = 1. - cosThetaAver; |
---|
379 | #ifdef G4EVERBOSE |
---|
380 | if( iverbose >= 2) G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver " << cosThetaAver << " thetaAver " << thetaAver << " QAver " << QAver << " stepLengthCm " << stepLengthCm << G4endl; |
---|
381 | #endif |
---|
382 | G4double TMSINT = thetaAver - sinThetaAver; |
---|
383 | #ifdef G4EVERBOSE |
---|
384 | if( iverbose >= 2 ) G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl; |
---|
385 | #endif |
---|
386 | |
---|
387 | G4ThreeVector vHUPre( -vHAverNorm.z() * vUPre.y(), |
---|
388 | vHAverNorm.z() * vUPre.x(), |
---|
389 | vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x() ); |
---|
390 | #ifdef G4EVERBOSE |
---|
391 | // if( iverbose >= 2 ) G4cout << "G4EP: HUPre(1) " << vHUPre.x() << " " << vHAverNorm.z() << " " << vUPre.y() << G4endl; |
---|
392 | #endif |
---|
393 | G4ThreeVector vHVPre( vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(), |
---|
394 | vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(), |
---|
395 | vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x() ); |
---|
396 | #ifdef G4EVERBOSE |
---|
397 | if( iverbose >= 2 ) G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl; |
---|
398 | #endif |
---|
399 | |
---|
400 | //------------------- COMPUTE MATRIX |
---|
401 | //---------- 1/P |
---|
402 | |
---|
403 | transf[0][0] = 1.-deltaPInv*pAver*(1.+(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())/stepLengthCm) |
---|
404 | +2.*deltaPInv*pAver; |
---|
405 | |
---|
406 | transf[0][1] = -deltaPInv/thetaAver* |
---|
407 | ( TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) + |
---|
408 | sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + |
---|
409 | OMcosThetaAver*(vHVPre.x()*vpPostNorm.x()+vHVPre.y()*vpPostNorm.y()+vHVPre.z()*vpPostNorm.z()) ); |
---|
410 | |
---|
411 | transf[0][2] = -sinpPre*deltaPInv/thetaAver* |
---|
412 | ( TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) + |
---|
413 | sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + |
---|
414 | OMcosThetaAver*(vHUPre.x()*vpPostNorm.x()+vHUPre.y()*vpPostNorm.y()+vHUPre.z()*vpPostNorm.z()) ); |
---|
415 | |
---|
416 | transf[0][3] = -deltaPInv/stepLengthCm*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ); |
---|
417 | |
---|
418 | transf[0][4] = -deltaPInv/stepLengthCm*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()); |
---|
419 | |
---|
420 | // *** Lambda |
---|
421 | transf[1][0] = -QP*ANV*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z()) |
---|
422 | *(1.+deltaPInv*pAver); |
---|
423 | #ifdef G4EVERBOSE |
---|
424 | if(iverbose >= 3) G4cout << "ctransf10= " << transf[1][0] << " " << -QP<< " " << ANV<< " " << vpPostNorm.x()<< " " << deltaPos.x()<< " " << vpPostNorm.y()<< " " << deltaPos.y()<< " " << vpPostNorm.z()<< " " << deltaPos.z() |
---|
425 | << " " << deltaPInv<< " " << pAver << G4endl; |
---|
426 | #endif |
---|
427 | |
---|
428 | transf[1][1] = cosThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) + |
---|
429 | sinThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) + |
---|
430 | OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())* |
---|
431 | (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) + |
---|
432 | ANV*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + |
---|
433 | OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) - |
---|
434 | TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) ); |
---|
435 | |
---|
436 | transf[1][2] = cosThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) + |
---|
437 | sinThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) + |
---|
438 | OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )* |
---|
439 | (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) + |
---|
440 | ANV*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + |
---|
441 | OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) - |
---|
442 | TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) ); |
---|
443 | transf[1][2] = sinpPre*transf[1][2]; |
---|
444 | |
---|
445 | transf[1][3] = -QAver*ANV*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ); |
---|
446 | |
---|
447 | transf[1][4] = -QAver*ANV*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()); |
---|
448 | |
---|
449 | // *** Phi |
---|
450 | |
---|
451 | transf[2][0] = -QP*ANU*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())*sinpPostInv |
---|
452 | *(1.+deltaPInv*pAver); |
---|
453 | #ifdef G4EVERBOSE |
---|
454 | if(iverbose >= 3)G4cout <<"ctransf20= " << transf[2][0] <<" "<< -QP<<" "<<ANU<<" "<<vpPostNorm.x()<<" "<<deltaPos.x()<<" "<<vpPostNorm.y()<<" "<<deltaPos.y()<<" "<<vpPostNorm.z()<<" "<<deltaPos.z()<<" "<<sinpPostInv |
---|
455 | <<" "<<deltaPInv<<" "<<pAver<< G4endl; |
---|
456 | #endif |
---|
457 | transf[2][1] = cosThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) + |
---|
458 | sinThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) + |
---|
459 | OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())* |
---|
460 | (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) + |
---|
461 | ANU*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) + |
---|
462 | OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) - |
---|
463 | TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) ); |
---|
464 | transf[2][1] = sinpPostInv*transf[2][1]; |
---|
465 | |
---|
466 | transf[2][2] = cosThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) + |
---|
467 | sinThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) + |
---|
468 | OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )* |
---|
469 | (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) + |
---|
470 | ANU*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) + |
---|
471 | OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) - |
---|
472 | TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) ); |
---|
473 | transf[2][2] = sinpPostInv*sinpPre*transf[2][2]; |
---|
474 | |
---|
475 | transf[2][3] = -QAver*ANU*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() )*sinpPostInv; |
---|
476 | #ifdef G4EVERBOSE |
---|
477 | if(iverbose >= 3)G4cout <<"ctransf23= " << transf[2][3] <<" "<< -QAver<<" "<<ANU<<" "<<vUPre.x()<<" "<<vpPostNorm.x()<<" "<< vUPre.y()<<" "<<vpPostNorm.y()<<" "<<sinpPostInv<<G4endl; |
---|
478 | #endif |
---|
479 | |
---|
480 | transf[2][4] = -QAver*ANU*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z())*sinpPostInv; |
---|
481 | |
---|
482 | // *** Yt |
---|
483 | |
---|
484 | transf[3][0] = pAver*(vUPost.x()*deltaPos.x()+vUPost.y()*deltaPos.y() ) |
---|
485 | *(1.+deltaPInv*pAver); |
---|
486 | #ifdef G4EVERBOSE |
---|
487 | if(iverbose >= 3) G4cout <<"ctransf30= " << transf[3][0] <<" "<< pAver<<" "<<vUPost.x()<<" "<<deltaPos.x()<<" "<<vUPost.y()<<" "<<deltaPos.y() |
---|
488 | <<" "<<deltaPInv<<" "<<pAver<<G4endl; |
---|
489 | #endif |
---|
490 | |
---|
491 | transf[3][1] = ( sinThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) + |
---|
492 | OMcosThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) + |
---|
493 | TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )* |
---|
494 | (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver; |
---|
495 | |
---|
496 | transf[3][2] = ( sinThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) + |
---|
497 | OMcosThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) + |
---|
498 | TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )* |
---|
499 | (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver; |
---|
500 | #ifdef G4EVERBOSE |
---|
501 | if(iverbose >= 3) G4cout <<"ctransf32= " << transf[3][2] <<" "<< sinThetaAver<<" "<<vUPre.x()<<" "<<vUPost.x()<<" "<<vUPre.y()<<" "<<vUPost.y() <<" "<< |
---|
502 | OMcosThetaAver<<" "<<vHUPre.x()<<" "<<vUPost.x()<<" "<<vHUPre.y()<<" "<<vUPost.y() <<" "<< |
---|
503 | TMSINT<<" "<<vHAverNorm.x()<<" "<<vUPost.x()<<" "<<vHAverNorm.y()<<" "<<vUPost.y() <<" "<< |
---|
504 | vHAverNorm.x()<<" "<<vUPre.x()<<" "<<vHAverNorm.y()<<" "<<vUPre.y() <<" "<<sinpPre<<" "<<QAver<<G4endl; |
---|
505 | #endif |
---|
506 | |
---|
507 | transf[3][3] = (vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ); |
---|
508 | |
---|
509 | transf[3][4] = (vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ); |
---|
510 | |
---|
511 | // *** Zt |
---|
512 | transf[4][0] = pAver*(vVPost.x()*deltaPos.x()+vVPost.y()*deltaPos.y()+vVPost.z()*deltaPos.z()) |
---|
513 | *(1.+deltaPInv*pAver); |
---|
514 | |
---|
515 | transf[4][1] = ( sinThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) + |
---|
516 | OMcosThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) + |
---|
517 | TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())* |
---|
518 | (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver; |
---|
519 | #ifdef G4EVERBOSE |
---|
520 | if(iverbose >= 3)G4cout <<"ctransf41= " << transf[4][1] <<" "<< sinThetaAver<<" "<< OMcosThetaAver <<" "<<TMSINT<<" "<< vVPre <<" "<<vVPost <<" "<<vHVPre<<" "<<vHAverNorm <<" "<< QAver<<G4endl; |
---|
521 | #endif |
---|
522 | |
---|
523 | transf[4][2] = ( sinThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) + |
---|
524 | OMcosThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) + |
---|
525 | TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())* |
---|
526 | (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver; |
---|
527 | |
---|
528 | transf[4][3] = (vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ); |
---|
529 | |
---|
530 | transf[4][4] = (vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()); |
---|
531 | // if(iverbose >= 3) G4cout <<"ctransf44= " << transf[4][4] <<" "<< vVPre.x() <<" "<<vVPost.x() <<" "<< vVPre.y() <<" "<< vVPost.y() <<" "<< vVPre.z() <<" "<< vVPost.z() << G4endl; |
---|
532 | |
---|
533 | |
---|
534 | #ifdef G4EVERBOSE |
---|
535 | if( iverbose >= 1 ) G4cout << "G4EP: transf matrix computed " << transf << G4endl; |
---|
536 | #endif |
---|
537 | /* for( G4int ii=0;ii<5;ii++){ |
---|
538 | for( G4int jj=0;jj<5;jj++){ |
---|
539 | G4cout << transf[ii][jj] << " "; |
---|
540 | } |
---|
541 | G4cout << G4endl; |
---|
542 | } */ |
---|
543 | } |
---|
544 | } |
---|
545 | // end of calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION |
---|
546 | /* if( iverbose >= 1 ) G4cout << "G4EP: transf not updated but initialized " << theFirstStep << G4endl; |
---|
547 | if( theFirstStep ) { |
---|
548 | theTransfMat = transf; |
---|
549 | theFirstStep = false; |
---|
550 | }else{ |
---|
551 | theTransfMat = theTransfMat * transf; |
---|
552 | if( iverbose >= 1 ) G4cout << "G4EP: transf matrix accumulated" << theTransfMat << G4endl; |
---|
553 | } |
---|
554 | */ |
---|
555 | theTransfMat = transf; |
---|
556 | #ifdef G4EVERBOSE |
---|
557 | if( iverbose >= 1 ) G4cout << "G4EP: error matrix before transformation " << fError << G4endl; |
---|
558 | if( iverbose >= 2 ) G4cout << " tf * err " << theTransfMat * fError << G4endl |
---|
559 | << " transf matrix " << theTransfMat.T() << G4endl; |
---|
560 | #endif |
---|
561 | |
---|
562 | fError = fError.similarity(theTransfMat).T(); |
---|
563 | //- fError = transf * fError * transf.T(); |
---|
564 | #ifdef G4EVERBOSE |
---|
565 | if( iverbose >= 1 ) G4cout << "G4EP: error matrix propagated " << fError << G4endl; |
---|
566 | #endif |
---|
567 | |
---|
568 | //? S = B*S*BT S.similarity(B) |
---|
569 | //? R = S |
---|
570 | // not needed * *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL VARIABLES; |
---|
571 | |
---|
572 | PropagateErrorMSC( aTrack ); |
---|
573 | |
---|
574 | PropagateErrorIoni( aTrack ); |
---|
575 | |
---|
576 | return 0; |
---|
577 | } |
---|
578 | |
---|
579 | |
---|
580 | //------------------------------------------------------------------------ |
---|
581 | G4int G4ErrorFreeTrajState::PropagateErrorMSC( const G4Track* aTrack ) |
---|
582 | { |
---|
583 | G4ThreeVector vpPre = aTrack->GetMomentum()/GeV; |
---|
584 | G4double pPre = vpPre.mag(); |
---|
585 | G4double pBeta = pPre*pPre / (aTrack->GetTotalEnergy()/GeV); |
---|
586 | G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm; |
---|
587 | |
---|
588 | G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial(); |
---|
589 | G4double effZ, effA; |
---|
590 | CalculateEffectiveZandA( mate, effZ, effA ); |
---|
591 | |
---|
592 | #ifdef G4EVERBOSE |
---|
593 | if( iverbose >= 4 ) G4cout << "material " << mate->GetName() |
---|
594 | //<< " " << mate->GetZ() << " " << mate->GetA() |
---|
595 | << " " << effZ << " " << effA |
---|
596 | << " " << mate->GetDensity()/g*mole << " " << mate->GetRadlen()/cm << " " << mate->GetNuclearInterLength()/cm << G4endl; |
---|
597 | #endif |
---|
598 | |
---|
599 | G4double RI = stepLengthCm / (mate->GetRadlen()/cm); |
---|
600 | #ifdef G4EVERBOSE |
---|
601 | if( iverbose >= 4 ) G4cout << std::setprecision(6) << std::setw(6) << "G4EP:MSC: RI " << RI << " stepLengthCm " << stepLengthCm << " radlen " << (mate->GetRadlen()/cm) << " " << RI*1.e10 << G4endl; |
---|
602 | #endif |
---|
603 | G4double charge = aTrack->GetDynamicParticle()->GetCharge(); |
---|
604 | G4double DD = 1.8496E-4*RI*(charge/pBeta * charge/pBeta ); |
---|
605 | #ifdef G4EVERBOSE |
---|
606 | if( iverbose >= 3 ) G4cout << "G4EP:MSC: D*1E6= " << DD*1.E6 <<" pBeta " << pBeta << G4endl; |
---|
607 | #endif |
---|
608 | G4double S1 = DD*stepLengthCm*stepLengthCm/3.; |
---|
609 | G4double S2 = DD; |
---|
610 | G4double S3 = DD*stepLengthCm/2.; |
---|
611 | |
---|
612 | G4double CLA = std::sqrt( vpPre.x() * vpPre.x() + vpPre.y() * vpPre.y() )/pPre; |
---|
613 | #ifdef G4EVERBOSE |
---|
614 | if( iverbose >= 2 ) G4cout << std::setw(6) << "G4EP:MSC: RI " << RI << " S1 " << S1 << " S2 " << S2 << " S3 " << S3 << " CLA " << CLA << G4endl; |
---|
615 | #endif |
---|
616 | fError[1][1] += S2; |
---|
617 | fError[1][4] -= S3; |
---|
618 | fError[2][2] += S2/CLA/CLA; |
---|
619 | fError[2][3] += S3/CLA; |
---|
620 | fError[3][3] += S1; |
---|
621 | fError[4][4] += S1; |
---|
622 | |
---|
623 | #ifdef G4EVERBOSE |
---|
624 | if( iverbose >= 2 ) G4cout << "G4EP:MSC: error matrix propagated msc " << fError << G4endl; |
---|
625 | #endif |
---|
626 | |
---|
627 | return 0; |
---|
628 | } |
---|
629 | |
---|
630 | |
---|
631 | //------------------------------------------------------------------------ |
---|
632 | void G4ErrorFreeTrajState::CalculateEffectiveZandA( const G4Material* mate, G4double& effZ, G4double& effA ) |
---|
633 | { |
---|
634 | effZ = 0.; |
---|
635 | effA = 0.; |
---|
636 | G4int ii, nelem = mate->GetNumberOfElements(); |
---|
637 | const G4double* fracVec = mate->GetFractionVector(); |
---|
638 | for(ii=0; ii < nelem; ii++ ) { |
---|
639 | effZ += mate->GetElement( ii )->GetZ() * fracVec[ii]; |
---|
640 | effA += mate->GetElement( ii )->GetA() * fracVec[ii] /g*mole; |
---|
641 | } |
---|
642 | |
---|
643 | } |
---|
644 | |
---|
645 | |
---|
646 | //------------------------------------------------------------------------ |
---|
647 | G4int G4ErrorFreeTrajState::PropagateErrorIoni( const G4Track* aTrack ) |
---|
648 | { |
---|
649 | G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm; |
---|
650 | G4double DEDX2; |
---|
651 | if( stepLengthCm < 1.E-7 ) { |
---|
652 | DEDX2=0.; |
---|
653 | } |
---|
654 | // * Calculate xi factor (KeV). |
---|
655 | G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial(); |
---|
656 | G4double effZ, effA; |
---|
657 | CalculateEffectiveZandA( mate, effZ, effA ); |
---|
658 | |
---|
659 | G4double Etot = aTrack->GetTotalEnergy()/GeV; |
---|
660 | G4double beta = aTrack->GetMomentum().mag()/GeV / Etot; |
---|
661 | G4double mass = aTrack->GetDynamicParticle()->GetMass() / GeV; |
---|
662 | G4double gamma = Etot / mass; |
---|
663 | |
---|
664 | // * Calculate xi factor (KeV). |
---|
665 | G4double XI = 153.5*effZ*stepLengthCm*(mate->GetDensity()/mg*mole) / |
---|
666 | (effA*beta*beta); |
---|
667 | |
---|
668 | #ifdef G4EVERBOSE |
---|
669 | if( iverbose >= 2 ){ |
---|
670 | G4cout << "G4EP:IONI: XI " << XI << " beta " << beta << " gamma " << gamma << G4endl; |
---|
671 | G4cout << " density " << (mate->GetDensity()/mg*mole) << " effA " << effA << " step " << stepLengthCm << G4endl; |
---|
672 | } |
---|
673 | #endif |
---|
674 | // * Maximum energy transfer to atomic electron (KeV). |
---|
675 | G4double eta = beta*gamma; |
---|
676 | G4double etasq = eta*eta; |
---|
677 | G4double eMass = 0.51099906/GeV; |
---|
678 | G4double massRatio = eMass / mass; |
---|
679 | G4double F1 = 2*eMass*etasq; |
---|
680 | G4double F2 = 1. + 2. * massRatio * gamma + massRatio * massRatio; |
---|
681 | G4double Emax = 1.E+6*F1/F2; |
---|
682 | |
---|
683 | // * *** and now sigma**2 in GeV |
---|
684 | G4double dedxSq = XI*Emax*(1.-(beta*beta/2.))*1.E-12; |
---|
685 | #ifdef G4EVERBOSE |
---|
686 | if( iverbose >= 2 ) G4cout << "G4EP:IONI: DEDX2 " << dedxSq << " emass " << eMass << " Emax " << Emax << G4endl; |
---|
687 | #endif |
---|
688 | |
---|
689 | // if( iverbose >= 2 ) G4cout << "G4EP:IONI: Etot " << Etot << " DEDX2 " << dedxSq << " emass " << eMass << G4endl; |
---|
690 | |
---|
691 | G4double pPre6 = (aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV).mag(); |
---|
692 | pPre6 = std::pow(pPre6, 6 ); |
---|
693 | //Apply it to error |
---|
694 | fError[0][0] += Etot*Etot*dedxSq / pPre6; |
---|
695 | #ifdef G4EVERBOSE |
---|
696 | if( iverbose >= 2 ) G4cout << "G4:IONI getot " << Etot << " dedx2 " << dedxSq << " p " << pPre6 << G4endl; |
---|
697 | if( iverbose >= 2 ) G4cout << "G4EP:IONI: error_from_ionisation " << (Etot*Etot*dedxSq) / pPre6 << G4endl; |
---|
698 | #endif |
---|
699 | |
---|
700 | return 0; |
---|
701 | } |
---|
702 | |
---|