Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Friends | List of all members
G4ErrorFreeTrajState Class Reference

#include <G4ErrorFreeTrajState.hh>

Inheritance diagram for G4ErrorFreeTrajState:
Inheritance graph
[legend]
Collaboration diagram for G4ErrorFreeTrajState:
Collaboration graph
[legend]

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))
 
virtual ~G4ErrorTrajState ()
 
void SetData (const G4String &partType, const G4Point3D &pos, const G4Vector3D &mom)
 
void BuildCharge ()
 
void UpdatePosMom (const G4Point3D &pos, const G4Vector3D &mom)
 
void DumpPosMomError (std::ostream &out=G4cout) const
 
const G4String & GetParticleType () const
 
void SetParticleType (const G4String &partType)
 
G4Point3D GetPosition () const
 
G4Vector3D GetMomentum () const
 
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
 
G4ErrorTrajErr fError
 
G4eTSType theTSType
 
G4Track * theG4Track
 
G4int iverbose
 

Detailed Description

Definition at line 65 of file G4ErrorFreeTrajState.hh.

Constructor & Destructor Documentation

G4ErrorFreeTrajState::G4ErrorFreeTrajState ( )
inline

Definition at line 69 of file G4ErrorFreeTrajState.hh.

69 : theFirstStep(true) {}
G4ErrorFreeTrajState::G4ErrorFreeTrajState ( const G4String &  partName,
const G4Point3D &  pos,
const G4Vector3D &  mom,
const G4ErrorTrajErr &  errmat = G4ErrorTrajErr(5,0) 
)

Definition at line 49 of file G4ErrorFreeTrajState.cc.

49  : G4ErrorTrajState( partName, pos, mom, errmat )
50 {
51  fTrajParam = G4ErrorFreeTrajParam( pos, mom );
52  Init();
53 }
G4ErrorFreeTrajState::G4ErrorFreeTrajState ( const G4ErrorSurfaceTrajState &  tpOS)

Definition at line 57 of file G4ErrorFreeTrajState.cc.

57  : G4ErrorTrajState( tpSD.GetParticleType(), tpSD.GetPosition(), tpSD.GetMomentum() )
58 {
59  // G4ThreeVector planeNormal = tpSD.GetPlaneNormal();
60  // G4double fPt = tpSD.GetMomentum()*planeNormal;//mom projected on normal to plane
61  // G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
62  // G4ThreeVector Psc = fPt * planeNormal + tpSDparam.GetPU()*tpSDparam.GetVectorU() + tpSD.GetPV()*tpSD.GetVectorW();
63 
64  fTrajParam = G4ErrorFreeTrajParam( fPosition, fMomentum );
65  Init();
66 
67  //----- Get the error matrix in SC coordinates
68  G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
69  G4double mom = fMomentum.mag();
70  G4double mom2 = fMomentum.mag2();
71  G4double TVW1 = std::sqrt( mom2 / ( mom2 + tpSDparam.GetPV()*tpSDparam.GetPV() + tpSDparam.GetPW()*tpSDparam.GetPW()) );
72  G4ThreeVector vTVW( TVW1, tpSDparam.GetPV()/mom * TVW1, tpSDparam.GetPW()/mom * TVW1 );
73  G4Vector3D vectorU = tpSDparam.GetVectorV().cross( tpSDparam.GetVectorW() );
74  G4Vector3D vTN = vTVW.x()*vectorU + vTVW.y()*tpSDparam.GetVectorV() + vTVW.z()*tpSDparam.GetVectorW();
75 
76 #ifdef G4EVERBOSE
77  if( iverbose >= 5){
78  G4double pc2 = std::asin( vTN.z() );
79  G4double pc3 = std::atan (vTN.y()/vTN.x());
80 
81  G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda() << " diff " << pc2-GetParameters().GetLambda() << G4endl;
82  G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi() << " diff " << pc3-GetParameters().GetPhi() << G4endl;
83  }
84 #endif
85 
86  //--- Get the unit vectors perp to P
87  G4double cosl = std::cos( GetParameters().GetLambda() );
88  if (cosl < 1.E-30) cosl = 1.E-30;
89  G4double cosl1 = 1./cosl;
90  G4Vector3D vUN(-vTN.y()*cosl1, vTN.x()*cosl1, 0. );
91  G4Vector3D vVN(-vTN.z()*vUN.y(), vTN.z()*vUN.x(), cosl );
92 
93  G4Vector3D vUperp = G4Vector3D( -fMomentum.y(), fMomentum.x(), 0.);
94  G4Vector3D vVperp = vUperp.cross( fMomentum );
95  vUperp *= 1./vUperp.mag();
96  vVperp *= 1./vVperp.mag();
97 
98 #ifdef G4EVERBOSE
99  if( iverbose >= 5 ){
100  G4cout << " CHECK: vUN " << vUN << " = " << vUperp << " diff " << (vUN-vUperp).mag() << G4endl;
101  G4cout << " CHECK: vVN " << vVN << " = " << vVperp << " diff " << (vVN-vVperp).mag() << G4endl;
102  }
103 #endif
104 
105  //get the dot products of vectors perpendicular to direction and vector defining SD plane
106  G4double dUU = vUperp * tpSD.GetVectorV();
107  G4double dUV = vUperp * tpSD.GetVectorW();
108  G4double dVU = vVperp * tpSD.GetVectorV();
109  G4double dVV = vVperp * tpSD.GetVectorW();
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  G4cout << " dir="<<dir<<" invCosTheta "<<invCosTheta << G4endl;
118 
119  if( fCharge != 0 && 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 }
tesla
static constexpr double tesla
Definition: G4SIunits.hh:268
G4ErrorTrajErr
G4ErrorSymMatrix G4ErrorTrajErr
Definition: G4ErrorTrajErr.hh:44
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Field::GetFieldValue
virtual void GetFieldValue(const double Point[4], double *fieldArr) const =0
CLHEP::Hep3Vector::x
double x() const
G4ErrorFreeTrajParam::GetPhi
G4double GetPhi() const
Definition: G4ErrorFreeTrajParam.hh:74
G4ErrorFreeTrajParam::GetLambda
G4double GetLambda() const
Definition: G4ErrorFreeTrajParam.hh:73
G4ErrorTrajState::fError
G4ErrorTrajErr fError
Definition: G4ErrorTrajState.hh:128
Q
static double Q[]
Definition: nf_gammaFunctions.cc:91
G4Vector3D
HepGeom::Vector3D< G4double > G4Vector3D
Definition: G4Vector3D.hh:35
CLHEP::Hep3Vector::z
double z() const
G4cout
G4GLOB_DLL std::ostream G4cout
G4ErrorFreeTrajParam::GetDirection
G4Vector3D GetDirection() const
Definition: G4ErrorFreeTrajParam.hh:70
cm
static constexpr double cm
Definition: G4SIunits.hh:119
CLHEP::Hep3Vector::theta
double theta() const
G4TransportationManager::GetTransportationManager
static G4TransportationManager * GetTransportationManager()
Definition: G4TransportationManager.cc:100
G4TransportationManager::GetFieldManager
G4FieldManager * GetFieldManager() const
CLHEP::c_light
static constexpr double c_light
Definition: PhysicalConstants.h:55
CLHEP::Hep3Vector::y
double y() const
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4FieldManager::GetDetectorField
const G4Field * GetDetectorField() const
CLHEP::Hep3Vector::mag
double mag() const
HepGeom::BasicVector3D::cross
BasicVector3D< T > cross(const BasicVector3D< T > &v) const
Definition: BasicVector3D.h:275
G4ErrorFreeTrajState::GetParameters
G4ErrorFreeTrajParam GetParameters() const
Definition: G4ErrorFreeTrajState.hh:108

Here is the call graph for this function:

G4ErrorFreeTrajState::~G4ErrorFreeTrajState ( )
inline

Definition at line 79 of file G4ErrorFreeTrajState.hh.

79 {}

Member Function Documentation

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

Implements G4ErrorTrajState.

Definition at line 170 of file G4ErrorFreeTrajState.cc.

171 {
172  out << *this;
173 }
G4ErrorFreeTrajParam G4ErrorFreeTrajState::GetParameters ( ) const
inline

Definition at line 108 of file G4ErrorFreeTrajState.hh.

109  { return fTrajParam; }

Here is the caller graph for this function:

G4ErrorMatrix G4ErrorFreeTrajState::GetTransfMat ( ) const
inline

Definition at line 111 of file G4ErrorFreeTrajState.hh.

112  { return theTransfMat; }
G4int G4ErrorFreeTrajState::PropagateError ( const G4Track *  aTrack)
virtual

Reimplemented from G4ErrorTrajState.

Definition at line 203 of file G4ErrorFreeTrajState.cc.

204 {
205  G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
206  if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetStage() == G4ErrorStage_Deflation ) stepLengthCm *= -1.;
207 
208  G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
209 
210  if( std::fabs(stepLengthCm) <= kCarTolerance/cm ) return 0;
211 
212 #ifdef G4EVERBOSE
213  if( iverbose >= 2 )G4cout << " G4ErrorFreeTrajState::PropagateError " << G4endl;
214  G4cout << "G4EP: iverbose="<< iverbose << G4endl;
215 #endif
216 
217  // * *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
218  G4Point3D vposPost = aTrack->GetPosition()/cm;
219  G4Vector3D vpPost = aTrack->GetMomentum()/GeV;
220  // G4Point3D vposPre = fPosition/cm;
221  // G4Vector3D vpPre = fMomentum/GeV;
222  G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition()/cm;
223  G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV;
224  //correct to avoid propagation along Z
225  if( vpPre.mag() == vpPre.z() ) vpPre.setX( 1.E-6*MeV );
226  if( vpPost.mag() == vpPost.z() ) vpPost.setX( 1.E-6*MeV );
227 
228  G4double pPre = vpPre.mag();
229  G4double pPost = vpPost.mag();
230 #ifdef G4EVERBOSE
231  if( iverbose >= 2 ) {
232  G4cout << "G4EP: vposPre " << vposPre << G4endl
233  << "G4EP: vposPost " << vposPost << G4endl;
234  G4cout << "G4EP: vpPre " << vpPre << G4endl
235  << "G4EP: vpPost " << vpPost << G4endl;
236  G4cout << " err start step " << fError << G4endl;
237  G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl;
238  }
239 #endif
240 
241  if( pPre == 0. || pPost == 0 ) return 2;
242  G4double pInvPre = 1./pPre;
243  G4double pInvPost = 1./pPost;
244  G4double deltaPInv = pInvPost - pInvPre;
245  if( iverbose >= 2 ) G4cout << "G4EP: pInvPre" << pInvPre<< " pInvPost:" << pInvPost<<" deltaPInv:" << deltaPInv<< G4endl;
246 
247 
248  G4Vector3D vpPreNorm = vpPre * pInvPre;
249  G4Vector3D vpPostNorm = vpPost * pInvPost;
250  if( iverbose >= 2 ) G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm << G4endl;
251  //return if propagation along Z??
252  if( 1. - std::fabs(vpPostNorm.z()) < kCarTolerance ) return 4;
253  G4double sinpPre = std::sin( vpPreNorm.theta() ); //cosine perpendicular to pPre = sine pPre
254  G4double sinpPost = std::sin( vpPostNorm.theta() ); //cosine perpendicular to pPost = sine pPost
255  G4double sinpPostInv = 1./std::sin( vpPreNorm.theta() );
256 
257 #ifdef G4EVERBOSE
258  if( iverbose >= 2 ) G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl;
259 #endif
260  //* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR
261  //* *** NEUTRAL PARTICLE OR FIELDFREE REGION
262  G4ErrorMatrix transf(5, 5, 0 );
263 
264  transf[3][2] = stepLengthCm * sinpPost;
265  transf[4][1] = stepLengthCm;
266  for( size_t ii=0;ii < 5; ii++ ){
267  transf[ii][ii] = 1.;
268  }
269 #ifdef G4EVERBOSE
270  if( iverbose >= 2 ) {
271  G4cout << "G4EP: transf matrix neutral " << transf;
272  }
273 #endif
274 
275  // charge X propagation direction
276  G4double charge = aTrack->GetDynamicParticle()->GetCharge();
277  if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() == G4ErrorMode_PropBackwards ) {
278  charge *= -1.;
279  }
280  // G4cout << " charge " << charge << G4endl;
281  //t check if particle has charge
282  //t if( charge == 0 ) goto 45;
283  // check if the magnetic field is = 0.
284 
285  //position is from geant4, it is assumed to be in mm (for debugging, eventually it will not be transformed)
286  //it is assumed vposPre[] is in cm and pos1[] is in mm.
287  G4double pos1[3]; pos1[0] = vposPre.x()*cm; pos1[1] = vposPre.y()*cm; pos1[2] = vposPre.z()*cm;
288  G4double pos2[3]; pos2[0] = vposPost.x()*cm; pos2[1] = vposPost.y()*cm; pos2[2] = vposPost.z()*cm;
289  G4double h1[3], h2[3];
290 
291  const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
292  if( !field ) return 0; //goto 45
293 
294 
295 
296  // calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION
297  if( charge != 0. && field ) {
298 
299  field->GetFieldValue( pos1, h1 ); //here pos1[], pos2[] are in mm, not changed
300  field->GetFieldValue( pos2, h2 );
301  G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.; //10. is to get same dimensions as GEANT3 (kilogauss)
302  G4ThreeVector HPost= G4ThreeVector( h2[0], h2[1], h2[2] ) / tesla *10.;
303  G4double magHPre = HPre.mag();
304  G4double magHPost = HPost.mag();
305 #ifdef G4EVERBOSE
306  if( iverbose >= 2 ) {
307  G4cout << "G4EP: h1 = "
308  << h1[0] << ", " << h1[1] << ", " << h1[2] << G4endl;
309  G4cout << "G4EP: pos1/mm = "
310  << pos1[0] << ", " << pos1[1] << ", " << pos1[2] << G4endl;
311  G4cout << "G4EP: pos2/mm = "
312  << pos2[0] << ", " << pos2[1] << ", " << pos2[2] << G4endl;
313  G4cout << "G4EP: B-filed in KGauss HPre " << HPre << G4endl
314  << "G4EP: in KGauss HPost " << HPost << G4endl;
315  }
316 #endif
317 
318  if( magHPre + magHPost != 0. ) {
319 
320  //* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2
321  G4double gam;
322  if( magHPost != 0. ){
323  gam = HPost * vpPostNorm / magHPost;
324  }else {
325  gam = HPre * vpPreNorm / magHPre;
326  }
327 
328  // G4eMagneticLimitsProcess will limit the step, but based on an straight line trajectory
329  G4double alphaSqr = 1. - gam * gam;
330  G4double diffHSqr = ( HPre * pInvPre - HPost * pInvPost ).mag2();
331  G4double delhp6Sqr = 300.*300.;
332 #ifdef G4EVERBOSE
333  if( iverbose >= 2 ) {
334  G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr
335  << " diffHSqr " << diffHSqr << G4endl;
336  G4cout << " alpha= " << std::sqrt(alphaSqr) << G4endl;
337  }
338 #endif
339  if( diffHSqr * alphaSqr > delhp6Sqr ) return 3;
340 
341 
342  //* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT
343  G4double pInvAver = 1./(pInvPre + pInvPost );
344  G4double CFACT8 = 2.997925E-4;
345  //G4double HAver
346  G4ThreeVector vHAverNorm( (HPre*pInvPre + HPost*pInvPost ) * pInvAver * charge * CFACT8 );
347  G4double HAver = vHAverNorm.mag();
348  G4double invHAver = 1./HAver;
349  vHAverNorm *= invHAver;
350 #ifdef G4EVERBOSE
351  if( iverbose >= 2 ) G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver << " charge " << charge<< G4endl;
352 #endif
353 
354  G4double pAver = (pPre+pPost)*0.5;
355  G4double QAver = -HAver/pAver;
356  G4double thetaAver = QAver * stepLengthCm;
357  G4double sinThetaAver = std::sin(thetaAver);
358  G4double cosThetaAver = std::cos(thetaAver);
359  G4double gamma = vHAverNorm * vpPostNorm;
360  G4ThreeVector AN2 = vHAverNorm.cross( vpPostNorm );
361 
362 #ifdef G4EVERBOSE
363  if( iverbose >= 2 ) G4cout << " G4EP: AN2 " << AN2 << " gamma:"<<gamma<< " theta="<< thetaAver<<G4endl;
364 #endif
365  G4double AU = 1./vpPreNorm.perp();
366  //t G4ThreeVector vU( vpPreNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU );
367  G4ThreeVector vUPre( -AU*vpPreNorm.y(),
368  AU*vpPreNorm.x(),
369  0. );
370  G4ThreeVector vVPre( -vpPreNorm.z()*vUPre.y(),
371  vpPreNorm.z()*vUPre.x(),
372  vpPreNorm.x()*vUPre.y() - vpPreNorm.y()*vUPre.x() );
373 
374  //
375  AU = 1./vpPostNorm.perp();
376  //t G4ThreeVector vU( vpPostNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU );
377  G4ThreeVector vUPost( -AU*vpPostNorm.y(),
378  AU*vpPostNorm.x(),
379  0. );
380  G4ThreeVector vVPost( -vpPostNorm.z()*vUPost.y(),
381  vpPostNorm.z()*vUPost.x(),
382  vpPostNorm.x()*vUPost.y() - vpPostNorm.y()*vUPost.x() );
383 #ifdef G4EVERBOSE
384  G4cout << " vpPostNorm " << vpPostNorm << G4endl;
385  if( iverbose >= 2 ) G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre << " vUPost " << vUPost << " vVPost " << vVPost << G4endl;
386 #endif
387  G4Point3D deltaPos( vposPre - vposPost );
388 
389  // * *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2
390  // * *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT
391  // * *** TAKEN INTO ACCOUNT
392 
393  G4double QP = QAver * pAver; // = -HAver
394 #ifdef G4EVERBOSE
395  if( iverbose >= 2) G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver << G4endl;
396 #endif
397  G4double ANV = -( vHAverNorm.x()*vUPost.x() + vHAverNorm.y()*vUPost.y() );
398  G4double ANU = ( vHAverNorm.x()*vVPost.x() + vHAverNorm.y()*vVPost.y() + vHAverNorm.z()*vVPost.z() );
399  G4double OMcosThetaAver = 1. - cosThetaAver;
400 #ifdef G4EVERBOSE
401  if( iverbose >= 2) G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver " << cosThetaAver << " thetaAver " << thetaAver << " QAver " << QAver << " stepLengthCm " << stepLengthCm << G4endl;
402 #endif
403  G4double TMSINT = thetaAver - sinThetaAver;
404 #ifdef G4EVERBOSE
405  if( iverbose >= 2 ) G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl;
406 #endif
407 
408  G4ThreeVector vHUPre( -vHAverNorm.z() * vUPre.y(),
409  vHAverNorm.z() * vUPre.x(),
410  vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x() );
411 #ifdef G4EVERBOSE
412  // if( iverbose >= 2 ) G4cout << "G4EP: HUPre(1) " << vHUPre.x() << " " << vHAverNorm.z() << " " << vUPre.y() << G4endl;
413 #endif
414  G4ThreeVector vHVPre( vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(),
415  vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(),
416  vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x() );
417 #ifdef G4EVERBOSE
418  if( iverbose >= 2 ) G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl;
419 #endif
420 
421  //------------------- COMPUTE MATRIX
422  //---------- 1/P
423 
424  transf[0][0] = 1.-deltaPInv*pAver*(1.+(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())/stepLengthCm)
425  +2.*deltaPInv*pAver;
426 
427  transf[0][1] = -deltaPInv/thetaAver*
428  ( TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) +
429  sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
430  OMcosThetaAver*(vHVPre.x()*vpPostNorm.x()+vHVPre.y()*vpPostNorm.y()+vHVPre.z()*vpPostNorm.z()) );
431 
432  transf[0][2] = -sinpPre*deltaPInv/thetaAver*
433  ( TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) +
434  sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) +
435  OMcosThetaAver*(vHUPre.x()*vpPostNorm.x()+vHUPre.y()*vpPostNorm.y()+vHUPre.z()*vpPostNorm.z()) );
436 
437  transf[0][3] = -deltaPInv/stepLengthCm*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() );
438 
439  transf[0][4] = -deltaPInv/stepLengthCm*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z());
440 
441  // *** Lambda
442  transf[1][0] = -QP*ANV*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())
443  *(1.+deltaPInv*pAver);
444 #ifdef G4EVERBOSE
445  if(iverbose >= 3) G4cout << "ctransf10= " << transf[1][0] << " " << -QP<< " " << ANV<< " " << vpPostNorm.x()<< " " << deltaPos.x()<< " " << vpPostNorm.y()<< " " << deltaPos.y()<< " " << vpPostNorm.z()<< " " << deltaPos.z()
446  << " " << deltaPInv<< " " << pAver << G4endl;
447 #endif
448 
449  transf[1][1] = cosThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) +
450  sinThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) +
451  OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())*
452  (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) +
453  ANV*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
454  OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) -
455  TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) );
456 
457  transf[1][2] = cosThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) +
458  sinThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) +
459  OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )*
460  (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) +
461  ANV*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) +
462  OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) -
463  TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) );
464  transf[1][2] = sinpPre*transf[1][2];
465 
466  transf[1][3] = -QAver*ANV*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() );
467 
468  transf[1][4] = -QAver*ANV*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z());
469 
470  // *** Phi
471 
472  transf[2][0] = -QP*ANU*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())*sinpPostInv
473  *(1.+deltaPInv*pAver);
474 #ifdef G4EVERBOSE
475  if(iverbose >= 3)G4cout <<"ctransf20= " << transf[2][0] <<" "<< -QP<<" "<<ANU<<" "<<vpPostNorm.x()<<" "<<deltaPos.x()<<" "<<vpPostNorm.y()<<" "<<deltaPos.y()<<" "<<vpPostNorm.z()<<" "<<deltaPos.z()<<" "<<sinpPostInv
476  <<" "<<deltaPInv<<" "<<pAver<< G4endl;
477 #endif
478  transf[2][1] = cosThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) +
479  sinThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) +
480  OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())*
481  (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) +
482  ANU*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
483  OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) -
484  TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) );
485  transf[2][1] = sinpPostInv*transf[2][1];
486 
487  transf[2][2] = cosThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) +
488  sinThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) +
489  OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )*
490  (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) +
491  ANU*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) +
492  OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) -
493  TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) );
494  transf[2][2] = sinpPostInv*sinpPre*transf[2][2];
495 
496  transf[2][3] = -QAver*ANU*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() )*sinpPostInv;
497 #ifdef G4EVERBOSE
498  if(iverbose >= 3)G4cout <<"ctransf23= " << transf[2][3] <<" "<< -QAver<<" "<<ANU<<" "<<vUPre.x()<<" "<<vpPostNorm.x()<<" "<< vUPre.y()<<" "<<vpPostNorm.y()<<" "<<sinpPostInv<<G4endl;
499 #endif
500 
501  transf[2][4] = -QAver*ANU*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z())*sinpPostInv;
502 
503  // *** Yt
504 
505  transf[3][0] = pAver*(vUPost.x()*deltaPos.x()+vUPost.y()*deltaPos.y() )
506  *(1.+deltaPInv*pAver);
507 #ifdef G4EVERBOSE
508  if(iverbose >= 3) G4cout <<"ctransf30= " << transf[3][0] <<" "<< pAver<<" "<<vUPost.x()<<" "<<deltaPos.x()<<" "<<vUPost.y()<<" "<<deltaPos.y()
509  <<" "<<deltaPInv<<" "<<pAver<<G4endl;
510 #endif
511 
512  transf[3][1] = ( sinThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) +
513  OMcosThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) +
514  TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )*
515  (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver;
516 
517  transf[3][2] = ( sinThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) +
518  OMcosThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) +
519  TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )*
520  (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver;
521 #ifdef G4EVERBOSE
522  if(iverbose >= 3) G4cout <<"ctransf32= " << transf[3][2] <<" "<< sinThetaAver<<" "<<vUPre.x()<<" "<<vUPost.x()<<" "<<vUPre.y()<<" "<<vUPost.y() <<" "<<
523  OMcosThetaAver<<" "<<vHUPre.x()<<" "<<vUPost.x()<<" "<<vHUPre.y()<<" "<<vUPost.y() <<" "<<
524  TMSINT<<" "<<vHAverNorm.x()<<" "<<vUPost.x()<<" "<<vHAverNorm.y()<<" "<<vUPost.y() <<" "<<
525  vHAverNorm.x()<<" "<<vUPre.x()<<" "<<vHAverNorm.y()<<" "<<vUPre.y() <<" "<<sinpPre<<" "<<QAver<<G4endl;
526 #endif
527 
528  transf[3][3] = (vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() );
529 
530  transf[3][4] = (vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() );
531 
532  // *** Zt
533  transf[4][0] = pAver*(vVPost.x()*deltaPos.x()+vVPost.y()*deltaPos.y()+vVPost.z()*deltaPos.z())
534  *(1.+deltaPInv*pAver);
535 
536  transf[4][1] = ( sinThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) +
537  OMcosThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) +
538  TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())*
539  (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver;
540 #ifdef G4EVERBOSE
541  if(iverbose >= 3)G4cout <<"ctransf41= " << transf[4][1] <<" "<< sinThetaAver<<" "<< OMcosThetaAver <<" "<<TMSINT<<" "<< vVPre <<" "<<vVPost <<" "<<vHVPre<<" "<<vHAverNorm <<" "<< QAver<<G4endl;
542 #endif
543 
544  transf[4][2] = ( sinThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) +
545  OMcosThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) +
546  TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())*
547  (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver;
548 
549  transf[4][3] = (vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() );
550 
551  transf[4][4] = (vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z());
552  // if(iverbose >= 3) G4cout <<"ctransf44= " << transf[4][4] <<" "<< vVPre.x() <<" "<<vVPost.x() <<" "<< vVPre.y() <<" "<< vVPost.y() <<" "<< vVPre.z() <<" "<< vVPost.z() << G4endl;
553 
554 
555 #ifdef G4EVERBOSE
556  if( iverbose >= 1 ) G4cout << "G4EP: transf matrix computed " << transf << G4endl;
557 #endif
558  /* for( G4int ii=0;ii<5;ii++){
559  for( G4int jj=0;jj<5;jj++){
560  G4cout << transf[ii][jj] << " ";
561  }
562  G4cout << G4endl;
563  } */
564  }
565  }
566  // end of calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION
567  /* if( iverbose >= 1 ) G4cout << "G4EP: transf not updated but initialized " << theFirstStep << G4endl;
568  if( theFirstStep ) {
569  theTransfMat = transf;
570  theFirstStep = false;
571  }else{
572  theTransfMat = theTransfMat * transf;
573  if( iverbose >= 1 ) G4cout << "G4EP: transf matrix accumulated" << theTransfMat << G4endl;
574  }
575  */
576  theTransfMat = transf;
577 #ifdef G4EVERBOSE
578  if( iverbose >= 1 ) G4cout << "G4EP: error matrix before transformation " << fError << G4endl;
579  if( iverbose >= 2 ) G4cout << " tf * err " << theTransfMat * fError << G4endl
580  << " transf matrix " << theTransfMat.T() << G4endl;
581 #endif
582 
583  fError = fError.similarity(theTransfMat).T();
584  //- fError = transf * fError * transf.T();
585 #ifdef G4EVERBOSE
586  if( iverbose >= 1 ) G4cout << "G4EP: error matrix propagated " << fError << G4endl;
587 #endif
588 
589  //? S = B*S*BT S.similarity(B)
590  //? R = S
591  // not needed * *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL VARIABLES;
592 
593  PropagateErrorMSC( aTrack );
594 
595  PropagateErrorIoni( aTrack );
596 
597  return 0;
598 }
tesla
static constexpr double tesla
Definition: G4SIunits.hh:268
G4ErrorSymMatrix::similarity
G4ErrorSymMatrix similarity(const G4ErrorMatrix &m1) const
Definition: G4ErrorSymMatrix.cc:589
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4Step::GetStepLength
G4double GetStepLength() const
G4Field::GetFieldValue
virtual void GetFieldValue(const double Point[4], double *fieldArr) const =0
CLHEP::Hep3Vector::x
double x() const
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4Track::GetPosition
const G4ThreeVector & GetPosition() const
G4GeometryTolerance::GetSurfaceTolerance
G4double GetSurfaceTolerance() const
Definition: G4GeometryTolerance.cc:87
G4ErrorTrajState::fError
G4ErrorTrajErr fError
Definition: G4ErrorTrajState.hh:128
G4StepPoint::GetMomentum
G4ThreeVector GetMomentum() const
G4ErrorSymMatrix::T
G4ErrorSymMatrix T() const
G4Track::GetStep
const G4Step * GetStep() const
CLHEP::Hep3Vector::z
double z() const
G4Step::GetPreStepPoint
G4StepPoint * GetPreStepPoint() const
G4ErrorMode_PropBackwards
Definition: G4ErrorPropagatorData.hh:48
G4cout
G4GLOB_DLL std::ostream G4cout
G4StepPoint::GetPosition
const G4ThreeVector & GetPosition() const
cm
static constexpr double cm
Definition: G4SIunits.hh:119
G4DynamicParticle::GetCharge
G4double GetCharge() const
kCarTolerance
const G4double kCarTolerance
Definition: G4BoundingEnvelope.cc:44
G4TransportationManager::GetTransportationManager
static G4TransportationManager * GetTransportationManager()
Definition: G4TransportationManager.cc:100
G4TransportationManager::GetFieldManager
G4FieldManager * GetFieldManager() const
G4ErrorMatrix::T
G4ErrorMatrix T() const
Definition: G4ErrorMatrix.cc:374
G4Track::GetMomentum
G4ThreeVector GetMomentum() const
CLHEP::Hep3Vector::y
double y() const
G4ErrorStage_Deflation
Definition: G4ErrorPropagatorData.hh:58
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
G4double
double G4double
Definition: G4Types.hh:76
G4ErrorPropagatorData::GetErrorPropagatorData
static G4ErrorPropagatorData * GetErrorPropagatorData()
Definition: G4ErrorPropagatorData.cc:54
G4FieldManager::GetDetectorField
const G4Field * GetDetectorField() const
CLHEP::Hep3Vector::perp
double perp() const
CLHEP::Hep3Vector::mag
double mag() const
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:73

Here is the call graph for this function:

Here is the caller graph for this function:

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

Reimplemented from G4ErrorTrajState.

Definition at line 98 of file G4ErrorFreeTrajState.hh.

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

Here is the call graph for this function:

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

Definition at line 101 of file G4ErrorFreeTrajState.hh.

102  {
103  fPosition = pos;
104  fMomentum = mom;
105  fTrajParam.SetParameters( pos, mom );
106  }
G4ErrorFreeTrajParam::SetParameters
void SetParameters(const G4Point3D &pos, const G4Vector3D &mom)
Definition: G4ErrorFreeTrajParam.cc:48
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66

Here is the call graph for this function:

Here is the caller graph for this function:

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

Reimplemented from G4ErrorTrajState.

Definition at line 95 of file G4ErrorFreeTrajState.hh.

96  { SetParameters( pos, fMomentum ); }
G4ErrorFreeTrajState::SetParameters
void SetParameters(const G4Point3D &pos, const G4Vector3D &mom)
Definition: G4ErrorFreeTrajState.hh:101

Here is the call graph for this function:

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

Reimplemented from G4ErrorTrajState.

Definition at line 176 of file G4ErrorFreeTrajState.cc.

177 {
178  G4int ierr = 0;
179  fTrajParam.Update( aTrack );
180  UpdatePosMom( aTrack->GetPosition(), aTrack->GetMomentum() );
181  return ierr;
182 }
G4Track::GetPosition
const G4ThreeVector & GetPosition() const
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorFreeTrajParam::Update
void Update(const G4Track *aTrack)
Definition: G4ErrorFreeTrajParam.cc:79
G4Track::GetMomentum
G4ThreeVector GetMomentum() const
G4ErrorTrajState::UpdatePosMom
void UpdatePosMom(const G4Point3D &pos, const G4Vector3D &mom)
Definition: G4ErrorTrajState.cc:65

Here is the call graph for this function:

Here is the caller graph for this function:

Friends And Related Function Documentation

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

Definition at line 186 of file G4ErrorFreeTrajState.cc.

187 {
188  std::ios::fmtflags orig_flags = out.flags();
189 
190  out.setf(std::ios::fixed,std::ios::floatfield);
191 
192  ts.DumpPosMomError( out );
193 
194  out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl;
195 
196  out.flags(orig_flags);
197 
198  return out;
199 }
G4ErrorTrajState::DumpPosMomError
void DumpPosMomError(std::ostream &out=G4cout) const
Definition: G4ErrorTrajState.cc:104
G4endl
#define G4endl
Definition: G4ios.hh:61
The documentation for this class was generated from the following files: