| Geant4 Cross Reference |
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,
49 const G4Point3D& pos,
50 const G4Vector3D& mom,
51 const G4ErrorTrajErr& errmat)
52 : G4ErrorTrajState(partName, pos, mom, errmat)
53 {
54 fTrajParam = G4ErrorFreeTrajParam(pos, mom);
55 Init();
56 }
57
58 //------------------------------------------------------------------------
59 G4ErrorFreeTrajState::G4ErrorFreeTrajState(const G4ErrorSurfaceTrajState& tpSD)
60 : G4ErrorTrajState(tpSD.GetParticleType(), tpSD.GetPosition(),
61 tpSD.GetMomentum())
62 {
63 // G4ThreeVector planeNormal = tpSD.GetPlaneNormal();
64 // G4double fPt = tpSD.GetMomentum()*planeNormal;//mom projected on normal to
65 // plane G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
66 // G4ThreeVector Psc = fPt * planeNormal +
67 // tpSDparam.GetPU()*tpSDparam.GetVectorU() + tpSD.GetPV()*tpSD.GetVectorW();
68
69 fTrajParam = G4ErrorFreeTrajParam(fPosition, fMomentum);
70 Init();
71
72 //----- Get the error matrix in SC coordinates
73 G4ErrorSurfaceTrajParam tpSDparam = tpSD.GetParameters();
74 G4double mom = fMomentum.mag();
75 G4double mom2 = fMomentum.mag2();
76 G4double TVW1 =
77 std::sqrt(mom2 / (mom2 + tpSDparam.GetPV() * tpSDparam.GetPV() +
78 tpSDparam.GetPW() * tpSDparam.GetPW()));
79 G4ThreeVector vTVW(TVW1, tpSDparam.GetPV() / mom * TVW1,
80 tpSDparam.GetPW() / mom * TVW1);
81 G4Vector3D vectorU = tpSDparam.GetVectorV().cross(tpSDparam.GetVectorW());
82 G4Vector3D vTN = vTVW.x() * vectorU + vTVW.y() * tpSDparam.GetVectorV() +
83 vTVW.z() * tpSDparam.GetVectorW();
84
85 #ifdef G4EVERBOSE
86 if(iverbose >= 5)
87 {
88 G4double pc2 = std::asin(vTN.z());
89 G4double pc3 = std::atan(vTN.y() / vTN.x());
90
91 G4cout << " CHECK: pc2 " << pc2 << " = " << GetParameters().GetLambda()
92 << " diff " << pc2 - GetParameters().GetLambda() << G4endl;
93 G4cout << " CHECK: pc3 " << pc3 << " = " << GetParameters().GetPhi()
94 << " diff " << pc3 - GetParameters().GetPhi() << G4endl;
95 }
96 #endif
97
98 //--- Get the unit vectors perp to P
99 G4double cosl = std::cos(GetParameters().GetLambda());
100 if(cosl < 1.E-30)
101 cosl = 1.E-30;
102 G4double cosl1 = 1. / cosl;
103 G4Vector3D vUN(-vTN.y() * cosl1, vTN.x() * cosl1, 0.);
104 G4Vector3D vVN(-vTN.z() * vUN.y(), vTN.z() * vUN.x(), cosl);
105
106 G4Vector3D vUperp = G4Vector3D(-fMomentum.y(), fMomentum.x(), 0.);
107 G4Vector3D vVperp = vUperp.cross(fMomentum);
108 vUperp *= 1. / vUperp.mag();
109 vVperp *= 1. / vVperp.mag();
110
111 #ifdef G4EVERBOSE
112 if(iverbose >= 5)
113 {
114 G4cout << " CHECK: vUN " << vUN << " = " << vUperp << " diff "
115 << (vUN - vUperp).mag() << G4endl;
116 G4cout << " CHECK: vVN " << vVN << " = " << vVperp << " diff "
117 << (vVN - vVperp).mag() << G4endl;
118 }
119 #endif
120
121 // get the dot products of vectors perpendicular to direction and vector
122 // defining SD plane
123 G4double dUU = vUperp * tpSD.GetVectorV();
124 G4double dUV = vUperp * tpSD.GetVectorW();
125 G4double dVU = vVperp * tpSD.GetVectorV();
126 G4double dVV = vVperp * tpSD.GetVectorW();
127
128 //--- Get transformation first
129 G4ErrorMatrix transfM(5, 5, 1);
130 //--- Get magnetic field
131 const G4Field* field = G4TransportationManager::GetTransportationManager()
132 ->GetFieldManager()
133 ->GetDetectorField();
134 G4ThreeVector dir = fTrajParam.GetDirection();
135 G4double invCosTheta = 1. / std::cos(dir.theta());
136 G4cout << " dir=" << dir << " invCosTheta " << invCosTheta << G4endl;
137
138 if(fCharge != 0 && field)
139 {
140 G4double pos1[3];
141 pos1[0] = fPosition.x() * cm;
142 pos1[1] = fPosition.y() * cm;
143 pos1[2] = fPosition.z() * cm;
144 G4double h1[3];
145 field->GetFieldValue(pos1, h1);
146 G4ThreeVector HPre = G4ThreeVector(h1[0], h1[1], h1[2]) / tesla * 10.;
147 G4double magHPre = HPre.mag();
148 G4double invP = 1. / fMomentum.mag();
149 G4double magHPreM = magHPre * invP;
150 if(magHPre != 0.)
151 {
152 G4double magHPreM2 = fCharge / magHPre;
153
154 G4double Q = -magHPreM * c_light;
155 G4double sinz = -HPre * vUperp * magHPreM2;
156 G4double cosz = HPre * vVperp * magHPreM2;
157
158 transfM[1][3] = -Q * dir.y() * sinz;
159 transfM[1][4] = -Q * dir.z() * sinz;
160 transfM[2][3] = -Q * dir.y() * cosz * invCosTheta;
161 transfM[2][4] = -Q * dir.z() * cosz * invCosTheta;
162 }
163 }
164
165 transfM[0][0] = 1.;
166 transfM[1][1] = dir.x() * dVU;
167 transfM[1][2] = dir.x() * dVV;
168 transfM[2][1] = dir.x() * dUU * invCosTheta;
169 transfM[2][2] = dir.x() * dUV * invCosTheta;
170 transfM[3][3] = dUU;
171 transfM[3][4] = dUV;
172 transfM[4][3] = dVU;
173 transfM[4][4] = dVV;
174
175 fError = G4ErrorTrajErr(tpSD.GetError().similarity(transfM));
176
177 #ifdef G4EVERBOSE
178 if(iverbose >= 1)
179 G4cout << "error matrix SD2SC " << fError << G4endl;
180 if(iverbose >= 4)
181 G4cout << "G4ErrorFreeTrajState from SD " << *this << G4endl;
182 #endif
183 }
184
185 //------------------------------------------------------------------------
186 void G4ErrorFreeTrajState::Init()
187 {
188 theTSType = G4eTS_FREE;
189 BuildCharge();
190 theTransfMat = G4ErrorMatrix(5, 5, 0);
191 theFirstStep = true;
192 }
193
194 //------------------------------------------------------------------------
195 void G4ErrorFreeTrajState::Dump(std::ostream& out) const { out << *this; }
196
197 //------------------------------------------------------------------------
198 G4int G4ErrorFreeTrajState::Update(const G4Track* aTrack)
199 {
200 G4int ierr = 0;
201 fTrajParam.Update(aTrack);
202 UpdatePosMom(aTrack->GetPosition(), aTrack->GetMomentum());
203 return ierr;
204 }
205
206 //------------------------------------------------------------------------
207 std::ostream& operator<<(std::ostream& out, const G4ErrorFreeTrajState& ts)
208 {
209 std::ios::fmtflags orig_flags = out.flags();
210
211 out.setf(std::ios::fixed, std::ios::floatfield);
212
213 ts.DumpPosMomError(out);
214
215 out << " G4ErrorFreeTrajState: Params: " << ts.fTrajParam << G4endl;
216
217 out.flags(orig_flags);
218
219 return out;
220 }
221
222 //------------------------------------------------------------------------
223 G4int G4ErrorFreeTrajState::PropagateError(const G4Track* aTrack)
224 {
225 G4double stepLengthCm = aTrack->GetStep()->GetStepLength() / cm;
226 if(G4ErrorPropagatorData::GetErrorPropagatorData()->GetStage() ==
227 G4ErrorStage_Deflation)
228 stepLengthCm *= -1.;
229
230 G4double kCarTolerance =
231 G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
232
233 if(std::fabs(stepLengthCm) <= kCarTolerance / cm)
234 return 0;
235
236 #ifdef G4EVERBOSE
237 if(iverbose >= 2)
238 G4cout << " G4ErrorFreeTrajState::PropagateError " << G4endl;
239 G4cout << "G4EP: iverbose=" << iverbose << G4endl;
240 #endif
241
242 // * *** ERROR PROPAGATION ON A HELIX ASSUMING SC VARIABLES
243 G4Point3D vposPost = aTrack->GetPosition() / cm;
244 G4Vector3D vpPost = aTrack->GetMomentum() / GeV;
245 // G4Point3D vposPre = fPosition/cm;
246 // G4Vector3D vpPre = fMomentum/GeV;
247 G4Point3D vposPre = aTrack->GetStep()->GetPreStepPoint()->GetPosition() / cm;
248 G4Vector3D vpPre = aTrack->GetStep()->GetPreStepPoint()->GetMomentum() / GeV;
249 // correct to avoid propagation along Z
250 if(vpPre.mag() == vpPre.z())
251 vpPre.setX(1.E-6 * MeV);
252 if(vpPost.mag() == vpPost.z())
253 vpPost.setX(1.E-6 * MeV);
254
255 G4double pPre = vpPre.mag();
256 G4double pPost = vpPost.mag();
257 #ifdef G4EVERBOSE
258 if(iverbose >= 2)
259 {
260 G4cout << "G4EP: vposPre " << vposPre << G4endl << "G4EP: vposPost "
261 << vposPost << G4endl;
262 G4cout << "G4EP: vpPre " << vpPre << G4endl << "G4EP: vpPost " << vpPost
263 << G4endl;
264 G4cout << " err start step " << fError << G4endl;
265 G4cout << "G4EP: stepLengthCm " << stepLengthCm << G4endl;
266 }
267 #endif
268
269 if(pPre == 0. || pPost == 0)
270 return 2;
271 G4double pInvPre = 1. / pPre;
272 G4double pInvPost = 1. / pPost;
273 G4double deltaPInv = pInvPost - pInvPre;
274 if(iverbose >= 2)
275 G4cout << "G4EP: pInvPre" << pInvPre << " pInvPost:" << pInvPost
276 << " deltaPInv:" << deltaPInv << G4endl;
277
278 G4Vector3D vpPreNorm = vpPre * pInvPre;
279 G4Vector3D vpPostNorm = vpPost * pInvPost;
280 if(iverbose >= 2)
281 G4cout << "G4EP: vpPreNorm " << vpPreNorm << " vpPostNorm " << vpPostNorm
282 << G4endl;
283 // return if propagation along Z??
284 if(1. - std::fabs(vpPreNorm.z()) < kCarTolerance)
285 return 4;
286 if(1. - std::fabs(vpPostNorm.z()) < kCarTolerance)
287 return 4;
288 G4double sinpPre =
289 std::sin(vpPreNorm.theta()); // cosine perpendicular to pPre = sine pPre
290 G4double sinpPost =
291 std::sin(vpPostNorm.theta()); // cosine perpendicular to pPost = sine pPost
292 G4double sinpPostInv = 1. / std::sin(vpPostNorm.theta());
293
294 #ifdef G4EVERBOSE
295 if(iverbose >= 2)
296 G4cout << "G4EP: cosl " << sinpPre << " cosl0 " << sinpPost << G4endl;
297 #endif
298 //* *** DEFINE TRANSFORMATION MATRIX BETWEEN X1 AND X2 FOR
299 //* *** NEUTRAL PARTICLE OR FIELDFREE REGION
300 G4ErrorMatrix transf(5, 5, 0);
301
302 transf[3][2] = stepLengthCm * sinpPost;
303 transf[4][1] = stepLengthCm;
304 for(auto ii = 0; ii < 5; ++ii)
305 {
306 transf[ii][ii] = 1.;
307 }
308 #ifdef G4EVERBOSE
309 if(iverbose >= 2)
310 {
311 G4cout << "G4EP: transf matrix neutral " << transf;
312 }
313 #endif
314
315 // charge X propagation direction
316 G4double charge = aTrack->GetDynamicParticle()->GetCharge();
317 if(G4ErrorPropagatorData::GetErrorPropagatorData()->GetMode() ==
318 G4ErrorMode_PropBackwards)
319 {
320 charge *= -1.;
321 }
322 // G4cout << " charge " << charge << G4endl;
323 // t check if particle has charge
324 // t if( charge == 0 ) goto 45;
325 // check if the magnetic field is = 0.
326
327 // position is from geant4, it is assumed to be in mm (for debugging,
328 // eventually it will not be transformed) it is assumed vposPre[] is in cm and
329 // pos1[] is in mm.
330 G4double pos1[3];
331 pos1[0] = vposPre.x() * cm;
332 pos1[1] = vposPre.y() * cm;
333 pos1[2] = vposPre.z() * cm;
334 G4double pos2[3];
335 pos2[0] = vposPost.x() * cm;
336 pos2[1] = vposPost.y() * cm;
337 pos2[2] = vposPost.z() * cm;
338 G4double h1[3], h2[3];
339
340 const G4Field* field = G4TransportationManager::GetTransportationManager()
341 ->GetFieldManager()
342 ->GetDetectorField();
343 if(!field)
344 return 0; // goto 45
345
346 // calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE REGION
347 if(charge != 0. && field)
348 {
349 field->GetFieldValue(pos1,
350 h1); // here pos1[], pos2[] are in mm, not changed
351 field->GetFieldValue(pos2, h2);
352 G4ThreeVector HPre =
353 G4ThreeVector(h1[0], h1[1], h1[2]) / tesla *
354 10.; // 10. is to get same dimensions as GEANT3 (kilogauss)
355 G4ThreeVector HPost = G4ThreeVector(h2[0], h2[1], h2[2]) / tesla * 10.;
356 G4double magHPre = HPre.mag();
357 G4double magHPost = HPost.mag();
358 #ifdef G4EVERBOSE
359 if(iverbose >= 2)
360 {
361 G4cout << "G4EP: h1 = " << h1[0] << ", " << h1[1] << ", " << h1[2]
362 << G4endl;
363 G4cout << "G4EP: pos1/mm = " << pos1[0] << ", " << pos1[1] << ", "
364 << pos1[2] << G4endl;
365 G4cout << "G4EP: pos2/mm = " << pos2[0] << ", " << pos2[1] << ", "
366 << pos2[2] << G4endl;
367 G4cout << "G4EP: B-filed in KGauss HPre " << HPre << G4endl
368 << "G4EP: in KGauss HPost " << HPost << G4endl;
369 }
370 #endif
371
372 if(magHPre + magHPost != 0.)
373 {
374 //* *** CHECK WHETHER H*ALFA/P IS TOO DIFFERENT AT X1 AND X2
375 G4double gam;
376 if(magHPost != 0.)
377 {
378 gam = HPost * vpPostNorm / magHPost;
379 }
380 else
381 {
382 gam = HPre * vpPreNorm / magHPre;
383 }
384
385 // G4eMagneticLimitsProcess will limit the step, but based on an straight
386 // line trajectory
387 G4double alphaSqr = 1. - gam * gam;
388 G4double diffHSqr = (HPre * pInvPre - HPost * pInvPost).mag2();
389 G4double delhp6Sqr = 300. * 300.;
390 #ifdef G4EVERBOSE
391 if(iverbose >= 2)
392 {
393 G4cout << " G4EP: gam " << gam << " alphaSqr " << alphaSqr
394 << " diffHSqr " << diffHSqr << G4endl;
395 G4cout << " alpha= " << std::sqrt(alphaSqr) << G4endl;
396 }
397 #endif
398 if(diffHSqr * alphaSqr > delhp6Sqr)
399 return 3;
400
401 //* *** DEFINE AVERAGE MAGNETIC FIELD AND GRADIENT
402 G4double pInvAver = 1. / (pInvPre + pInvPost);
403 G4double CFACT8 = 2.997925E-4;
404 // G4double HAver
405 G4ThreeVector vHAverNorm((HPre * pInvPre + HPost * pInvPost) * pInvAver *
406 charge * CFACT8);
407 G4double HAver = vHAverNorm.mag();
408 G4double invHAver = 1. / HAver;
409 vHAverNorm *= invHAver;
410 #ifdef G4EVERBOSE
411 if(iverbose >= 2)
412 G4cout << " G4EP: HaverNorm " << vHAverNorm << " magHAver " << HAver
413 << " charge " << charge << G4endl;
414 #endif
415
416 G4double pAver = (pPre + pPost) * 0.5;
417 G4double QAver = -HAver / pAver;
418 G4double thetaAver = QAver * stepLengthCm;
419 G4double sinThetaAver = std::sin(thetaAver);
420 G4double cosThetaAver = std::cos(thetaAver);
421 G4double gamma = vHAverNorm * vpPostNorm;
422 G4ThreeVector AN2 = vHAverNorm.cross(vpPostNorm);
423
424 #ifdef G4EVERBOSE
425 if(iverbose >= 2)
426 G4cout << " G4EP: AN2 " << AN2 << " gamma:" << gamma
427 << " theta=" << thetaAver << G4endl;
428 #endif
429 G4double AU = 1. / vpPreNorm.perp();
430 // t G4ThreeVector vU( vpPreNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU );
431 G4ThreeVector vUPre(-AU * vpPreNorm.y(), AU * vpPreNorm.x(), 0.);
432 G4ThreeVector vVPre(-vpPreNorm.z() * vUPre.y(), vpPreNorm.z() * vUPre.x(),
433 vpPreNorm.x() * vUPre.y() -
434 vpPreNorm.y() * vUPre.x());
435
436 //
437 AU = 1. / vpPostNorm.perp();
438 // t G4ThreeVector vU( vpPostNorm.cross( G4ThreeVector(0.,0.,1.) ) * AU
439 // );
440 G4ThreeVector vUPost(-AU * vpPostNorm.y(), AU * vpPostNorm.x(), 0.);
441 G4ThreeVector vVPost(
442 -vpPostNorm.z() * vUPost.y(), vpPostNorm.z() * vUPost.x(),
443 vpPostNorm.x() * vUPost.y() - vpPostNorm.y() * vUPost.x());
444 #ifdef G4EVERBOSE
445 G4cout << " vpPostNorm " << vpPostNorm << G4endl;
446 if(iverbose >= 2)
447 G4cout << " G4EP: AU " << AU << " vUPre " << vUPre << " vVPre " << vVPre
448 << " vUPost " << vUPost << " vVPost " << vVPost << G4endl;
449 #endif
450 G4Point3D deltaPos(vposPre - vposPost);
451
452 // * *** COMPLETE TRANSFORMATION MATRIX BETWEEN ERRORS AT X1 AND X2
453 // * *** FIELD GRADIENT PERPENDICULAR TO TRACK IS PRESENTLY NOT
454 // * *** TAKEN INTO ACCOUNT
455
456 G4double QP = QAver * pAver; // = -HAver
457 #ifdef G4EVERBOSE
458 if(iverbose >= 2)
459 G4cout << " G4EP: QP " << QP << " QAver " << QAver << " pAver " << pAver
460 << G4endl;
461 #endif
462 G4double ANV =
463 -(vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y());
464 G4double ANU =
465 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
466 vHAverNorm.z() * vVPost.z());
467 G4double OMcosThetaAver = 1. - cosThetaAver;
468 #ifdef G4EVERBOSE
469 if(iverbose >= 2)
470 G4cout << "G4EP: OMcosThetaAver " << OMcosThetaAver << " cosThetaAver "
471 << cosThetaAver << " thetaAver " << thetaAver << " QAver "
472 << QAver << " stepLengthCm " << stepLengthCm << G4endl;
473 #endif
474 G4double TMSINT = thetaAver - sinThetaAver;
475 #ifdef G4EVERBOSE
476 if(iverbose >= 2)
477 G4cout << " G4EP: ANV " << ANV << " ANU " << ANU << G4endl;
478 #endif
479
480 G4ThreeVector vHUPre(
481 -vHAverNorm.z() * vUPre.y(), vHAverNorm.z() * vUPre.x(),
482 vHAverNorm.x() * vUPre.y() - vHAverNorm.y() * vUPre.x());
483 #ifdef G4EVERBOSE
484 // if( iverbose >= 2 ) G4cout << "G4EP: HUPre(1) " << vHUPre.x() << " "
485 // << vHAverNorm.z() << " " << vUPre.y() << G4endl;
486 #endif
487 G4ThreeVector vHVPre(
488 vHAverNorm.y() * vVPre.z() - vHAverNorm.z() * vVPre.y(),
489 vHAverNorm.z() * vVPre.x() - vHAverNorm.x() * vVPre.z(),
490 vHAverNorm.x() * vVPre.y() - vHAverNorm.y() * vVPre.x());
491 #ifdef G4EVERBOSE
492 if(iverbose >= 2)
493 G4cout << " G4EP: HUPre " << vHUPre << " HVPre " << vHVPre << G4endl;
494 #endif
495
496 //------------------- COMPUTE MATRIX
497 //---------- 1/P
498
499 transf[0][0] =
500 1. -
501 deltaPInv * pAver *
502 (1. + (vpPostNorm.x() * deltaPos.x() + vpPostNorm.y() * deltaPos.y() +
503 vpPostNorm.z() * deltaPos.z()) /
504 stepLengthCm) +
505 2. * deltaPInv * pAver;
506
507 transf[0][1] =
508 -deltaPInv / thetaAver *
509 (TMSINT * gamma *
510 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
511 vHAverNorm.z() * vVPre.z()) +
512 sinThetaAver *
513 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
514 vVPre.z() * vpPostNorm.z()) +
515 OMcosThetaAver *
516 (vHVPre.x() * vpPostNorm.x() + vHVPre.y() * vpPostNorm.y() +
517 vHVPre.z() * vpPostNorm.z()));
518
519 transf[0][2] =
520 -sinpPre * deltaPInv / thetaAver *
521 (TMSINT * gamma *
522 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()) +
523 sinThetaAver *
524 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) +
525 OMcosThetaAver *
526 (vHUPre.x() * vpPostNorm.x() + vHUPre.y() * vpPostNorm.y() +
527 vHUPre.z() * vpPostNorm.z()));
528
529 transf[0][3] = -deltaPInv / stepLengthCm *
530 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y());
531
532 transf[0][4] = -deltaPInv / stepLengthCm *
533 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
534 vVPre.z() * vpPostNorm.z());
535
536 // *** Lambda
537 transf[1][0] =
538 -QP * ANV *
539 (vpPostNorm.x() * deltaPos.x() + vpPostNorm.y() * deltaPos.y() +
540 vpPostNorm.z() * deltaPos.z()) *
541 (1. + deltaPInv * pAver);
542 #ifdef G4EVERBOSE
543 if(iverbose >= 3)
544 G4cout << "ctransf10= " << transf[1][0] << " " << -QP << " " << ANV
545 << " " << vpPostNorm.x() << " " << deltaPos.x() << " "
546 << vpPostNorm.y() << " " << deltaPos.y() << " " << vpPostNorm.z()
547 << " " << deltaPos.z() << " " << deltaPInv << " " << pAver
548 << G4endl;
549 #endif
550
551 transf[1][1] =
552 cosThetaAver * (vVPre.x() * vVPost.x() + vVPre.y() * vVPost.y() +
553 vVPre.z() * vVPost.z()) +
554 sinThetaAver * (vHVPre.x() * vVPost.x() + vHVPre.y() * vVPost.y() +
555 vHVPre.z() * vVPost.z()) +
556 OMcosThetaAver *
557 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
558 vHAverNorm.z() * vVPre.z()) *
559 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
560 vHAverNorm.z() * vVPost.z()) +
561 ANV * (-sinThetaAver *
562 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
563 vVPre.z() * vpPostNorm.z()) +
564 OMcosThetaAver * (vVPre.x() * AN2.x() + vVPre.y() * AN2.y() +
565 vVPre.z() * AN2.z()) -
566 TMSINT * gamma *
567 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
568 vHAverNorm.z() * vVPre.z()));
569
570 transf[1][2] =
571 cosThetaAver * (vUPre.x() * vVPost.x() + vUPre.y() * vVPost.y()) +
572 sinThetaAver * (vHUPre.x() * vVPost.x() + vHUPre.y() * vVPost.y() +
573 vHUPre.z() * vVPost.z()) +
574 OMcosThetaAver *
575 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()) *
576 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
577 vHAverNorm.z() * vVPost.z()) +
578 ANV * (-sinThetaAver *
579 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) +
580 OMcosThetaAver * (vUPre.x() * AN2.x() + vUPre.y() * AN2.y()) -
581 TMSINT * gamma *
582 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()));
583 transf[1][2] = sinpPre * transf[1][2];
584
585 transf[1][3] = -QAver * ANV *
586 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y());
587
588 transf[1][4] = -QAver * ANV *
589 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
590 vVPre.z() * vpPostNorm.z());
591
592 // *** Phi
593
594 transf[2][0] =
595 -QP * ANU *
596 (vpPostNorm.x() * deltaPos.x() + vpPostNorm.y() * deltaPos.y() +
597 vpPostNorm.z() * deltaPos.z()) *
598 sinpPostInv * (1. + deltaPInv * pAver);
599 #ifdef G4EVERBOSE
600 if(iverbose >= 3)
601 G4cout << "ctransf20= " << transf[2][0] << " " << -QP << " " << ANU
602 << " " << vpPostNorm.x() << " " << deltaPos.x() << " "
603 << vpPostNorm.y() << " " << deltaPos.y() << " " << vpPostNorm.z()
604 << " " << deltaPos.z() << " " << sinpPostInv << " " << deltaPInv
605 << " " << pAver << G4endl;
606 #endif
607 transf[2][1] =
608 cosThetaAver * (vVPre.x() * vUPost.x() + vVPre.y() * vUPost.y()) +
609 sinThetaAver * (vHVPre.x() * vUPost.x() + vHVPre.y() * vUPost.y()) +
610 OMcosThetaAver *
611 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
612 vHAverNorm.z() * vVPre.z()) *
613 (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) +
614 ANU * (-sinThetaAver *
615 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
616 vVPre.z() * vpPostNorm.z()) +
617 OMcosThetaAver * (vVPre.x() * AN2.x() + vVPre.y() * AN2.y() +
618 vVPre.z() * AN2.z()) -
619 TMSINT * gamma *
620 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
621 vHAverNorm.z() * vVPre.z()));
622 transf[2][1] = sinpPostInv * transf[2][1];
623
624 transf[2][2] =
625 cosThetaAver * (vUPre.x() * vUPost.x() + vUPre.y() * vUPost.y()) +
626 sinThetaAver * (vHUPre.x() * vUPost.x() + vHUPre.y() * vUPost.y()) +
627 OMcosThetaAver *
628 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()) *
629 (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) +
630 ANU * (-sinThetaAver *
631 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) +
632 OMcosThetaAver * (vUPre.x() * AN2.x() + vUPre.y() * AN2.y()) -
633 TMSINT * gamma *
634 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y()));
635 transf[2][2] = sinpPostInv * sinpPre * transf[2][2];
636
637 transf[2][3] = -QAver * ANU *
638 (vUPre.x() * vpPostNorm.x() + vUPre.y() * vpPostNorm.y()) *
639 sinpPostInv;
640 #ifdef G4EVERBOSE
641 if(iverbose >= 3)
642 G4cout << "ctransf23= " << transf[2][3] << " " << -QAver << " " << ANU
643 << " " << vUPre.x() << " " << vpPostNorm.x() << " " << vUPre.y()
644 << " " << vpPostNorm.y() << " " << sinpPostInv << G4endl;
645 #endif
646
647 transf[2][4] = -QAver * ANU *
648 (vVPre.x() * vpPostNorm.x() + vVPre.y() * vpPostNorm.y() +
649 vVPre.z() * vpPostNorm.z()) *
650 sinpPostInv;
651
652 // *** Yt
653
654 transf[3][0] = pAver *
655 (vUPost.x() * deltaPos.x() + vUPost.y() * deltaPos.y()) *
656 (1. + deltaPInv * pAver);
657 #ifdef G4EVERBOSE
658 if(iverbose >= 3)
659 G4cout << "ctransf30= " << transf[3][0] << " " << pAver << " "
660 << vUPost.x() << " " << deltaPos.x() << " " << vUPost.y() << " "
661 << deltaPos.y() << " " << deltaPInv << " " << pAver << G4endl;
662 #endif
663
664 transf[3][1] =
665 (sinThetaAver * (vVPre.x() * vUPost.x() + vVPre.y() * vUPost.y()) +
666 OMcosThetaAver * (vHVPre.x() * vUPost.x() + vHVPre.y() * vUPost.y()) +
667 TMSINT * (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) *
668 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
669 vHAverNorm.z() * vVPre.z())) /
670 QAver;
671
672 transf[3][2] =
673 (sinThetaAver * (vUPre.x() * vUPost.x() + vUPre.y() * vUPost.y()) +
674 OMcosThetaAver * (vHUPre.x() * vUPost.x() + vHUPre.y() * vUPost.y()) +
675 TMSINT * (vHAverNorm.x() * vUPost.x() + vHAverNorm.y() * vUPost.y()) *
676 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y())) *
677 sinpPre / QAver;
678 #ifdef G4EVERBOSE
679 if(iverbose >= 3)
680 G4cout << "ctransf32= " << transf[3][2] << " " << sinThetaAver << " "
681 << vUPre.x() << " " << vUPost.x() << " " << vUPre.y() << " "
682 << vUPost.y() << " " << OMcosThetaAver << " " << vHUPre.x()
683 << " " << vUPost.x() << " " << vHUPre.y() << " " << vUPost.y()
684 << " " << TMSINT << " " << vHAverNorm.x() << " " << vUPost.x()
685 << " " << vHAverNorm.y() << " " << vUPost.y() << " "
686 << vHAverNorm.x() << " " << vUPre.x() << " " << vHAverNorm.y()
687 << " " << vUPre.y() << " " << sinpPre << " " << QAver << G4endl;
688 #endif
689
690 transf[3][3] = (vUPre.x() * vUPost.x() + vUPre.y() * vUPost.y());
691
692 transf[3][4] = (vVPre.x() * vUPost.x() + vVPre.y() * vUPost.y());
693
694 // *** Zt
695 transf[4][0] = pAver *
696 (vVPost.x() * deltaPos.x() + vVPost.y() * deltaPos.y() +
697 vVPost.z() * deltaPos.z()) *
698 (1. + deltaPInv * pAver);
699
700 transf[4][1] =
701 (sinThetaAver * (vVPre.x() * vVPost.x() + vVPre.y() * vVPost.y() +
702 vVPre.z() * vVPost.z()) +
703 OMcosThetaAver * (vHVPre.x() * vVPost.x() + vHVPre.y() * vVPost.y() +
704 vHVPre.z() * vVPost.z()) +
705 TMSINT *
706 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
707 vHAverNorm.z() * vVPost.z()) *
708 (vHAverNorm.x() * vVPre.x() + vHAverNorm.y() * vVPre.y() +
709 vHAverNorm.z() * vVPre.z())) /
710 QAver;
711 #ifdef G4EVERBOSE
712 if(iverbose >= 3)
713 G4cout << "ctransf41= " << transf[4][1] << " " << sinThetaAver << " "
714 << OMcosThetaAver << " " << TMSINT << " " << vVPre << " "
715 << vVPost << " " << vHVPre << " " << vHAverNorm << " " << QAver
716 << G4endl;
717 #endif
718
719 transf[4][2] =
720 (sinThetaAver * (vUPre.x() * vVPost.x() + vUPre.y() * vVPost.y()) +
721 OMcosThetaAver * (vHUPre.x() * vVPost.x() + vHUPre.y() * vVPost.y() +
722 vHUPre.z() * vVPost.z()) +
723 TMSINT *
724 (vHAverNorm.x() * vVPost.x() + vHAverNorm.y() * vVPost.y() +
725 vHAverNorm.z() * vVPost.z()) *
726 (vHAverNorm.x() * vUPre.x() + vHAverNorm.y() * vUPre.y())) *
727 sinpPre / QAver;
728
729 transf[4][3] = (vUPre.x() * vVPost.x() + vUPre.y() * vVPost.y());
730
731 transf[4][4] = (vVPre.x() * vVPost.x() + vVPre.y() * vVPost.y() +
732 vVPre.z() * vVPost.z());
733 // if(iverbose >= 3) G4cout <<"ctransf44= " << transf[4][4] <<" "<<
734 // vVPre.x() <<" "<<vVPost.x() <<" "<< vVPre.y() <<" "<< vVPost.y() <<"
735 // "<< vVPre.z() <<" "<< vVPost.z() << G4endl;
736
737 #ifdef G4EVERBOSE
738 if(iverbose >= 1)
739 G4cout << "G4EP: transf matrix computed " << transf << G4endl;
740 #endif
741 /* for( G4int ii=0;ii<5;ii++){
742 for( G4int jj=0;jj<5;jj++){
743 G4cout << transf[ii][jj] << " ";
744 }
745 G4cout << G4endl;
746 } */
747 }
748 }
749 // end of calculate transformation except it NEUTRAL PARTICLE OR FIELDFREE
750 // REGION
751 /* if( iverbose >= 1 ) G4cout << "G4EP: transf not updated but initialized "
752 << theFirstStep << G4endl; if( theFirstStep ) { theTransfMat = transf;
753 theFirstStep = false;
754 }else{
755 theTransfMat = theTransfMat * transf;
756 if( iverbose >= 1 ) G4cout << "G4EP: transf matrix accumulated" <<
757 theTransfMat << G4endl;
758 }
759 */
760 theTransfMat = std::move(transf);
761 #ifdef G4EVERBOSE
762 if(iverbose >= 1)
763 G4cout << "G4EP: error matrix before transformation " << fError << G4endl;
764 if(iverbose >= 2)
765 G4cout << " tf * err " << theTransfMat * fError << G4endl
766 << " transf matrix " << theTransfMat.T() << G4endl;
767 #endif
768
769 fError = fError.similarity(theTransfMat).T();
770 //- fError = transf * fError * transf.T();
771 #ifdef G4EVERBOSE
772 if(iverbose >= 1)
773 G4cout << "G4EP: error matrix propagated " << fError << G4endl;
774 #endif
775
776 //? S = B*S*BT S.similarity(B)
777 //? R = S
778 // not needed * *** TRANSFORM ERROR MATRIX FROM INTERNAL TO EXTERNAL
779 // VARIABLES;
780
781 PropagateErrorMSC(aTrack);
782
783 PropagateErrorIoni(aTrack);
784
785 return 0;
786 }
787
788 //------------------------------------------------------------------------
789 G4int G4ErrorFreeTrajState::PropagateErrorMSC(const G4Track* aTrack)
790 {
791 G4ThreeVector vpPre = aTrack->GetMomentum() / GeV;
792 G4double pPre = vpPre.mag();
793 G4double pBeta = pPre * pPre / (aTrack->GetTotalEnergy() / GeV);
794 G4double stepLengthCm = aTrack->GetStep()->GetStepLength() / cm;
795
796 G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
797 G4double effZ, effA;
798 CalculateEffectiveZandA(mate, effZ, effA);
799
800 #ifdef G4EVERBOSE
801 if(iverbose >= 4)
802 G4cout << "material "
803 << mate->GetName()
804 //<< " " << mate->GetZ() << " " << mate->GetA()
805 << " effZ:" << effZ << " effA:" << effA
806 << " dens(g/mole):" << mate->GetDensity() / g * mole
807 << " Radlen/cm:" << mate->GetRadlen() / cm << " nuclLen/cm"
808 << mate->GetNuclearInterLength() / cm << G4endl;
809 #endif
810
811 G4double RI = stepLengthCm / (mate->GetRadlen() / cm);
812 #ifdef G4EVERBOSE
813 if(iverbose >= 4)
814 G4cout << std::setprecision(6) << std::setw(6) << "G4EP:MSC: RI=X/X0 " << RI
815 << " stepLengthCm " << stepLengthCm << " radlen/cm "
816 << (mate->GetRadlen() / cm) << " RI*1.e10:" << RI * 1.e10 << G4endl;
817 #endif
818 G4double charge = aTrack->GetDynamicParticle()->GetCharge();
819 G4double DD = 1.8496E-4 * RI * (charge / pBeta * charge / pBeta);
820 #ifdef G4EVERBOSE
821 if(iverbose >= 3)
822 G4cout << "G4EP:MSC: D*1E6= " << DD * 1.E6 << " pBeta " << pBeta << G4endl;
823 #endif
824 G4double S1 = DD * stepLengthCm * stepLengthCm / 3.;
825 G4double S2 = DD;
826 G4double S3 = DD * stepLengthCm / 2.;
827
828 G4double CLA =
829 std::sqrt(vpPre.x() * vpPre.x() + vpPre.y() * vpPre.y()) / pPre;
830 #ifdef G4EVERBOSE
831 if(iverbose >= 2)
832 G4cout << std::setw(6) << "G4EP:MSC: RI " << RI << " S1 " << S1 << " S2 "
833 << S2 << " S3 " << S3 << " CLA " << CLA << G4endl;
834 #endif
835 fError[1][1] += S2;
836 fError[1][4] -= S3;
837 fError[2][2] += S2 / CLA / CLA;
838 fError[2][3] += S3 / CLA;
839 fError[3][3] += S1;
840 fError[4][4] += S1;
841
842 #ifdef G4EVERBOSE
843 if(iverbose >= 2)
844 G4cout << "G4EP:MSC: error matrix propagated msc " << fError << G4endl;
845 #endif
846
847 return 0;
848 }
849
850 //------------------------------------------------------------------------
851 void G4ErrorFreeTrajState::CalculateEffectiveZandA(const G4Material* mate,
852 G4double& effZ,
853 G4double& effA)
854 {
855 effZ = 0.;
856 effA = 0.;
857 auto nelem = mate->GetNumberOfElements();
858 const G4double* fracVec = mate->GetFractionVector();
859 for(G4int ii = 0; ii < (G4int)nelem; ++ii)
860 {
861 effZ += mate->GetElement(ii)->GetZ() * fracVec[ii];
862 effA += mate->GetElement(ii)->GetA() * fracVec[ii] / g * mole;
863 }
864 }
865
866 //------------------------------------------------------------------------
867 G4int G4ErrorFreeTrajState::PropagateErrorIoni(const G4Track* aTrack)
868 {
869 G4double stepLengthCm = aTrack->GetStep()->GetStepLength() / cm;
870 #ifdef G4EVERBOSE
871 G4double DEDX2;
872 if(stepLengthCm < 1.E-7)
873 {
874 DEDX2 = 0.;
875 }
876 #endif
877 // * Calculate xi factor (KeV).
878 G4Material* mate = aTrack->GetVolume()->GetLogicalVolume()->GetMaterial();
879 G4double effZ, effA;
880 CalculateEffectiveZandA(mate, effZ, effA);
881
882 G4double Etot = aTrack->GetTotalEnergy() / GeV;
883 G4double beta = aTrack->GetMomentum().mag() / GeV / Etot;
884 G4double mass = aTrack->GetDynamicParticle()->GetMass() / GeV;
885 G4double gamma = Etot / mass;
886
887 // * Calculate xi factor (keV).
888 G4double XI = 153.5 * effZ * stepLengthCm * (mate->GetDensity() / mg * mole) /
889 (effA * beta * beta);
890
891 #ifdef G4EVERBOSE
892 if(iverbose >= 2)
893 {
894 G4cout << "G4EP:IONI: XI/keV " << XI << " beta " << beta << " gamma "
895 << gamma << G4endl;
896 G4cout << " density " << (mate->GetDensity() / mg * mole) << " effA "
897 << effA << " step " << stepLengthCm << G4endl;
898 }
899 #endif
900 // * Maximum energy transfer to atomic electron (KeV).
901 G4double eta = beta * gamma;
902 G4double etasq = eta * eta;
903 G4double eMass = 0.51099906 / GeV;
904 G4double massRatio = eMass / mass;
905 G4double F1 = 2 * eMass * etasq;
906 G4double F2 = 1. + 2. * massRatio * gamma + massRatio * massRatio;
907 G4double Emax = 1.E+6 * F1 / F2; // now in keV
908
909 // * *** and now sigma**2 in GeV
910 G4double dedxSq =
911 XI * Emax * (1. - (beta * beta / 2.)) * 1.E-12; // now in GeV^2
912 /*The above formula for var(1/p) good for dens scatterers. However, for MIPS
913 passing through a gas it leads to overestimation. Further more for incident
914 electrons the Emax is almost equal to incident energy. This leads to
915 k=Xi/Emax as small as e-6 and gradually the cov matrix explodes.
916
917 http://www2.pv.infn.it/~rotondi/kalman_1.pdf
918
919 Since I do not have enough info at the moment to implement Landau &
920 sub-Landau models for k=Xi/Emax <0.01 I'll saturate k at this value for now
921 */
922
923 if(XI / Emax < 0.01)
924 dedxSq *=
925 XI / Emax * 100; // Quench for low Elos, see above: newVar=odVar *k/0.01
926
927 #ifdef G4EVERBOSE
928 if(iverbose >= 2)
929 G4cout << "G4EP:IONI: DEDX^2(GeV^2) " << dedxSq << " emass/GeV: " << eMass
930 << " Emax/keV: " << Emax << " k=Xi/Emax=" << XI / Emax << G4endl;
931
932 #endif
933
934 G4double pPre6 =
935 (aTrack->GetStep()->GetPreStepPoint()->GetMomentum() / GeV).mag();
936 pPre6 = std::pow(pPre6, 6);
937 // Apply it to error
938 fError[0][0] += Etot * Etot * dedxSq / pPre6;
939 #ifdef G4EVERBOSE
940 if(iverbose >= 2)
941 G4cout << "G4:IONI Etot/GeV: " << Etot << " err_dedx^2/GeV^2: " << dedxSq
942 << " p^6: " << pPre6 << G4endl;
943 if(iverbose >= 2)
944 G4cout << "G4EP:IONI: error2_from_ionisation "
945 << (Etot * Etot * dedxSq) / pPre6 << G4endl;
946 #endif
947
948 return 0;
949 }
950