ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ErrorFreeTrajState.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
27 // ------------------------------------------------------------
28 // GEANT 4 class implementation file
29 // ------------------------------------------------------------
30 //
31 
32 #include <iomanip>
33 
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4Field.hh"
37 #include "G4FieldManager.hh"
38 #include "G4TransportationManager.hh"
39 #include "G4GeometryTolerance.hh"
40 #include "G4Material.hh"
41 #include "G4ErrorPropagatorData.hh"
42 #include "G4ErrorFreeTrajState.hh"
43 #include "G4ErrorFreeTrajParam.hh"
44 #include "G4ErrorSurfaceTrajState.hh"
45 #include "G4ErrorMatrix.hh"
46 
47 //------------------------------------------------------------------------
48 G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4String& partName, const G4Point3D& pos, const G4Vector3D& mom, const G4ErrorTrajErr& errmat) : G4ErrorTrajState( partName, pos, mom, errmat )
49 {
50  fTrajParam = G4ErrorFreeTrajParam( pos, mom );
51  Init();
52 }
53 
54 
55 //------------------------------------------------------------------------
56 G4ErrorFreeTrajState::G4ErrorFreeTrajState( const G4ErrorSurfaceTrajState& tpSD ) : G4ErrorTrajState( tpSD.GetParticleType(), tpSD.GetPosition(), tpSD.GetMomentum() )
57 {
58  // G4ThreeVector planeNormal = tpSD.GetPlaneNormal();
59  // G4double fPt = tpSD.GetMomentum()*planeNormal;//mom projected on normal to plane
60  // G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
61  // G4ThreeVector Psc = fPt * planeNormal + tpSDparam.GetPU()*tpSDparam.GetVectorU() + tpSD.GetPV()*tpSD.GetVectorW();
62 
63  fTrajParam = G4ErrorFreeTrajParam( fPosition, fMomentum );
64  Init();
65 
66  //----- Get the error matrix in SC coordinates
67  G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
68  G4double mom = fMomentum.mag();
69  G4double mom2 = fMomentum.mag2();
70  G4double TVW1 = std::sqrt( mom2 / ( mom2 + tpSDparam.GetPV()*tpSDparam.GetPV() + tpSDparam.GetPW()*tpSDparam.GetPW()) );
71  G4ThreeVector vTVW( TVW1, tpSDparam.GetPV()/mom * TVW1, tpSDparam.GetPW()/mom * TVW1 );
72  G4Vector3D vectorU = tpSDparam.GetVectorV().cross( tpSDparam.GetVectorW() );
73  G4Vector3D vTN = vTVW.x()*vectorU + vTVW.y()*tpSDparam.GetVectorV() + vTVW.z()*tpSDparam.GetVectorW();
74 
75 #ifdef G4EVERBOSE
76  if( iverbose >= 5){
77  G4double pc2 = std::asin( vTN.z() );
78  G4double pc3 = std::atan (vTN.y()/vTN.x());
79 
80  G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda() << " diff " << pc2-GetParameters().GetLambda() << G4endl;
81  G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi() << " diff " << pc3-GetParameters().GetPhi() << G4endl;
82  }
83 #endif
84 
85  //--- Get the unit vectors perp to P
86  G4double cosl = std::cos( GetParameters().GetLambda() );
87  if (cosl < 1.E-30) cosl = 1.E-30;
88  G4double cosl1 = 1./cosl;
89  G4Vector3D vUN(-vTN.y()*cosl1, vTN.x()*cosl1, 0. );
90  G4Vector3D vVN(-vTN.z()*vUN.y(), vTN.z()*vUN.x(), cosl );
91 
92  G4Vector3D vUperp = G4Vector3D( -fMomentum.y(), fMomentum.x(), 0.);
93  G4Vector3D vVperp = vUperp.cross( fMomentum );
94  vUperp *= 1./vUperp.mag();
95  vVperp *= 1./vVperp.mag();
96 
97 #ifdef G4EVERBOSE
98  if( iverbose >= 5 ){
99  G4cout << " CHECK: vUN " << vUN << " = " << vUperp << " diff " << (vUN-vUperp).mag() << G4endl;
100  G4cout << " CHECK: vVN " << vVN << " = " << vVperp << " diff " << (vVN-vVperp).mag() << G4endl;
101  }
102 #endif
103 
104  //get the dot products of vectors perpendicular to direction and vector defining SD plane
105  G4double dUU = vUperp * tpSD.GetVectorV();
106  G4double dUV = vUperp * tpSD.GetVectorW();
107  G4double dVU = vVperp * tpSD.GetVectorV();
108  G4double dVV = vVperp * tpSD.GetVectorW();
109 
110  //--- Get transformation first
111  G4ErrorMatrix transfM(5, 5, 1 );
112  //--- Get magnetic field
113  const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
114  G4ThreeVector dir = fTrajParam.GetDirection();
115  G4double invCosTheta = 1./std::cos( dir.theta() );
116  G4cout << " dir="<<dir<<" invCosTheta "<<invCosTheta << G4endl;
117 
118  if( fCharge != 0 && field ) {
119  G4double pos1[3]; pos1[0] = fPosition.x()*cm; pos1[1] = fPosition.y()*cm; pos1[2] = fPosition.z()*cm;
120  G4double h1[3];
121  field->GetFieldValue( pos1, h1 );
122  G4ThreeVector HPre = G4ThreeVector( h1[0], h1[1], h1[2] ) / tesla *10.;
123  G4double magHPre = HPre.mag();
124  G4double invP = 1./fMomentum.mag();
125  G4double magHPreM = magHPre * invP;
126  if( magHPre != 0. ) {
127  G4double magHPreM2 = fCharge / magHPre;
128 
129  G4double Q = -magHPreM * c_light;
130  G4double sinz = -HPre*vUperp * magHPreM2;
131  G4double cosz = HPre*vVperp * magHPreM2;
132 
133  transfM[1][3] = -Q*dir.y()*sinz;
134  transfM[1][4] = -Q*dir.z()*sinz;
135  transfM[2][3] = -Q*dir.y()*cosz*invCosTheta;
136  transfM[2][4] = -Q*dir.z()*cosz*invCosTheta;
137  }
138  }
139 
140  transfM[0][0] = 1.;
141  transfM[1][1] = dir.x()*dVU;
142  transfM[1][2] = dir.x()*dVV;
143  transfM[2][1] = dir.x()*dUU*invCosTheta;
144  transfM[2][2] = dir.x()*dUV*invCosTheta;
145  transfM[3][3] = dUU;
146  transfM[3][4] = dUV;
147  transfM[4][3] = dVU;
148  transfM[4][4] = dVV;
149 
150  fError = G4ErrorTrajErr( tpSD.GetError().similarity( transfM ) );
151 
152 #ifdef G4EVERBOSE
153  if( iverbose >= 1) G4cout << "error matrix SD2SC " << fError << G4endl;
154  if( iverbose >= 4) G4cout << "G4ErrorFreeTrajState from SD " << *this << G4endl;
155 #endif
156 }
157 
158 
159 //------------------------------------------------------------------------
160 void G4ErrorFreeTrajState::Init()
161 {
162  theTSType = G4eTS_FREE;
163  BuildCharge();
164  theTransfMat = G4ErrorMatrix(5,5,0);
165  theFirstStep = true;
166 }
167 
168 //------------------------------------------------------------------------
169 void G4ErrorFreeTrajState::Dump( std::ostream& out ) const
170 {
171  out << *this;
172 }
173 
174 //------------------------------------------------------------------------
175 G4int G4ErrorFreeTrajState::Update( const G4Track* aTrack )
176 {
177  G4int ierr = 0;
178  fTrajParam.Update( aTrack );
179  UpdatePosMom( aTrack->GetPosition(), aTrack->GetMomentum() );
180  return ierr;
181 }
182 
183 
184 //------------------------------------------------------------------------
185 std::ostream& operator<<(std::ostream& out, const G4ErrorFreeTrajState& ts)
186 {
187  std::ios::fmtflags orig_flags = out.flags();
188 
189  out.setf(std::ios::fixed,std::ios::floatfield);
190 
191  ts.DumpPosMomError( out );
192 
193  out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl;
194 
195  out.flags(orig_flags);
196 
197  return out;
198 }
199 
200 
201 //------------------------------------------------------------------------
202 G4int G4ErrorFreeTrajState::PropagateError( const G4Track* aTrack )
203 {
204  G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
205  if( G4ErrorPropagatorData::GetErrorPropagatorData()->GetStage() == G4ErrorStage_Deflation ) stepLengthCm *= -1.;
206 
207  G4double kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
208 
209  if( std::fabs(stepLengthCm) <= kCarTolerance/cm ) return 0;
210 
211 #ifdef G4EVERBOSE
212  if( iverbose >= 2 )G4cout << " G4ErrorFreeTrajState::PropagateError " << G4endl;
213  G4cout << "G4EP: iverbose="<< iverbose << G4endl;
214 #endif
215 
216  // * *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
217  G4Point3D vposPost = aTrack->GetPosition()/cm;
218  G4Vector3D vpPost = aTrack->GetMomentum()/GeV;
219  // G4Point3D vposPre = fPosition/cm;
220  // G4Vector3D vpPre = fMomentum/GeV;
221  G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition()/cm;
222  G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV;
223  //correct to avoid propagation along Z
224  if( vpPre.mag() == vpPre.z() ) vpPre.setX( 1.E-6*MeV );
225  if( vpPost.mag() == vpPost.z() ) vpPost.setX( 1.E-6*MeV );
226 
227  G4double pPre = vpPre.mag();
228  G4double pPost = vpPost.mag();
229 #ifdef G4EVERBOSE
230  if( iverbose >= 2 ) {
231  G4cout << "G4EP: vposPre " << vposPre << G4endl
232  << "G4EP: vposPost " << vposPost << G4endl;
233  G4cout << "G4EP: vpPre " << vpPre << G4endl
234  << "G4EP: vpPost " << vpPost << G4endl;
235  G4cout << " err start step " << fError << G4endl;
236  G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl;
237  }
238 #endif
239 
240  if( pPre == 0. || pPost == 0 ) return 2;
241  G4double pInvPre = 1./pPre;
242  G4double pInvPost = 1./pPost;
243  G4double deltaPInv = pInvPost - pInvPre;
244  if( iverbose >= 2 ) G4cout << "G4EP: pInvPre" << pInvPre<< " pInvPost:" << pInvPost<<" deltaPInv:" << deltaPInv<< G4endl;
245 
246 
247  G4Vector3D vpPreNorm = vpPre * pInvPre;
248  G4Vector3D vpPostNorm = vpPost * pInvPost;
249  if( iverbose >= 2 ) G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm << G4endl;
250  //return if propagation along Z??
251  if( 1. - std::fabs(vpPreNorm.z()) < kCarTolerance ) return 4;
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( vpPostNorm.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 }
599 
600 
601 //------------------------------------------------------------------------
602 G4int G4ErrorFreeTrajState::PropagateErrorMSC( const G4Track* aTrack )
603 {
604  G4ThreeVector vpPre = aTrack->GetMomentum()/GeV;
605  G4double pPre = vpPre.mag();
606  G4double pBeta = pPre*pPre / (aTrack->GetTotalEnergy()/GeV);
607  G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
608 
609  G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
610  G4double effZ, effA;
611  CalculateEffectiveZandA( mate, effZ, effA );
612 
613 #ifdef G4EVERBOSE
614  if( iverbose >= 4 ) G4cout << "material " << mate->GetName()
615  //<< " " << mate->GetZ() << " " << mate->GetA()
616  << " effZ:" << effZ << " effA:" << effA
617  << " dens(g/mole):" << mate->GetDensity()/g*mole << " Radlen/cm:" << mate->GetRadlen()/cm << " nuclLen/cm" << mate->GetNuclearInterLength()/cm << G4endl;
618 #endif
619 
620  G4double RI = stepLengthCm / (mate->GetRadlen()/cm);
621 #ifdef G4EVERBOSE
622  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;
623 #endif
624  G4double charge = aTrack->GetDynamicParticle()->GetCharge();
625  G4double DD = 1.8496E-4*RI*(charge/pBeta * charge/pBeta );
626 #ifdef G4EVERBOSE
627  if( iverbose >= 3 ) G4cout << "G4EP:MSC: D*1E6= " << DD*1.E6 <<" pBeta " << pBeta << G4endl;
628 #endif
629  G4double S1 = DD*stepLengthCm*stepLengthCm/3.;
630  G4double S2 = DD;
631  G4double S3 = DD*stepLengthCm/2.;
632 
633  G4double CLA = std::sqrt( vpPre.x() * vpPre.x() + vpPre.y() * vpPre.y() )/pPre;
634 #ifdef G4EVERBOSE
635  if( iverbose >= 2 ) G4cout << std::setw(6) << "G4EP:MSC: RI " << RI << " S1 " << S1 << " S2 " << S2 << " S3 " << S3 << " CLA " << CLA << G4endl;
636 #endif
637  fError[1][1] += S2;
638  fError[1][4] -= S3;
639  fError[2][2] += S2/CLA/CLA;
640  fError[2][3] += S3/CLA;
641  fError[3][3] += S1;
642  fError[4][4] += S1;
643 
644 #ifdef G4EVERBOSE
645  if( iverbose >= 2 ) G4cout << "G4EP:MSC: error matrix propagated msc " << fError << G4endl;
646 #endif
647 
648  return 0;
649 }
650 
651 
652 //------------------------------------------------------------------------
653 void G4ErrorFreeTrajState::CalculateEffectiveZandA( const G4Material* mate, G4double& effZ, G4double& effA )
654 {
655  effZ = 0.;
656  effA = 0.;
657  G4int ii, nelem = mate->GetNumberOfElements();
658  const G4double* fracVec = mate->GetFractionVector();
659  for(ii=0; ii < nelem; ii++ ) {
660  effZ += mate->GetElement( ii )->GetZ() * fracVec[ii];
661  effA += mate->GetElement( ii )->GetA() * fracVec[ii] /g*mole;
662  }
663 
664 }
665 
666 
667 //------------------------------------------------------------------------
668 G4int G4ErrorFreeTrajState::PropagateErrorIoni( const G4Track* aTrack )
669 {
670  G4double stepLengthCm = aTrack->GetStep()->GetStepLength()/cm;
671 #ifdef G4EVERBOSE
672  G4double DEDX2;
673  if( stepLengthCm < 1.E-7 ) {
674  DEDX2=0.;
675  }
676 #endif
677  // * Calculate xi factor (KeV).
678  G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
679  G4double effZ, effA;
680  CalculateEffectiveZandA( mate, effZ, effA );
681 
682  G4double Etot = aTrack->GetTotalEnergy()/GeV;
683  G4double beta = aTrack->GetMomentum().mag()/GeV / Etot;
684  G4double mass = aTrack->GetDynamicParticle()->GetMass() / GeV;
685  G4double gamma = Etot / mass;
686 
687  // * Calculate xi factor (keV).
688  G4double XI = 153.5*effZ*stepLengthCm*(mate->GetDensity()/mg*mole) /
689  (effA*beta*beta);
690 
691 #ifdef G4EVERBOSE
692  if( iverbose >= 2 ){
693  G4cout << "G4EP:IONI: XI/keV " << XI << " beta " << beta << " gamma " << gamma << G4endl;
694  G4cout << " density " << (mate->GetDensity()/mg*mole) << " effA " << effA << " step " << stepLengthCm << G4endl;
695  }
696 #endif
697  // * Maximum energy transfer to atomic electron (KeV).
698  G4double eta = beta*gamma;
699  G4double etasq = eta*eta;
700  G4double eMass = 0.51099906/GeV;
701  G4double massRatio = eMass / mass;
702  G4double F1 = 2*eMass*etasq;
703  G4double F2 = 1. + 2. * massRatio * gamma + massRatio * massRatio;
704  G4double Emax = 1.E+6*F1/F2; // now in keV
705 
706  // * *** and now sigma**2 in GeV
707  G4double dedxSq = XI*Emax*(1.-(beta*beta/2.))*1.E-12; // now in GeV^2
708  /*The above formula for var(1/p) good for dens scatterers. However, for MIPS passing
709  through a gas it leads to overestimation. Further more for incident electrons
710  the Emax is almost equal to incident energy.
711  This leads to k=Xi/Emax as small as e-6 and gradually the cov matrix explodes.
712 
713  http://www2.pv.infn.it/~rotondi/kalman_1.pdf
714 
715  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
716  */
717 
718  if (XI/Emax<0.01) dedxSq *=XI/Emax*100 ;// Quench for low Elos, see above: newVar=odVar *k/0.01
719 
720 #ifdef G4EVERBOSE
721  if( iverbose >= 2 ) G4cout << "G4EP:IONI: DEDX^2(GeV^2) " << dedxSq << " emass/GeV: " << eMass << " Emax/keV: " << Emax
722  <<" k=Xi/Emax="<< XI/Emax<< G4endl;
723 
724 #endif
725 
726  G4double pPre6 = (aTrack->GetStep()->GetPreStepPoint()->GetMomentum()/GeV).mag();
727  pPre6 = std::pow(pPre6, 6 );
728  //Apply it to error
729  fError[0][0] += Etot*Etot*dedxSq / pPre6;
730 #ifdef G4EVERBOSE
731  if( iverbose >= 2 ) G4cout << "G4:IONI Etot/GeV: " << Etot << " err_dedx^2/GeV^2: " << dedxSq << " p^6: " << pPre6 << G4endl;
732  if( iverbose >= 2 ) G4cout << "G4EP:IONI: error2_from_ionisation " << (Etot*Etot*dedxSq) / pPre6 << G4endl;
733 #endif
734 
735  return 0;
736 }