Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ErrorFreeTrajState.cc
이 파일의 문서화 페이지로 가기
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 107409 2017-11-10 13:36:41Z gcosmo $
27 //
28 // ------------------------------------------------------------
29 // GEANT 4 class implementation file
30 // ------------------------------------------------------------
31 //
32 
33 #include <iomanip>
34 
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4Field.hh"
38 #include "G4FieldManager.hh"
40 #include "G4GeometryTolerance.hh"
41 #include "G4Material.hh"
42 #include "G4ErrorPropagatorData.hh"
43 #include "G4ErrorFreeTrajState.hh"
44 #include "G4ErrorFreeTrajParam.hh"
46 #include "G4ErrorMatrix.hh"
47 
48 //------------------------------------------------------------------------
49 G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4String& partName, const G4Point3D& pos, const G4Vector3D& mom, const G4ErrorTrajErr& errmat) : G4ErrorTrajState( partName, pos, mom, errmat )
50 {
51  fTrajParam = G4ErrorFreeTrajParam( pos, mom );
52  Init();
53 }
54 
55 
56 //------------------------------------------------------------------------
57 G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4ErrorSurfaceTrajState& tpSD ) : 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 
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
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 }
158 
159 
160 //------------------------------------------------------------------------
162 {
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 //------------------------------------------------------------------------
177 {
178  G4int ierr = 0;
179  fTrajParam.Update( aTrack );
180  UpdatePosMom( aTrack->GetPosition(), aTrack->GetMomentum() );
181  return ierr;
182 }
183 
184 
185 //------------------------------------------------------------------------
186 std::ostream& operator<<(std::ostream& out, const G4ErrorFreeTrajState& ts)
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 }
200 
201 
202 //------------------------------------------------------------------------
204 {
205  G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
206  if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetStage() == G4ErrorStage_Deflation ) stepLengthCm *= -1.;
207 
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(vpPreNorm.z()) < kCarTolerance ) return 4;
253  if( 1. - std::fabs(vpPostNorm.z()) < kCarTolerance ) return 4;
254  G4double sinpPre = std::sin( vpPreNorm.theta() ); //cosine perpendicular to pPre = sine pPre
255  G4double sinpPost = std::sin( vpPostNorm.theta() ); //cosine perpendicular to pPost = sine pPost
256  G4double sinpPostInv = 1./std::sin( vpPostNorm.theta() );
257 
258 #ifdef G4EVERBOSE
259  if( iverbose >= 2 ) G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl;
260 #endif
261  //* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR
262  //* *** NEUTRAL PARTICLE OR FIELDFREE REGION
263  G4ErrorMatrix transf(5, 5, 0 );
264 
265  transf[3][2] = stepLengthCm * sinpPost;
266  transf[4][1] = stepLengthCm;
267  for( size_t ii=0;ii < 5; ii++ ){
268  transf[ii][ii] = 1.;
269  }
270 #ifdef G4EVERBOSE
271  if( iverbose >= 2 ) {
272  G4cout << "G4EP: transf matrix neutral " << transf;
273  }
274 #endif
275 
276  // charge X propagation direction
277  G4double charge = aTrack->GetDynamicParticle()->GetCharge();
279  charge *= -1.;
280  }
281  // G4cout << " charge " << charge << G4endl;
282  //t check if particle has charge
283  //t if( charge == 0 ) goto 45;
284  // check if the magnetic field is = 0.
285 
286  //position is from geant4, it is assumed to be in mm (for debugging, eventually it will not be transformed)
287  //it is assumed vposPre[] is in cm and pos1[] is in mm.
288  G4double pos1[3]; pos1[0] = vposPre.x()*cm; pos1[1] = vposPre.y()*cm; pos1[2] = vposPre.z()*cm;
289  G4double pos2[3]; pos2[0] = vposPost.x()*cm; pos2[1] = vposPost.y()*cm; pos2[2] = vposPost.z()*cm;
290  G4double h1[3], h2[3];
291 
293  if( !field ) return 0; //goto 45
294 
295 
296 
297  // calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION
298  if( charge != 0. && field ) {
299 
300  field->GetFieldValue( pos1, h1 ); //here pos1[], pos2[] are in mm, not changed
301  field->GetFieldValue( pos2, h2 );
302  G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.; //10. is to get same dimensions as GEANT3 (kilogauss)
303  G4ThreeVector HPost= G4ThreeVector( h2[0], h2[1], h2[2] ) / tesla *10.;
304  G4double magHPre = HPre.mag();
305  G4double magHPost = HPost.mag();
306 #ifdef G4EVERBOSE
307  if( iverbose >= 2 ) {
308  G4cout << "G4EP: h1 = "
309  << h1[0] << ", " << h1[1] << ", " << h1[2] << G4endl;
310  G4cout << "G4EP: pos1/mm = "
311  << pos1[0] << ", " << pos1[1] << ", " << pos1[2] << G4endl;
312  G4cout << "G4EP: pos2/mm = "
313  << pos2[0] << ", " << pos2[1] << ", " << pos2[2] << G4endl;
314  G4cout << "G4EP: B-filed in KGauss HPre " << HPre << G4endl
315  << "G4EP: in KGauss HPost " << HPost << G4endl;
316  }
317 #endif
318 
319  if( magHPre + magHPost != 0. ) {
320 
321  //* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2
322  G4double gam;
323  if( magHPost != 0. ){
324  gam = HPost * vpPostNorm / magHPost;
325  }else {
326  gam = HPre * vpPreNorm / magHPre;
327  }
328 
329  // G4eMagneticLimitsProcess will limit the step, but based on an straight line trajectory
330  G4double alphaSqr = 1. - gam * gam;
331  G4double diffHSqr = ( HPre * pInvPre - HPost * pInvPost ).mag2();
332  G4double delhp6Sqr = 300.*300.;
333 #ifdef G4EVERBOSE
334  if( iverbose >= 2 ) {
335  G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr
336  << " diffHSqr " << diffHSqr << G4endl;
337  G4cout << " alpha= " << std::sqrt(alphaSqr) << G4endl;
338  }
339 #endif
340  if( diffHSqr * alphaSqr > delhp6Sqr ) return 3;
341 
342 
343  //* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT
344  G4double pInvAver = 1./(pInvPre + pInvPost );
345  G4double CFACT8 = 2.997925E-4;
346  //G4double HAver
347  G4ThreeVector vHAverNorm( (HPre*pInvPre + HPost*pInvPost ) * pInvAver * charge * CFACT8 );
348  G4double HAver = vHAverNorm.mag();
349  G4double invHAver = 1./HAver;
350  vHAverNorm *= invHAver;
351 #ifdef G4EVERBOSE
352  if( iverbose >= 2 ) G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver << " charge " << charge<< G4endl;
353 #endif
354 
355  G4double pAver = (pPre+pPost)*0.5;
356  G4double QAver = -HAver/pAver;
357  G4double thetaAver = QAver * stepLengthCm;
358  G4double sinThetaAver = std::sin(thetaAver);
359  G4double cosThetaAver = std::cos(thetaAver);
360  G4double gamma = vHAverNorm * vpPostNorm;
361  G4ThreeVector AN2 = vHAverNorm.cross( vpPostNorm );
362 
363 #ifdef G4EVERBOSE
364  if( iverbose >= 2 ) G4cout << " G4EP: AN2 " << AN2 << " gamma:"<<gamma<< " theta="<< thetaAver<<G4endl;
365 #endif
366  G4double AU = 1./vpPreNorm.perp();
367  //t G4ThreeVector vU( vpPreNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU );
368  G4ThreeVector vUPre( -AU*vpPreNorm.y(),
369  AU*vpPreNorm.x(),
370  0. );
371  G4ThreeVector vVPre( -vpPreNorm.z()*vUPre.y(),
372  vpPreNorm.z()*vUPre.x(),
373  vpPreNorm.x()*vUPre.y() - vpPreNorm.y()*vUPre.x() );
374 
375  //
376  AU = 1./vpPostNorm.perp();
377  //t G4ThreeVector vU( vpPostNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU );
378  G4ThreeVector vUPost( -AU*vpPostNorm.y(),
379  AU*vpPostNorm.x(),
380  0. );
381  G4ThreeVector vVPost( -vpPostNorm.z()*vUPost.y(),
382  vpPostNorm.z()*vUPost.x(),
383  vpPostNorm.x()*vUPost.y() - vpPostNorm.y()*vUPost.x() );
384 #ifdef G4EVERBOSE
385  G4cout << " vpPostNorm " << vpPostNorm << G4endl;
386  if( iverbose >= 2 ) G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre << " vUPost " << vUPost << " vVPost " << vVPost << G4endl;
387 #endif
388  G4Point3D deltaPos( vposPre - vposPost );
389 
390  // * *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2
391  // * *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT
392  // * *** TAKEN INTO ACCOUNT
393 
394  G4double QP = QAver * pAver; // = -HAver
395 #ifdef G4EVERBOSE
396  if( iverbose >= 2) G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver << G4endl;
397 #endif
398  G4double ANV = -( vHAverNorm.x()*vUPost.x() + vHAverNorm.y()*vUPost.y() );
399  G4double ANU = ( vHAverNorm.x()*vVPost.x() + vHAverNorm.y()*vVPost.y() + vHAverNorm.z()*vVPost.z() );
400  G4double OMcosThetaAver = 1. - cosThetaAver;
401 #ifdef G4EVERBOSE
402  if( iverbose >= 2) G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver " << cosThetaAver << " thetaAver " << thetaAver << " QAver " << QAver << " stepLengthCm " << stepLengthCm << G4endl;
403 #endif
404  G4double TMSINT = thetaAver - sinThetaAver;
405 #ifdef G4EVERBOSE
406  if( iverbose >= 2 ) G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl;
407 #endif
408 
409  G4ThreeVector vHUPre( -vHAverNorm.z() * vUPre.y(),
410  vHAverNorm.z() * vUPre.x(),
411  vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x() );
412 #ifdef G4EVERBOSE
413  // if( iverbose >= 2 ) G4cout << "G4EP: HUPre(1) " << vHUPre.x() << " " << vHAverNorm.z() << " " << vUPre.y() << G4endl;
414 #endif
415  G4ThreeVector vHVPre( vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(),
416  vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(),
417  vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x() );
418 #ifdef G4EVERBOSE
419  if( iverbose >= 2 ) G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl;
420 #endif
421 
422  //------------------- COMPUTE MATRIX
423  //---------- 1/P
424 
425  transf[0][0] = 1.-deltaPInv*pAver*(1.+(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())/stepLengthCm)
426  +2.*deltaPInv*pAver;
427 
428  transf[0][1] = -deltaPInv/thetaAver*
429  ( TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) +
430  sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
431  OMcosThetaAver*(vHVPre.x()*vpPostNorm.x()+vHVPre.y()*vpPostNorm.y()+vHVPre.z()*vpPostNorm.z()) );
432 
433  transf[0][2] = -sinpPre*deltaPInv/thetaAver*
434  ( TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) +
435  sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) +
436  OMcosThetaAver*(vHUPre.x()*vpPostNorm.x()+vHUPre.y()*vpPostNorm.y()+vHUPre.z()*vpPostNorm.z()) );
437 
438  transf[0][3] = -deltaPInv/stepLengthCm*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() );
439 
440  transf[0][4] = -deltaPInv/stepLengthCm*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z());
441 
442  // *** Lambda
443  transf[1][0] = -QP*ANV*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())
444  *(1.+deltaPInv*pAver);
445 #ifdef G4EVERBOSE
446  if(iverbose >= 3) G4cout << "ctransf10= " << transf[1][0] << " " << -QP<< " " << ANV<< " " << vpPostNorm.x()<< " " << deltaPos.x()<< " " << vpPostNorm.y()<< " " << deltaPos.y()<< " " << vpPostNorm.z()<< " " << deltaPos.z()
447  << " " << deltaPInv<< " " << pAver << G4endl;
448 #endif
449 
450  transf[1][1] = cosThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) +
451  sinThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) +
452  OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())*
453  (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) +
454  ANV*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
455  OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) -
456  TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) );
457 
458  transf[1][2] = cosThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) +
459  sinThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) +
460  OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )*
461  (vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z()) +
462  ANV*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) +
463  OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) -
464  TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) );
465  transf[1][2] = sinpPre*transf[1][2];
466 
467  transf[1][3] = -QAver*ANV*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() );
468 
469  transf[1][4] = -QAver*ANV*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z());
470 
471  // *** Phi
472 
473  transf[2][0] = -QP*ANU*(vpPostNorm.x()*deltaPos.x()+vpPostNorm.y()*deltaPos.y()+vpPostNorm.z()*deltaPos.z())*sinpPostInv
474  *(1.+deltaPInv*pAver);
475 #ifdef G4EVERBOSE
476  if(iverbose >= 3)G4cout <<"ctransf20= " << transf[2][0] <<" "<< -QP<<" "<<ANU<<" "<<vpPostNorm.x()<<" "<<deltaPos.x()<<" "<<vpPostNorm.y()<<" "<<deltaPos.y()<<" "<<vpPostNorm.z()<<" "<<deltaPos.z()<<" "<<sinpPostInv
477  <<" "<<deltaPInv<<" "<<pAver<< G4endl;
478 #endif
479  transf[2][1] = cosThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) +
480  sinThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) +
481  OMcosThetaAver*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z())*
482  (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) +
483  ANU*( -sinThetaAver*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z()) +
484  OMcosThetaAver*(vVPre.x()*AN2.x()+vVPre.y()*AN2.y()+vVPre.z()*AN2.z()) -
485  TMSINT*gamma*(vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) );
486  transf[2][1] = sinpPostInv*transf[2][1];
487 
488  transf[2][2] = cosThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) +
489  sinThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) +
490  OMcosThetaAver*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() )*
491  (vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() ) +
492  ANU*( -sinThetaAver*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() ) +
493  OMcosThetaAver*(vUPre.x()*AN2.x()+vUPre.y()*AN2.y() ) -
494  TMSINT*gamma*(vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) );
495  transf[2][2] = sinpPostInv*sinpPre*transf[2][2];
496 
497  transf[2][3] = -QAver*ANU*(vUPre.x()*vpPostNorm.x()+vUPre.y()*vpPostNorm.y() )*sinpPostInv;
498 #ifdef G4EVERBOSE
499  if(iverbose >= 3)G4cout <<"ctransf23= " << transf[2][3] <<" "<< -QAver<<" "<<ANU<<" "<<vUPre.x()<<" "<<vpPostNorm.x()<<" "<< vUPre.y()<<" "<<vpPostNorm.y()<<" "<<sinpPostInv<<G4endl;
500 #endif
501 
502  transf[2][4] = -QAver*ANU*(vVPre.x()*vpPostNorm.x()+vVPre.y()*vpPostNorm.y()+vVPre.z()*vpPostNorm.z())*sinpPostInv;
503 
504  // *** Yt
505 
506  transf[3][0] = pAver*(vUPost.x()*deltaPos.x()+vUPost.y()*deltaPos.y() )
507  *(1.+deltaPInv*pAver);
508 #ifdef G4EVERBOSE
509  if(iverbose >= 3) G4cout <<"ctransf30= " << transf[3][0] <<" "<< pAver<<" "<<vUPost.x()<<" "<<deltaPos.x()<<" "<<vUPost.y()<<" "<<deltaPos.y()
510  <<" "<<deltaPInv<<" "<<pAver<<G4endl;
511 #endif
512 
513  transf[3][1] = ( sinThetaAver*(vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() ) +
514  OMcosThetaAver*(vHVPre.x()*vUPost.x()+vHVPre.y()*vUPost.y() ) +
515  TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )*
516  (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver;
517 
518  transf[3][2] = ( sinThetaAver*(vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() ) +
519  OMcosThetaAver*(vHUPre.x()*vUPost.x()+vHUPre.y()*vUPost.y() ) +
520  TMSINT*(vHAverNorm.x()*vUPost.x()+vHAverNorm.y()*vUPost.y() )*
521  (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver;
522 #ifdef G4EVERBOSE
523  if(iverbose >= 3) G4cout <<"ctransf32= " << transf[3][2] <<" "<< sinThetaAver<<" "<<vUPre.x()<<" "<<vUPost.x()<<" "<<vUPre.y()<<" "<<vUPost.y() <<" "<<
524  OMcosThetaAver<<" "<<vHUPre.x()<<" "<<vUPost.x()<<" "<<vHUPre.y()<<" "<<vUPost.y() <<" "<<
525  TMSINT<<" "<<vHAverNorm.x()<<" "<<vUPost.x()<<" "<<vHAverNorm.y()<<" "<<vUPost.y() <<" "<<
526  vHAverNorm.x()<<" "<<vUPre.x()<<" "<<vHAverNorm.y()<<" "<<vUPre.y() <<" "<<sinpPre<<" "<<QAver<<G4endl;
527 #endif
528 
529  transf[3][3] = (vUPre.x()*vUPost.x()+vUPre.y()*vUPost.y() );
530 
531  transf[3][4] = (vVPre.x()*vUPost.x()+vVPre.y()*vUPost.y() );
532 
533  // *** Zt
534  transf[4][0] = pAver*(vVPost.x()*deltaPos.x()+vVPost.y()*deltaPos.y()+vVPost.z()*deltaPos.z())
535  *(1.+deltaPInv*pAver);
536 
537  transf[4][1] = ( sinThetaAver*(vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z()) +
538  OMcosThetaAver*(vHVPre.x()*vVPost.x()+vHVPre.y()*vVPost.y()+vHVPre.z()*vVPost.z()) +
539  TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())*
540  (vHAverNorm.x()*vVPre.x()+vHAverNorm.y()*vVPre.y()+vHAverNorm.z()*vVPre.z()) )/QAver;
541 #ifdef G4EVERBOSE
542  if(iverbose >= 3)G4cout <<"ctransf41= " << transf[4][1] <<" "<< sinThetaAver<<" "<< OMcosThetaAver <<" "<<TMSINT<<" "<< vVPre <<" "<<vVPost <<" "<<vHVPre<<" "<<vHAverNorm <<" "<< QAver<<G4endl;
543 #endif
544 
545  transf[4][2] = ( sinThetaAver*(vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() ) +
546  OMcosThetaAver*(vHUPre.x()*vVPost.x()+vHUPre.y()*vVPost.y()+vHUPre.z()*vVPost.z()) +
547  TMSINT*(vHAverNorm.x()*vVPost.x()+vHAverNorm.y()*vVPost.y()+vHAverNorm.z()*vVPost.z())*
548  (vHAverNorm.x()*vUPre.x()+vHAverNorm.y()*vUPre.y() ) )*sinpPre/QAver;
549 
550  transf[4][3] = (vUPre.x()*vVPost.x()+vUPre.y()*vVPost.y() );
551 
552  transf[4][4] = (vVPre.x()*vVPost.x()+vVPre.y()*vVPost.y()+vVPre.z()*vVPost.z());
553  // if(iverbose >= 3) G4cout <<"ctransf44= " << transf[4][4] <<" "<< vVPre.x() <<" "<<vVPost.x() <<" "<< vVPre.y() <<" "<< vVPost.y() <<" "<< vVPre.z() <<" "<< vVPost.z() << G4endl;
554 
555 
556 #ifdef G4EVERBOSE
557  if( iverbose >= 1 ) G4cout << "G4EP: transf matrix computed " << transf << G4endl;
558 #endif
559  /* for( G4int ii=0;ii<5;ii++){
560  for( G4int jj=0;jj<5;jj++){
561  G4cout << transf[ii][jj] << " ";
562  }
563  G4cout << G4endl;
564  } */
565  }
566  }
567  // end of calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION
568  /* if( iverbose >= 1 ) G4cout << "G4EP: transf not updated but initialized " << theFirstStep << G4endl;
569  if( theFirstStep ) {
570  theTransfMat = transf;
571  theFirstStep = false;
572  }else{
573  theTransfMat = theTransfMat * transf;
574  if( iverbose >= 1 ) G4cout << "G4EP: transf matrix accumulated" << theTransfMat << G4endl;
575  }
576  */
577  theTransfMat = transf;
578 #ifdef G4EVERBOSE
579  if( iverbose >= 1 ) G4cout << "G4EP: error matrix before transformation " << fError << G4endl;
580  if( iverbose >= 2 ) G4cout << " tf * err " << theTransfMat * fError << G4endl
581  << " transf matrix " << theTransfMat.T() << G4endl;
582 #endif
583 
585  //- fError = transf * fError * transf.T();
586 #ifdef G4EVERBOSE
587  if( iverbose >= 1 ) G4cout << "G4EP: error matrix propagated " << fError << G4endl;
588 #endif
589 
590  //? S = B*S*BT S.similarity(B)
591  //? R = S
592  // not needed * *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL VARIABLES;
593 
594  PropagateErrorMSC( aTrack );
595 
596  PropagateErrorIoni( aTrack );
597 
598  return 0;
599 }
600 
601 
602 //------------------------------------------------------------------------
604 {
605  G4ThreeVector vpPre = aTrack->GetMomentum()/GeV;
606  G4double pPre = vpPre.mag();
607  G4double pBeta = pPre*pPre / (aTrack->GetTotalEnergy()/GeV);
608  G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
609 
610  G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
611  G4double effZ, effA;
612  CalculateEffectiveZandA( mate, effZ, effA );
613 
614 #ifdef G4EVERBOSE
615  if( iverbose >= 4 ) G4cout << "material " << mate->GetName()
616  //<< " " << mate->GetZ() << " " << mate->GetA()
617  << " effZ:" << effZ << " effA:" << effA
618  << " dens(g/mole):" << mate->GetDensity()/g*mole << " Radlen/cm:" << mate->GetRadlen()/cm << " nuclLen/cm" << mate->GetNuclearInterLength()/cm << G4endl;
619 #endif
620 
621  G4double RI = stepLengthCm / (mate->GetRadlen()/cm);
622 #ifdef G4EVERBOSE
623  if( iverbose >= 4 ) G4cout << std::setprecision(6) << std::setw(6) << "G4EP:MSC: RI=X/X0 " << RI << " stepLengthCm " << stepLengthCm << " radlen/cm " << (mate->GetRadlen()/cm) << " RI*1.e10:" << RI*1.e10 << G4endl;
624 #endif
625  G4double charge = aTrack->GetDynamicParticle()->GetCharge();
626  G4double DD = 1.8496E-4*RI*(charge/pBeta * charge/pBeta );
627 #ifdef G4EVERBOSE
628  if( iverbose >= 3 ) G4cout << "G4EP:MSC: D*1E6= " << DD*1.E6 <<" pBeta " << pBeta << G4endl;
629 #endif
630  G4double S1 = DD*stepLengthCm*stepLengthCm/3.;
631  G4double S2 = DD;
632  G4double S3 = DD*stepLengthCm/2.;
633 
634  G4double CLA = std::sqrt( vpPre.x() * vpPre.x() + vpPre.y() * vpPre.y() )/pPre;
635 #ifdef G4EVERBOSE
636  if( iverbose >= 2 ) G4cout << std::setw(6) << "G4EP:MSC: RI " << RI << " S1 " << S1 << " S2 " << S2 << " S3 " << S3 << " CLA " << CLA << G4endl;
637 #endif
638  fError[1][1] += S2;
639  fError[1][4] -= S3;
640  fError[2][2] += S2/CLA/CLA;
641  fError[2][3] += S3/CLA;
642  fError[3][3] += S1;
643  fError[4][4] += S1;
644 
645 #ifdef G4EVERBOSE
646  if( iverbose >= 2 ) G4cout << "G4EP:MSC: error matrix propagated msc " << fError << G4endl;
647 #endif
648 
649  return 0;
650 }
651 
652 
653 //------------------------------------------------------------------------
655 {
656  effZ = 0.;
657  effA = 0.;
658  G4int ii, nelem = mate->GetNumberOfElements();
659  const G4double* fracVec = mate->GetFractionVector();
660  for(ii=0; ii < nelem; ii++ ) {
661  effZ += mate->GetElement( ii )->GetZ() * fracVec[ii];
662  effA += mate->GetElement( ii )->GetA() * fracVec[ii] /g*mole;
663  }
664 
665 }
666 
667 
668 //------------------------------------------------------------------------
670 {
671  G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
672 #ifdef G4EVERBOSE
673  G4double DEDX2;
674  if( stepLengthCm < 1.E-7 ) {
675  DEDX2=0.;
676  }
677 #endif
678  // * Calculate xi factor (KeV).
679  G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
680  G4double effZ, effA;
681  CalculateEffectiveZandA( mate, effZ, effA );
682 
683  G4double Etot = aTrack->GetTotalEnergy()/GeV;
684  G4double beta = aTrack->GetMomentum().mag()/GeV / Etot;
685  G4double mass = aTrack->GetDynamicParticle()->GetMass() / GeV;
686  G4double gamma = Etot / mass;
687 
688  // * Calculate xi factor (keV).
689  G4double XI = 153.5*effZ*stepLengthCm*(mate->GetDensity()/mg*mole) /
690  (effA*beta*beta);
691 
692 #ifdef G4EVERBOSE
693  if( iverbose >= 2 ){
694  G4cout << "G4EP:IONI: XI/keV " << XI << " beta " << beta << " gamma " << gamma << G4endl;
695  G4cout << " density " << (mate->GetDensity()/mg*mole) << " effA " << effA << " step " << stepLengthCm << G4endl;
696  }
697 #endif
698  // * Maximum energy transfer to atomic electron (KeV).
699  G4double eta = beta*gamma;
700  G4double etasq = eta*eta;
701  G4double eMass = 0.51099906/GeV;
702  G4double massRatio = eMass / mass;
703  G4double F1 = 2*eMass*etasq;
704  G4double F2 = 1. + 2. * massRatio * gamma + massRatio * massRatio;
705  G4double Emax = 1.E+6*F1/F2; // now in keV
706 
707  // * *** and now sigma**2 in GeV
708  G4double dedxSq = XI*Emax*(1.-(beta*beta/2.))*1.E-12; // now in GeV^2
709  /*The above formula for var(1/p) good for dens scatterers. However, for MIPS passing
710  through a gas it leads to overestimation. Further more for incident electrons
711  the Emax is almost equal to incident energy.
712  This leads to k=Xi/Emax as small as e-6 and gradually the cov matrix explodes.
713 
714  http://www2.pv.infn.it/~rotondi/kalman_1.pdf
715 
716  Since I do not have enough info at the moment to implement Landau & sub-Landau models for k=Xi/Emax <0.01 I'll saturate k at this value for now
717  */
718 
719  if (XI/Emax<0.01) dedxSq *=XI/Emax*100 ;// Quench for low Elos, see above: newVar=odVar *k/0.01
720 
721 #ifdef G4EVERBOSE
722  if( iverbose >= 2 ) G4cout << "G4EP:IONI: DEDX^2(GeV^2) " << dedxSq << " emass/GeV: " << eMass << " Emax/keV: " << Emax
723  <<" k=Xi/Emax="<< XI/Emax<< G4endl;
724 
725 #endif
726 
727  G4double pPre6 = (aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV).mag();
728  pPre6 = std::pow(pPre6, 6 );
729  //Apply it to error
730  fError[0][0] += Etot*Etot*dedxSq / pPre6;
731 #ifdef G4EVERBOSE
732  if( iverbose >= 2 ) G4cout << "G4:IONI Etot/GeV: " << Etot << " err_dedx^2/GeV^2: " << dedxSq << " p^6: " << pPre6 << G4endl;
733  if( iverbose >= 2 ) G4cout << "G4EP:IONI: error2_from_ionisation " << (Etot*Etot*dedxSq) / pPre6 << G4endl;
734 #endif
735 
736  return 0;
737 }
double perp() const
virtual G4int Update(const G4Track *aTrack)
CLHEP::Hep3Vector G4ThreeVector
const G4double * GetFractionVector() const
Definition: G4Material.hh:195
static const G4double pos
G4LogicalVolume * GetLogicalVolume() const
static constexpr double MeV
Definition: G4SIunits.hh:214
std::ostream & operator<<(std::ostream &, const BasicVector3D< float > &)
G4StepPoint * GetPreStepPoint() const
G4int PropagateErrorMSC(const G4Track *aTrack)
#define G4endl
Definition: G4ios.hh:61
G4double GetRadlen() const
Definition: G4Material.hh:221
G4ErrorMatrix T() const
void Update(const G4Track *aTrack)
BasicVector3D< T > cross(const BasicVector3D< T > &v) const
G4double GetLambda() const
G4ErrorSymMatrix G4ErrorTrajErr
double z() const
virtual void GetFieldValue(const double Point[4], double *fieldArr) const =0
const G4double kCarTolerance
Double_t beta
double theta() const
G4Material * GetMaterial() const
const G4String & GetName() const
Definition: G4Material.hh:179
static constexpr double g
Definition: G4SIunits.hh:183
void CalculateEffectiveZandA(const G4Material *mate, double &effZ, double &effA)
G4ThreeVector GetMomentum() const
double G4double
Definition: G4Types.hh:76
virtual G4int PropagateError(const G4Track *aTrack)
G4FieldManager * GetFieldManager() const
G4ErrorSurfaceTrajParam GetParameters() const
G4double GetStepLength() const
G4ErrorFreeTrajParam fTrajParam
const G4ThreeVector & GetPosition() const
static constexpr double mg
Definition: G4SIunits.hh:184
void DumpPosMomError(std::ostream &out=G4cout) const
G4double GetCharge() const
G4ErrorFreeTrajParam GetParameters() const
const G4ThreeVector & GetPosition() const
const G4Field * GetDetectorField() const
G4double GetPhi() const
G4ErrorSymMatrix similarity(const G4ErrorMatrix &m1) const
const G4Step * GetStep() const
G4ErrorTrajErr fError
HepGeom::Vector3D< G4double > G4Vector3D
Definition: G4Vector3D.hh:35
Hep3Vector cross(const Hep3Vector &) const
static G4TransportationManager * GetTransportationManager()
TH1F * h2
double mag() const
G4int PropagateErrorIoni(const G4Track *aTrack)
int G4int
Definition: G4Types.hh:78
static constexpr double c_light
G4ErrorTrajErr GetError() const
G4double GetA() const
Definition: G4Element.hh:139
TDirectory * dir
void UpdatePosMom(const G4Point3D &pos, const G4Vector3D &mom)
static const G4double Emax
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
double x() const
TH1F * h1
static G4GeometryTolerance * GetInstance()
G4ThreeVector GetMomentum() const
G4Vector3D GetDirection() const
G4double GetZ() const
Definition: G4Element.hh:131
static double Q[]
G4double GetMass() const
static constexpr double mole
Definition: G4SIunits.hh:286
double y() const
virtual void Dump(std::ostream &out=G4cout) const
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:203
G4ErrorSymMatrix T() const
static constexpr double GeV
Definition: G4SIunits.hh:217
const G4DynamicParticle * GetDynamicParticle() const
G4double GetNuclearInterLength() const
Definition: G4Material.hh:224
G4double GetSurfaceTolerance() const
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
static constexpr double tesla
Definition: G4SIunits.hh:268
G4VPhysicalVolume * GetVolume() const
G4double GetDensity() const
Definition: G4Material.hh:181
G4double GetTotalEnergy() const
static G4ErrorPropagatorData * GetErrorPropagatorData()