Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Friends | List of all members
G4ErrorFreeTrajState Class Reference

#include <G4ErrorFreeTrajState.hh>

+ Inheritance diagram for G4ErrorFreeTrajState:

Public Member Functions

 G4ErrorFreeTrajState ()
 
 G4ErrorFreeTrajState (const G4String &partName, const G4Point3D &pos, const G4Vector3D &mom, const G4ErrorTrajErr &errmat=G4ErrorTrajErr(5, 0))
 
 G4ErrorFreeTrajState (const G4ErrorSurfaceTrajState &tpOS)
 
 ~G4ErrorFreeTrajState ()
 
virtual G4int Update (const G4Track *aTrack)
 
virtual G4int PropagateError (const G4Track *aTrack)
 
virtual void Dump (std::ostream &out=G4cout) const
 
virtual void SetPosition (const G4Point3D pos)
 
virtual void SetMomentum (const G4Vector3D &mom)
 
void SetParameters (const G4Point3D &pos, const G4Vector3D &mom)
 
G4ErrorFreeTrajParam GetParameters () const
 
G4ErrorMatrix GetTransfMat () const
 
- Public Member Functions inherited from G4ErrorTrajState
 G4ErrorTrajState ()
 
 G4ErrorTrajState (const G4String &partType, const G4Point3D &pos, const G4Vector3D &mom, const G4ErrorTrajErr &errmat=G4ErrorTrajErr(5, 0))
 
 G4ErrorTrajState (const G4ErrorTrajState &)
 
 G4ErrorTrajState (G4ErrorTrajState &&)
 
virtual ~G4ErrorTrajState ()
 
G4ErrorTrajState & operator= (const G4ErrorTrajState &)
 
G4ErrorTrajState & operator= (G4ErrorTrajState &&)
 
void SetData (const G4String &partType, const G4Point3D &pos, const G4Vector3D &mom)
 
void BuildCharge ()
 
virtual G4int PropagateError (const G4Track *)
 
virtual G4int Update (const G4Track *)
 
void UpdatePosMom (const G4Point3D &pos, const G4Vector3D &mom)
 
void DumpPosMomError (std::ostream &out=G4cout) const
 
virtual void Dump (std::ostream &out=G4cout) const =0
 
const G4String & GetParticleType () const
 
void SetParticleType (const G4String &partType)
 
G4Point3D GetPosition () const
 
virtual void SetPosition (const G4Point3D pos)
 
G4Vector3D GetMomentum () const
 
virtual void SetMomentum (const G4Vector3D &mom)
 
G4ErrorTrajErr GetError () const
 
virtual void SetError (G4ErrorTrajErr em)
 
G4Track * GetG4Track () const
 
void SetG4Track (G4Track *trk)
 
G4double GetCharge () const
 
void SetCharge (G4double ch)
 
virtual G4eTSType GetTSType () const
 

Friends

std::ostream & operator<< (std::ostream &, const G4ErrorFreeTrajState &ts)
 

Additional Inherited Members

- Protected Attributes inherited from G4ErrorTrajState
G4String fParticleType
 
G4Point3D fPosition
 
G4Vector3D fMomentum
 
G4double fCharge = 0.
 
G4ErrorTrajErr fError
 
G4eTSType theTSType
 
G4Track * theG4Track = nullptr
 
G4int iverbose = 0
 

Detailed Description

Definition at line 64 of file G4ErrorFreeTrajState.hh.

Constructor & Destructor Documentation

◆ G4ErrorFreeTrajState() [1/3]

G4ErrorFreeTrajState::G4ErrorFreeTrajState ( )
inline

Definition at line 67 of file G4ErrorFreeTrajState.hh.

68 : theFirstStep(true)
69 {}

◆ G4ErrorFreeTrajState() [2/3]

G4ErrorFreeTrajState::G4ErrorFreeTrajState ( const G4String &  partName,
const G4Point3D &  pos,
const G4Vector3D &  mom,
const G4ErrorTrajErr &  errmat = G4ErrorTrajErr(5, 0) 
)

Definition at line 48 of file G4ErrorFreeTrajState.cc.

52 : G4ErrorTrajState(partName, pos, mom, errmat)
53{
54 fTrajParam = G4ErrorFreeTrajParam(pos, mom);
55 Init();
56}

◆ G4ErrorFreeTrajState() [3/3]

G4ErrorFreeTrajState::G4ErrorFreeTrajState ( const G4ErrorSurfaceTrajState &  tpOS)

Definition at line 59 of file G4ErrorFreeTrajState.cc.

60 : G4ErrorTrajState(tpSD.GetParticleType(), tpSD.GetPosition(),
61 tpSD.GetMomentum())
62{
63 // G4ThreeVector planeNormal = tpSD.GetPlaneNormal();
64 // G4double fPt = tpSD.GetMomentum()*planeNormal;//mom projected on normal to
65 // plane G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
66 // G4ThreeVector Psc = fPt * planeNormal +
67 // tpSDparam.GetPU()*tpSDparam.GetVectorU() + tpSD.GetPV()*tpSD.GetVectorW();
68
69 fTrajParam = G4ErrorFreeTrajParam(fPosition, fMomentum);
70 Init();
71
72 //----- Get the error matrix in SC coordinates
73 G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
74 G4double mom = fMomentum.mag();
75 G4double mom2 = fMomentum.mag2();
76 G4double TVW1 =
77 std::sqrt(mom2 / (mom2 + tpSDparam.GetPV() * tpSDparam.GetPV() +
78 tpSDparam.GetPW() * tpSDparam.GetPW()));
79 G4ThreeVector vTVW(TVW1, tpSDparam.GetPV() / mom * TVW1,
80 tpSDparam.GetPW() / mom * TVW1);
81 G4Vector3D vectorU = tpSDparam.GetVectorV().cross(tpSDparam.GetVectorW());
82 G4Vector3D vTN = vTVW.x() * vectorU + vTVW.y() * tpSDparam.GetVectorV() +
83 vTVW.z() * tpSDparam.GetVectorW();
84
85#ifdef G4EVERBOSE
86 if(iverbose >= 5)
87 {
88 G4double pc2 = std::asin(vTN.z());
89 G4double pc3 = std::atan(vTN.y() / vTN.x());
90
91 G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda()
92 << " diff " << pc2 - GetParameters().GetLambda() << G4endl;
93 G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi()
94 << " diff " << pc3 - GetParameters().GetPhi() << G4endl;
95 }
96#endif
97
98 //--- Get the unit vectors perp to P
99 G4double cosl = std::cos(GetParameters().GetLambda());
100 if(cosl < 1.E-30)
101 cosl = 1.E-30;
102 G4double cosl1 = 1. / cosl;
103 G4Vector3D vUN(-vTN.y() * cosl1, vTN.x() * cosl1, 0.);
104 G4Vector3D vVN(-vTN.z() * vUN.y(), vTN.z() * vUN.x(), cosl);
105
106 G4Vector3D vUperp = G4Vector3D(-fMomentum.y(), fMomentum.x(), 0.);
107 G4Vector3D vVperp = vUperp.cross(fMomentum);
108 vUperp *= 1. / vUperp.mag();
109 vVperp *= 1. / vVperp.mag();
110
111#ifdef G4EVERBOSE
112 if(iverbose >= 5)
113 {
114 G4cout << " CHECK: vUN " << vUN << " = " << vUperp << " diff "
115 << (vUN - vUperp).mag() << G4endl;
116 G4cout << " CHECK: vVN " << vVN << " = " << vVperp << " diff "
117 << (vVN - vVperp).mag() << G4endl;
118 }
119#endif
120
121 // get the dot products of vectors perpendicular to direction and vector
122 // defining SD plane
123 G4double dUU = vUperp * tpSD.GetVectorV();
124 G4double dUV = vUperp * tpSD.GetVectorW();
125 G4double dVU = vVperp * tpSD.GetVectorV();
126 G4double dVV = vVperp * tpSD.GetVectorW();
127
128 //--- Get transformation first
129 G4ErrorMatrix transfM(5, 5, 1);
130 //--- Get magnetic field
131 const G4Field* field = G4TransportationManager::GetTransportationManager()
132 ->GetFieldManager()
133 ->GetDetectorField();
134 G4ThreeVector dir = fTrajParam.GetDirection();
135 G4double invCosTheta = 1. / std::cos(dir.theta());
136 G4cout << " dir=" << dir << " invCosTheta " << invCosTheta << G4endl;
137
138 if(fCharge != 0 && field)
139 {
140 G4double pos1[3];
141 pos1[0] = fPosition.x() * cm;
142 pos1[1] = fPosition.y() * cm;
143 pos1[2] = fPosition.z() * cm;
144 G4double h1[3];
145 field->GetFieldValue(pos1, h1);
146 G4ThreeVector HPre = G4ThreeVector(h1[0], h1[1], h1[2]) / tesla * 10.;
147 G4double magHPre = HPre.mag();
148 G4double invP = 1. / fMomentum.mag();
149 G4double magHPreM = magHPre * invP;
150 if(magHPre != 0.)
151 {
152 G4double magHPreM2 = fCharge / magHPre;
153
154 G4double Q = -magHPreM * c_light;
155 G4double sinz = -HPre * vUperp * magHPreM2;
156 G4double cosz = HPre * vVperp * magHPreM2;
157
158 transfM[1][3] = -Q * dir.y() * sinz;
159 transfM[1][4] = -Q * dir.z() * sinz;
160 transfM[2][3] = -Q * dir.y() * cosz * invCosTheta;
161 transfM[2][4] = -Q * dir.z() * cosz * invCosTheta;
162 }
163 }
164
165 transfM[0][0] = 1.;
166 transfM[1][1] = dir.x() * dVU;
167 transfM[1][2] = dir.x() * dVV;
168 transfM[2][1] = dir.x() * dUU * invCosTheta;
169 transfM[2][2] = dir.x() * dUV * invCosTheta;
170 transfM[3][3] = dUU;
171 transfM[3][4] = dUV;
172 transfM[4][3] = dVU;
173 transfM[4][4] = dVV;
174
175 fError = G4ErrorTrajErr(tpSD.GetError().similarity(transfM));
176
177#ifdef G4EVERBOSE
178 if(iverbose >= 1)
179 G4cout << "error matrix SD2SC " << fError << G4endl;
180 if(iverbose >= 4)
181 G4cout << "G4ErrorFreeTrajState from SD " << *this << G4endl;
182#endif
183}
G4ErrorTrajErr
G4ErrorSymMatrix G4ErrorTrajErr
Definition: G4ErrorTrajErr.hh:43
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4double
double G4double
Definition: G4Types.hh:83
G4Vector3D
HepGeom::Vector3D< G4double > G4Vector3D
Definition: G4Vector3D.hh:34
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::theta
double theta() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::mag
double mag() const
G4ErrorFreeTrajParam::GetLambda
G4double GetLambda() const
Definition: G4ErrorFreeTrajParam.hh:84
G4ErrorFreeTrajParam::GetDirection
G4Vector3D GetDirection() const
Definition: G4ErrorFreeTrajParam.hh:81
G4ErrorFreeTrajState::GetParameters
G4ErrorFreeTrajParam GetParameters() const
Definition: G4ErrorFreeTrajState.hh:111
G4ErrorTrajState::fError
G4ErrorTrajErr fError
Definition: G4ErrorTrajState.hh:127
G4FieldManager::GetDetectorField
const G4Field * GetDetectorField() const
G4Field::GetFieldValue
virtual void GetFieldValue(const G4double Point[4], G4double *fieldArr) const =0
G4TransportationManager::GetTransportationManager
static G4TransportationManager * GetTransportationManager()
Definition: G4TransportationManager.cc:104
G4TransportationManager::GetFieldManager
G4FieldManager * GetFieldManager() const
HepGeom::BasicVector3D::cross
BasicVector3D< T > cross(const BasicVector3D< T > &v) const
Definition: BasicVector3D.h:281

◆ ~G4ErrorFreeTrajState()

G4ErrorFreeTrajState::~G4ErrorFreeTrajState ( )
inline

Definition at line 78 of file G4ErrorFreeTrajState.hh.

78{}

Member Function Documentation

◆ Dump()

void G4ErrorFreeTrajState::Dump ( std::ostream &  out = G4cout) const
virtual

Implements G4ErrorTrajState.

Definition at line 195 of file G4ErrorFreeTrajState.cc.

195{ out << *this; }

◆ GetParameters()

G4ErrorFreeTrajParam G4ErrorFreeTrajState::GetParameters ( ) const
inline

Definition at line 111 of file G4ErrorFreeTrajState.hh.

111{ return fTrajParam; }

Referenced by G4ErrorSurfaceTrajState::BuildErrorMatrix(), and G4ErrorFreeTrajState().

◆ GetTransfMat()

G4ErrorMatrix G4ErrorFreeTrajState::GetTransfMat ( ) const
inline

Definition at line 113 of file G4ErrorFreeTrajState.hh.

113{ return theTransfMat; }

◆ PropagateError()

G4int G4ErrorFreeTrajState::PropagateError ( const G4Track *  aTrack)
virtual

Reimplemented from G4ErrorTrajState.

Definition at line 223 of file G4ErrorFreeTrajState.cc.

224{
225 G4double stepLengthCm = aTrack->GetStep()->GetStepLength() / cm;
226 if(G4ErrorPropagatorData::GetErrorPropagatorData()->GetStage() ==
227 G4ErrorStage_Deflation)
228 stepLengthCm *= -1.;
229
230 G4double kCarTolerance =
231 G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
232
233 if(std::fabs(stepLengthCm) <= kCarTolerance / cm)
234 return 0;
235
236#ifdef G4EVERBOSE
237 if(iverbose >= 2)
238 G4cout << " G4ErrorFreeTrajState::PropagateError " << G4endl;
239 G4cout << "G4EP: iverbose=" << iverbose << G4endl;
240#endif
241
242 // * *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
243 G4Point3D vposPost = aTrack->GetPosition() / cm;
244 G4Vector3D vpPost = aTrack->GetMomentum() / GeV;
245 // G4Point3D vposPre = fPosition/cm;
246 // G4Vector3D vpPre = fMomentum/GeV;
247 G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition() / cm;
248 G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum() / GeV;
249 // correct to avoid propagation along Z
250 if(vpPre.mag() == vpPre.z())
251 vpPre.setX(1.E-6 * MeV);
252 if(vpPost.mag() == vpPost.z())
253 vpPost.setX(1.E-6 * MeV);
254
255 G4double pPre = vpPre.mag();
256 G4double pPost = vpPost.mag();
257#ifdef G4EVERBOSE
258 if(iverbose >= 2)
259 {
260 G4cout << "G4EP: vposPre " << vposPre << G4endl << "G4EP: vposPost "
261 << vposPost << G4endl;
262 G4cout << "G4EP: vpPre " << vpPre << G4endl << "G4EP: vpPost " << vpPost
263 << G4endl;
264 G4cout << " err start step " << fError << G4endl;
265 G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl;
266 }
267#endif
268
269 if(pPre == 0. || pPost == 0)
270 return 2;
271 G4double pInvPre = 1. / pPre;
272 G4double pInvPost = 1. / pPost;
273 G4double deltaPInv = pInvPost - pInvPre;
274 if(iverbose >= 2)
275 G4cout << "G4EP: pInvPre" << pInvPre << " pInvPost:" << pInvPost
276 << " deltaPInv:" << deltaPInv << G4endl;
277
278 G4Vector3D vpPreNorm = vpPre * pInvPre;
279 G4Vector3D vpPostNorm = vpPost * pInvPost;
280 if(iverbose >= 2)
281 G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm
282 << G4endl;
283 // return if propagation along Z??
284 if(1. - std::fabs(vpPreNorm.z()) < kCarTolerance)
285 return 4;
286 if(1. - std::fabs(vpPostNorm.z()) < kCarTolerance)
287 return 4;
288 G4double sinpPre =
289 std::sin(vpPreNorm.theta()); // cosine perpendicular to pPre = sine pPre
290 G4double sinpPost =
291 std::sin(vpPostNorm.theta()); // cosine perpendicular to pPost = sine pPost
292 G4double sinpPostInv = 1. / std::sin(vpPostNorm.theta());
293
294#ifdef G4EVERBOSE
295 if(iverbose >= 2)
296 G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl;
297#endif
298 //* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR
299 //* *** NEUTRAL PARTICLE OR FIELDFREE REGION
300 G4ErrorMatrix transf(5, 5, 0);
301
302 transf[3][2] = stepLengthCm * sinpPost;
303 transf[4][1] = stepLengthCm;
304 for(auto ii = 0; ii < 5; ++ii)
305 {
306 transf[ii][ii] = 1.;
307 }
308#ifdef G4EVERBOSE
309 if(iverbose >= 2)
310 {
311 G4cout << "G4EP: transf matrix neutral " << transf;
312 }
313#endif
314
315 // charge X propagation direction
316 G4double charge = aTrack->GetDynamicParticle()->GetCharge();
317 if(G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() ==
318 G4ErrorMode_PropBackwards)
319 {
320 charge *= -1.;
321 }
322 // G4cout << " charge " << charge << G4endl;
323 // t check if particle has charge
324 // t if( charge == 0 ) goto 45;
325 // check if the magnetic field is = 0.
326
327 // position is from geant4, it is assumed to be in mm (for debugging,
328 // eventually it will not be transformed) it is assumed vposPre[] is in cm and
329 // pos1[] is in mm.
330 G4double pos1[3];
331 pos1[0] = vposPre.x() * cm;
332 pos1[1] = vposPre.y() * cm;
333 pos1[2] = vposPre.z() * cm;
334 G4double pos2[3];
335 pos2[0] = vposPost.x() * cm;
336 pos2[1] = vposPost.y() * cm;
337 pos2[2] = vposPost.z() * cm;
338 G4double h1[3], h2[3];
339
340 const G4Field* field = G4TransportationManager::GetTransportationManager()
341 ->GetFieldManager()
342 ->GetDetectorField();
343 if(!field)
344 return 0; // goto 45
345
346 // calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION
347 if(charge != 0. && field)
348 {
349 field->GetFieldValue(pos1,
350 h1); // here pos1[], pos2[] are in mm, not changed
351 field->GetFieldValue(pos2, h2);
352 G4ThreeVector HPre =
353 G4ThreeVector(h1[0], h1[1], h1[2]) / tesla *
354 10.; // 10. is to get same dimensions as GEANT3 (kilogauss)
355 G4ThreeVector HPost = G4ThreeVector(h2[0], h2[1], h2[2]) / tesla * 10.;
356 G4double magHPre = HPre.mag();
357 G4double magHPost = HPost.mag();
358#ifdef G4EVERBOSE
359 if(iverbose >= 2)
360 {
361 G4cout << "G4EP: h1 = " << h1[0] << ", " << h1[1] << ", " << h1[2]
362 << G4endl;
363 G4cout << "G4EP: pos1/mm = " << pos1[0] << ", " << pos1[1] << ", "
364 << pos1[2] << G4endl;
365 G4cout << "G4EP: pos2/mm = " << pos2[0] << ", " << pos2[1] << ", "
366 << pos2[2] << G4endl;
367 G4cout << "G4EP: B-filed in KGauss HPre " << HPre << G4endl
368 << "G4EP: in KGauss HPost " << HPost << G4endl;
369 }
370#endif
371
372 if(magHPre + magHPost != 0.)
373 {
374 //* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2
375 G4double gam;
376 if(magHPost != 0.)
377 {
378 gam = HPost * vpPostNorm / magHPost;
379 }
380 else
381 {
382 gam = HPre * vpPreNorm / magHPre;
383 }
384
385 // G4eMagneticLimitsProcess will limit the step, but based on an straight
386 // line trajectory
387 G4double alphaSqr = 1. - gam * gam;
388 G4double diffHSqr = (HPre * pInvPre - HPost * pInvPost).mag2();
389 G4double delhp6Sqr = 300. * 300.;
390#ifdef G4EVERBOSE
391 if(iverbose >= 2)
392 {
393 G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr
394 << " diffHSqr " << diffHSqr << G4endl;
395 G4cout << " alpha= " << std::sqrt(alphaSqr) << G4endl;
396 }
397#endif
398 if(diffHSqr * alphaSqr > delhp6Sqr)
399 return 3;
400
401 //* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT
402 G4double pInvAver = 1. / (pInvPre + pInvPost);
403 G4double CFACT8 = 2.997925E-4;
404 // G4double HAver
405 G4ThreeVector vHAverNorm((HPre * pInvPre + HPost * pInvPost) * pInvAver *
406 charge * CFACT8);
407 G4double HAver = vHAverNorm.mag();
408 G4double invHAver = 1. / HAver;
409 vHAverNorm *= invHAver;
410#ifdef G4EVERBOSE
411 if(iverbose >= 2)
412 G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver
413 << " charge " << charge << G4endl;
414#endif
415
416 G4double pAver = (pPre + pPost) * 0.5;
417 G4double QAver = -HAver / pAver;
418 G4double thetaAver = QAver * stepLengthCm;
419 G4double sinThetaAver = std::sin(thetaAver);
420 G4double cosThetaAver = std::cos(thetaAver);
421 G4double gamma = vHAverNorm * vpPostNorm;
422 G4ThreeVector AN2 = vHAverNorm.cross(vpPostNorm);
423
424#ifdef G4EVERBOSE
425 if(iverbose >= 2)
426 G4cout << " G4EP: AN2 " << AN2 << " gamma:" << gamma
427 << " theta=" << thetaAver << G4endl;
428#endif
429 G4double AU = 1. / vpPreNorm.perp();
430 // t G4ThreeVector vU( vpPreNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU );
431 G4ThreeVector vUPre(-AU * vpPreNorm.y(), AU * vpPreNorm.x(), 0.);
432 G4ThreeVector vVPre(-vpPreNorm.z() * vUPre.y(), vpPreNorm.z() * vUPre.x(),
433 vpPreNorm.x() * vUPre.y() -
434 vpPreNorm.y() * vUPre.x());
435
436 //
437 AU = 1. / vpPostNorm.perp();
438 // t G4ThreeVector vU( vpPostNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU
439 // );
440 G4ThreeVector vUPost(-AU * vpPostNorm.y(), AU * vpPostNorm.x(), 0.);
441 G4ThreeVector vVPost(
442 -vpPostNorm.z() * vUPost.y(), vpPostNorm.z() * vUPost.x(),
443 vpPostNorm.x() * vUPost.y() - vpPostNorm.y() * vUPost.x());
444#ifdef G4EVERBOSE
445 G4cout << " vpPostNorm " << vpPostNorm << G4endl;
446 if(iverbose >= 2)
447 G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre
448 << " vUPost " << vUPost << " vVPost " << vVPost << G4endl;
449#endif
450 G4Point3D deltaPos(vposPre - vposPost);
451
452 // * *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2
453 // * *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT
454 // * *** TAKEN INTO ACCOUNT
455
456 G4double QP = QAver * pAver; // = -HAver
457#ifdef G4EVERBOSE
458 if(iverbose >= 2)
459 G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver
460 << G4endl;
461#endif
462 G4double ANV =
463 -(vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y());
464 G4double ANU =
465 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
466 vHAverNorm.z() * vVPost.z());
467 G4double OMcosThetaAver = 1. - cosThetaAver;
468#ifdef G4EVERBOSE
469 if(iverbose >= 2)
470 G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver "
471 << cosThetaAver << " thetaAver " << thetaAver << " QAver "
472 << QAver << " stepLengthCm " << stepLengthCm << G4endl;
473#endif
474 G4double TMSINT = thetaAver - sinThetaAver;
475#ifdef G4EVERBOSE
476 if(iverbose >= 2)
477 G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl;
478#endif
479
480 G4ThreeVector vHUPre(
481 -vHAverNorm.z() * vUPre.y(), vHAverNorm.z() * vUPre.x(),
482 vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x());
483#ifdef G4EVERBOSE
484 // if( iverbose >= 2 ) G4cout << "G4EP: HUPre(1) " << vHUPre.x() << " "
485 // << vHAverNorm.z() << " " << vUPre.y() << G4endl;
486#endif
487 G4ThreeVector vHVPre(
488 vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(),
489 vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(),
490 vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x());
491#ifdef G4EVERBOSE
492 if(iverbose >= 2)
493 G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl;
494#endif
495
496 //------------------- COMPUTE MATRIX
497 //---------- 1/P
498
499 transf[0][0] =
500 1. -
501 deltaPInv * pAver *
502 (1. + (vpPostNorm.x() * deltaPos.x() + vpPostNorm.y() * deltaPos.y() +
503 vpPostNorm.z() * deltaPos.z()) /
504 stepLengthCm) +
505 2. * deltaPInv * pAver;
506
507 transf[0][1] =
508 -deltaPInv / thetaAver *
509 (TMSINT * gamma *
510 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
511 vHAverNorm.z() * vVPre.z()) +
512 sinThetaAver *
513 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
514 vVPre.z() * vpPostNorm.z()) +
515 OMcosThetaAver *
516 (vHVPre.x() * vpPostNorm.x() + vHVPre.y() * vpPostNorm.y() +
517 vHVPre.z() * vpPostNorm.z()));
518
519 transf[0][2] =
520 -sinpPre * deltaPInv / thetaAver *
521 (TMSINT * gamma *
522 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()) +
523 sinThetaAver *
524 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) +
525 OMcosThetaAver *
526 (vHUPre.x() * vpPostNorm.x() + vHUPre.y() * vpPostNorm.y() +
527 vHUPre.z() * vpPostNorm.z()));
528
529 transf[0][3] = -deltaPInv / stepLengthCm *
530 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y());
531
532 transf[0][4] = -deltaPInv / stepLengthCm *
533 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
534 vVPre.z() * vpPostNorm.z());
535
536 // *** Lambda
537 transf[1][0] =
538 -QP * ANV *
539 (vpPostNorm.x() * deltaPos.x() + vpPostNorm.y() * deltaPos.y() +
540 vpPostNorm.z() * deltaPos.z()) *
541 (1. + deltaPInv * pAver);
542#ifdef G4EVERBOSE
543 if(iverbose >= 3)
544 G4cout << "ctransf10= " << transf[1][0] << " " << -QP << " " << ANV
545 << " " << vpPostNorm.x() << " " << deltaPos.x() << " "
546 << vpPostNorm.y() << " " << deltaPos.y() << " " << vpPostNorm.z()
547 << " " << deltaPos.z() << " " << deltaPInv << " " << pAver
548 << G4endl;
549#endif
550
551 transf[1][1] =
552 cosThetaAver * (vVPre.x() * vVPost.x() + vVPre.y() * vVPost.y() +
553 vVPre.z() * vVPost.z()) +
554 sinThetaAver * (vHVPre.x() * vVPost.x() + vHVPre.y() * vVPost.y() +
555 vHVPre.z() * vVPost.z()) +
556 OMcosThetaAver *
557 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
558 vHAverNorm.z() * vVPre.z()) *
559 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
560 vHAverNorm.z() * vVPost.z()) +
561 ANV * (-sinThetaAver *
562 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
563 vVPre.z() * vpPostNorm.z()) +
564 OMcosThetaAver * (vVPre.x() * AN2.x() + vVPre.y() * AN2.y() +
565 vVPre.z() * AN2.z()) -
566 TMSINT * gamma *
567 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
568 vHAverNorm.z() * vVPre.z()));
569
570 transf[1][2] =
571 cosThetaAver * (vUPre.x() * vVPost.x() + vUPre.y() * vVPost.y()) +
572 sinThetaAver * (vHUPre.x() * vVPost.x() + vHUPre.y() * vVPost.y() +
573 vHUPre.z() * vVPost.z()) +
574 OMcosThetaAver *
575 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()) *
576 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
577 vHAverNorm.z() * vVPost.z()) +
578 ANV * (-sinThetaAver *
579 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) +
580 OMcosThetaAver * (vUPre.x() * AN2.x() + vUPre.y() * AN2.y()) -
581 TMSINT * gamma *
582 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()));
583 transf[1][2] = sinpPre * transf[1][2];
584
585 transf[1][3] = -QAver * ANV *
586 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y());
587
588 transf[1][4] = -QAver * ANV *
589 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
590 vVPre.z() * vpPostNorm.z());
591
592 // *** Phi
593
594 transf[2][0] =
595 -QP * ANU *
596 (vpPostNorm.x() * deltaPos.x() + vpPostNorm.y() * deltaPos.y() +
597 vpPostNorm.z() * deltaPos.z()) *
598 sinpPostInv * (1. + deltaPInv * pAver);
599#ifdef G4EVERBOSE
600 if(iverbose >= 3)
601 G4cout << "ctransf20= " << transf[2][0] << " " << -QP << " " << ANU
602 << " " << vpPostNorm.x() << " " << deltaPos.x() << " "
603 << vpPostNorm.y() << " " << deltaPos.y() << " " << vpPostNorm.z()
604 << " " << deltaPos.z() << " " << sinpPostInv << " " << deltaPInv
605 << " " << pAver << G4endl;
606#endif
607 transf[2][1] =
608 cosThetaAver * (vVPre.x() * vUPost.x() + vVPre.y() * vUPost.y()) +
609 sinThetaAver * (vHVPre.x() * vUPost.x() + vHVPre.y() * vUPost.y()) +
610 OMcosThetaAver *
611 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
612 vHAverNorm.z() * vVPre.z()) *
613 (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) +
614 ANU * (-sinThetaAver *
615 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
616 vVPre.z() * vpPostNorm.z()) +
617 OMcosThetaAver * (vVPre.x() * AN2.x() + vVPre.y() * AN2.y() +
618 vVPre.z() * AN2.z()) -
619 TMSINT * gamma *
620 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
621 vHAverNorm.z() * vVPre.z()));
622 transf[2][1] = sinpPostInv * transf[2][1];
623
624 transf[2][2] =
625 cosThetaAver * (vUPre.x() * vUPost.x() + vUPre.y() * vUPost.y()) +
626 sinThetaAver * (vHUPre.x() * vUPost.x() + vHUPre.y() * vUPost.y()) +
627 OMcosThetaAver *
628 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()) *
629 (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) +
630 ANU * (-sinThetaAver *
631 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) +
632 OMcosThetaAver * (vUPre.x() * AN2.x() + vUPre.y() * AN2.y()) -
633 TMSINT * gamma *
634 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()));
635 transf[2][2] = sinpPostInv * sinpPre * transf[2][2];
636
637 transf[2][3] = -QAver * ANU *
638 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) *
639 sinpPostInv;
640#ifdef G4EVERBOSE
641 if(iverbose >= 3)
642 G4cout << "ctransf23= " << transf[2][3] << " " << -QAver << " " << ANU
643 << " " << vUPre.x() << " " << vpPostNorm.x() << " " << vUPre.y()
644 << " " << vpPostNorm.y() << " " << sinpPostInv << G4endl;
645#endif
646
647 transf[2][4] = -QAver * ANU *
648 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
649 vVPre.z() * vpPostNorm.z()) *
650 sinpPostInv;
651
652 // *** Yt
653
654 transf[3][0] = pAver *
655 (vUPost.x() * deltaPos.x() + vUPost.y() * deltaPos.y()) *
656 (1. + deltaPInv * pAver);
657#ifdef G4EVERBOSE
658 if(iverbose >= 3)
659 G4cout << "ctransf30= " << transf[3][0] << " " << pAver << " "
660 << vUPost.x() << " " << deltaPos.x() << " " << vUPost.y() << " "
661 << deltaPos.y() << " " << deltaPInv << " " << pAver << G4endl;
662#endif
663
664 transf[3][1] =
665 (sinThetaAver * (vVPre.x() * vUPost.x() + vVPre.y() * vUPost.y()) +
666 OMcosThetaAver * (vHVPre.x() * vUPost.x() + vHVPre.y() * vUPost.y()) +
667 TMSINT * (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) *
668 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
669 vHAverNorm.z() * vVPre.z())) /
670 QAver;
671
672 transf[3][2] =
673 (sinThetaAver * (vUPre.x() * vUPost.x() + vUPre.y() * vUPost.y()) +
674 OMcosThetaAver * (vHUPre.x() * vUPost.x() + vHUPre.y() * vUPost.y()) +
675 TMSINT * (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) *
676 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y())) *
677 sinpPre / QAver;
678#ifdef G4EVERBOSE
679 if(iverbose >= 3)
680 G4cout << "ctransf32= " << transf[3][2] << " " << sinThetaAver << " "
681 << vUPre.x() << " " << vUPost.x() << " " << vUPre.y() << " "
682 << vUPost.y() << " " << OMcosThetaAver << " " << vHUPre.x()
683 << " " << vUPost.x() << " " << vHUPre.y() << " " << vUPost.y()
684 << " " << TMSINT << " " << vHAverNorm.x() << " " << vUPost.x()
685 << " " << vHAverNorm.y() << " " << vUPost.y() << " "
686 << vHAverNorm.x() << " " << vUPre.x() << " " << vHAverNorm.y()
687 << " " << vUPre.y() << " " << sinpPre << " " << QAver << G4endl;
688#endif
689
690 transf[3][3] = (vUPre.x() * vUPost.x() + vUPre.y() * vUPost.y());
691
692 transf[3][4] = (vVPre.x() * vUPost.x() + vVPre.y() * vUPost.y());
693
694 // *** Zt
695 transf[4][0] = pAver *
696 (vVPost.x() * deltaPos.x() + vVPost.y() * deltaPos.y() +
697 vVPost.z() * deltaPos.z()) *
698 (1. + deltaPInv * pAver);
699
700 transf[4][1] =
701 (sinThetaAver * (vVPre.x() * vVPost.x() + vVPre.y() * vVPost.y() +
702 vVPre.z() * vVPost.z()) +
703 OMcosThetaAver * (vHVPre.x() * vVPost.x() + vHVPre.y() * vVPost.y() +
704 vHVPre.z() * vVPost.z()) +
705 TMSINT *
706 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
707 vHAverNorm.z() * vVPost.z()) *
708 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
709 vHAverNorm.z() * vVPre.z())) /
710 QAver;
711#ifdef G4EVERBOSE
712 if(iverbose >= 3)
713 G4cout << "ctransf41= " << transf[4][1] << " " << sinThetaAver << " "
714 << OMcosThetaAver << " " << TMSINT << " " << vVPre << " "
715 << vVPost << " " << vHVPre << " " << vHAverNorm << " " << QAver
716 << G4endl;
717#endif
718
719 transf[4][2] =
720 (sinThetaAver * (vUPre.x() * vVPost.x() + vUPre.y() * vVPost.y()) +
721 OMcosThetaAver * (vHUPre.x() * vVPost.x() + vHUPre.y() * vVPost.y() +
722 vHUPre.z() * vVPost.z()) +
723 TMSINT *
724 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
725 vHAverNorm.z() * vVPost.z()) *
726 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y())) *
727 sinpPre / QAver;
728
729 transf[4][3] = (vUPre.x() * vVPost.x() + vUPre.y() * vVPost.y());
730
731 transf[4][4] = (vVPre.x() * vVPost.x() + vVPre.y() * vVPost.y() +
732 vVPre.z() * vVPost.z());
733 // if(iverbose >= 3) G4cout <<"ctransf44= " << transf[4][4] <<" "<<
734 // vVPre.x() <<" "<<vVPost.x() <<" "<< vVPre.y() <<" "<< vVPost.y() <<"
735 // "<< vVPre.z() <<" "<< vVPost.z() << G4endl;
736
737#ifdef G4EVERBOSE
738 if(iverbose >= 1)
739 G4cout << "G4EP: transf matrix computed " << transf << G4endl;
740#endif
741 /* for( G4int ii=0;ii<5;ii++){
742 for( G4int jj=0;jj<5;jj++){
743 G4cout << transf[ii][jj] << " ";
744 }
745 G4cout << G4endl;
746 } */
747 }
748 }
749 // end of calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE
750 // REGION
751 /* if( iverbose >= 1 ) G4cout << "G4EP: transf not updated but initialized "
752 << theFirstStep << G4endl; if( theFirstStep ) { theTransfMat = transf;
753 theFirstStep = false;
754 }else{
755 theTransfMat = theTransfMat * transf;
756 if( iverbose >= 1 ) G4cout << "G4EP: transf matrix accumulated" <<
757 theTransfMat << G4endl;
758 }
759 */
760 theTransfMat = transf;
761#ifdef G4EVERBOSE
762 if(iverbose >= 1)
763 G4cout << "G4EP: error matrix before transformation " << fError << G4endl;
764 if(iverbose >= 2)
765 G4cout << " tf * err " << theTransfMat * fError << G4endl
766 << " transf matrix " << theTransfMat.T() << G4endl;
767#endif
768
769 fError = fError.similarity(theTransfMat).T();
770 //- fError = transf * fError * transf.T();
771#ifdef G4EVERBOSE
772 if(iverbose >= 1)
773 G4cout << "G4EP: error matrix propagated " << fError << G4endl;
774#endif
775
776 //? S = B*S*BT S.similarity(B)
777 //? R = S
778 // not needed * *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL
779 // VARIABLES;
780
781 PropagateErrorMSC(aTrack);
782
783 PropagateErrorIoni(aTrack);
784
785 return 0;
786}
kCarTolerance
const G4double kCarTolerance
Definition: G4BoundingEnvelope.cc:38
G4ErrorMode_PropBackwards
@ G4ErrorMode_PropBackwards
Definition: G4ErrorPropagatorData.hh:44
G4ErrorStage_Deflation
@ G4ErrorStage_Deflation
Definition: G4ErrorPropagatorData.hh:60
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
G4DynamicParticle::GetCharge
G4double GetCharge() const
G4ErrorMatrix::T
G4ErrorMatrix T() const
Definition: G4ErrorMatrix.cc:399
G4ErrorPropagatorData::GetErrorPropagatorData
static G4ErrorPropagatorData * GetErrorPropagatorData()
Definition: G4ErrorPropagatorData.cc:46
G4ErrorSymMatrix::similarity
G4ErrorSymMatrix similarity(const G4ErrorMatrix &m1) const
Definition: G4ErrorSymMatrix.cc:629
G4ErrorSymMatrix::T
G4ErrorSymMatrix T() const
G4GeometryTolerance::GetSurfaceTolerance
G4double GetSurfaceTolerance() const
Definition: G4GeometryTolerance.cc:72
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:58
G4StepPoint::GetMomentum
G4ThreeVector GetMomentum() const
G4StepPoint::GetPosition
const G4ThreeVector & GetPosition() const
G4Step::GetPreStepPoint
G4StepPoint * GetPreStepPoint() const
G4Step::GetStepLength
G4double GetStepLength() const
G4Track::GetPosition
const G4ThreeVector & GetPosition() const
G4Track::GetMomentum
G4ThreeVector GetMomentum() const
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4Track::GetStep
const G4Step * GetStep() const

Referenced by G4ErrorPropagator::MakeOneStep().

◆ SetMomentum()

virtual void G4ErrorFreeTrajState::SetMomentum ( const G4Vector3D &  mom)
inlinevirtual

Reimplemented from G4ErrorTrajState.

Definition at line 99 of file G4ErrorFreeTrajState.hh.

100 {
101 SetParameters(fPosition, mom);
102 }
G4ErrorFreeTrajState::SetParameters
void SetParameters(const G4Point3D &pos, const G4Vector3D &mom)
Definition: G4ErrorFreeTrajState.hh:104

◆ SetParameters()

void G4ErrorFreeTrajState::SetParameters ( const G4Point3D &  pos,
const G4Vector3D &  mom 
)
inline

Definition at line 104 of file G4ErrorFreeTrajState.hh.

105 {
106 fPosition = pos;
107 fMomentum = mom;
108 fTrajParam.SetParameters(pos, mom);
109 }
G4ErrorFreeTrajParam::SetParameters
void SetParameters(const G4Point3D &pos, const G4Vector3D &mom)
Definition: G4ErrorFreeTrajParam.cc:46

Referenced by SetMomentum(), and SetPosition().

◆ SetPosition()

virtual void G4ErrorFreeTrajState::SetPosition ( const G4Point3D  pos)
inlinevirtual

Reimplemented from G4ErrorTrajState.

Definition at line 94 of file G4ErrorFreeTrajState.hh.

95 {
96 SetParameters(pos, fMomentum);
97 }

◆ Update()

G4int G4ErrorFreeTrajState::Update ( const G4Track *  aTrack)
virtual

Reimplemented from G4ErrorTrajState.

Definition at line 198 of file G4ErrorFreeTrajState.cc.

199{
200 G4int ierr = 0;
201 fTrajParam.Update(aTrack);
202 UpdatePosMom(aTrack->GetPosition(), aTrack->GetMomentum());
203 return ierr;
204}
G4int
int G4int
Definition: G4Types.hh:85
G4ErrorFreeTrajParam::Update
void Update(const G4Track *aTrack)
Definition: G4ErrorFreeTrajParam.cc:81
G4ErrorTrajState::UpdatePosMom
void UpdatePosMom(const G4Point3D &pos, const G4Vector3D &mom)
Definition: G4ErrorTrajState.cc:121

Referenced by G4ErrorPropagator::MakeOneStep().

Friends And Related Function Documentation

◆ operator<<

std::ostream & operator<< ( std::ostream &  out,
const G4ErrorFreeTrajState &  ts 
)
friend

Definition at line 207 of file G4ErrorFreeTrajState.cc.

208{
209 std::ios::fmtflags orig_flags = out.flags();
210
211 out.setf(std::ios::fixed, std::ios::floatfield);
212
213 ts.DumpPosMomError(out);
214
215 out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl;
216
217 out.flags(orig_flags);
218
219 return out;
220}
G4ErrorTrajState::DumpPosMomError
void DumpPosMomError(std::ostream &out=G4cout) const
Definition: G4ErrorTrajState.cc:156
The documentation for this class was generated from the following files: