| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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(con
49 con
50 con
51 con
52 : G4ErrorTrajState(partName, pos, mom, errma
53 {
54 fTrajParam = G4ErrorFreeTrajParam(pos, mom);
55 Init();
56 }
57
58 //--------------------------------------------
59 G4ErrorFreeTrajState::G4ErrorFreeTrajState(con
60 : G4ErrorTrajState(tpSD.GetParticleType(), t
61 tpSD.GetMomentum())
62 {
63 // G4ThreeVector planeNormal = tpSD.GetPlane
64 // G4double fPt = tpSD.GetMomentum()*planeNo
65 // plane G4ErrorSurfaceTrajParam tpSDparam =
66 // G4ThreeVector Psc = fPt * planeNormal +
67 // tpSDparam.GetPU()*tpSDparam.GetVectorU()
68
69 fTrajParam = G4ErrorFreeTrajParam(fPosition,
70 Init();
71
72 //----- Get the error matrix in SC coordinat
73 G4ErrorSurfaceTrajParam tpSDparam = tpSD.Get
74 G4double mom = fMomentu
75 G4double mom2 = fMomentu
76 G4double TVW1 =
77 std::sqrt(mom2 / (mom2 + tpSDparam.GetPV()
78 tpSDparam.GetPW() * tpSD
79 G4ThreeVector vTVW(TVW1, tpSDparam.GetPV() /
80 tpSDparam.GetPW() / mom *
81 G4Vector3D vectorU = tpSDparam.GetVectorV().
82 G4Vector3D vTN = vTVW.x() * vectorU + vT
83 vTVW.z() * tpSDparam.GetVec
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 << " = " <
92 << " diff " << pc2 - GetParameters(
93 G4cout << " CHECK: pc3 " << pc3 << " = " <
94 << " diff " << pc3 - GetParameters(
95 }
96 #endif
97
98 //--- Get the unit vectors perp to P
99 G4double cosl = std::cos(GetParameters().Get
100 if(cosl < 1.E-30)
101 cosl = 1.E-30;
102 G4double cosl1 = 1. / cosl;
103 G4Vector3D vUN(-vTN.y() * cosl1, vTN.x() * c
104 G4Vector3D vVN(-vTN.z() * vUN.y(), vTN.z() *
105
106 G4Vector3D vUperp = G4Vector3D(-fMomentum.y(
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 << " = " <
115 << (vUN - vUperp).mag() << G4endl;
116 G4cout << " CHECK: vVN " << vVN << " = " <
117 << (vVN - vVperp).mag() << G4endl;
118 }
119 #endif
120
121 // get the dot products of vectors perpendic
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 = G4TransportationManag
132 ->GetFieldManager()
133 ->GetDetectorField(
134 G4ThreeVector dir = fTrajParam.GetDirecti
135 G4double invCosTheta = 1. / std::cos(dir.the
136 G4cout << " dir=" << dir << " invCosTheta "
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],
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 * magHPre
156 G4double cosz = HPre * vVperp * magHPreM
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 * in
161 transfM[2][4] = -Q * dir.z() * cosz * in
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().simi
176
177 #ifdef G4EVERBOSE
178 if(iverbose >= 1)
179 G4cout << "error matrix SD2SC " << fError
180 if(iverbose >= 4)
181 G4cout << "G4ErrorFreeTrajState from SD "
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&
196
197 //--------------------------------------------
198 G4int G4ErrorFreeTrajState::Update(const G4Tra
199 {
200 G4int ierr = 0;
201 fTrajParam.Update(aTrack);
202 UpdatePosMom(aTrack->GetPosition(), aTrack->
203 return ierr;
204 }
205
206 //--------------------------------------------
207 std::ostream& operator<<(std::ostream& out, co
208 {
209 std::ios::fmtflags orig_flags = out.flags();
210
211 out.setf(std::ios::fixed, std::ios::floatfie
212
213 ts.DumpPosMomError(out);
214
215 out << " G4ErrorFreeTrajState: Params: " <<
216
217 out.flags(orig_flags);
218
219 return out;
220 }
221
222 //--------------------------------------------
223 G4int G4ErrorFreeTrajState::PropagateError(con
224 {
225 G4double stepLengthCm = aTrack->GetStep()->G
226 if(G4ErrorPropagatorData::GetErrorPropagator
227 G4ErrorStage_Deflation)
228 stepLengthCm *= -1.;
229
230 G4double kCarTolerance =
231 G4GeometryTolerance::GetInstance()->GetSur
232
233 if(std::fabs(stepLengthCm) <= kCarTolerance
234 return 0;
235
236 #ifdef G4EVERBOSE
237 if(iverbose >= 2)
238 G4cout << " G4ErrorFreeTrajState::Propaga
239 G4cout << "G4EP: iverbose=" << iverbose << G
240 #endif
241
242 // * *** ERROR PROPAGATION ON A HELIX ASSUMI
243 G4Point3D vposPost = aTrack->GetPosition() /
244 G4Vector3D vpPost = aTrack->GetMomentum() /
245 // G4Point3D vposPre = fPosition/cm;
246 // G4Vector3D vpPre = fMomentum/GeV;
247 G4Point3D vposPre = aTrack->GetStep()->GetPr
248 G4Vector3D vpPre = aTrack->GetStep()->GetPr
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 << G
261 << vposPost << G4endl;
262 G4cout << "G4EP: vpPre " << vpPre << G4end
263 << G4endl;
264 G4cout << " err start step " << fError <<
265 G4cout << "G4EP: stepLengthCm " << stepLen
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 << "
276 << " deltaPInv:" << deltaPInv << G
277
278 G4Vector3D vpPreNorm = vpPre * pInvPre;
279 G4Vector3D vpPostNorm = vpPost * pInvPost;
280 if(iverbose >= 2)
281 G4cout << "G4EP: vpPreNorm " << vpPreNorm
282 << G4endl;
283 // return if propagation along Z??
284 if(1. - std::fabs(vpPreNorm.z()) < kCarToler
285 return 4;
286 if(1. - std::fabs(vpPostNorm.z()) < kCarTole
287 return 4;
288 G4double sinpPre =
289 std::sin(vpPreNorm.theta()); // cosine pe
290 G4double sinpPost =
291 std::sin(vpPostNorm.theta()); // cosine p
292 G4double sinpPostInv = 1. / std::sin(vpPostN
293
294 #ifdef G4EVERBOSE
295 if(iverbose >= 2)
296 G4cout << "G4EP: cosl " << sinpPre << " co
297 #endif
298 //* *** DEFINE TRANSFORMATION MATRIX BETWEEN
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 " <
312 }
313 #endif
314
315 // charge X propagation direction
316 G4double charge = aTrack->GetDynamicParticle
317 if(G4ErrorPropagatorData::GetErrorPropagator
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
328 // eventually it will not be transformed) it
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 = G4TransportationManag
341 ->GetFieldManager()
342 ->GetDetectorField(
343 if(!field)
344 return 0; // goto 45
345
346 // calculate transformation except it NEUTRA
347 if(charge != 0. && field)
348 {
349 field->GetFieldValue(pos1,
350 h1); // here pos1[],
351 field->GetFieldValue(pos2, h2);
352 G4ThreeVector HPre =
353 G4ThreeVector(h1[0], h1[1], h1[2]) / tes
354 10.; // 10. is to get same dimensions a
355 G4ThreeVector HPost = G4ThreeVector(h2[0],
356 G4double magHPre = HPre.mag();
357 G4double magHPost = HPost.mag();
358 #ifdef G4EVERBOSE
359 if(iverbose >= 2)
360 {
361 G4cout << "G4EP: h1 = " << h1[0] << ", "
362 << G4endl;
363 G4cout << "G4EP: pos1/mm = " << pos1[0]
364 << pos1[2] << G4endl;
365 G4cout << "G4EP: pos2/mm = " << pos2[0]
366 << pos2[2] << G4endl;
367 G4cout << "G4EP: B-filed in KGauss HPre
368 << "G4EP: in KGauss HPost " << HP
369 }
370 #endif
371
372 if(magHPre + magHPost != 0.)
373 {
374 //* *** CHECK WHETHER H*ALFA/P IS TOO DI
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 t
386 // line trajectory
387 G4double alphaSqr = 1. - gam * gam;
388 G4double diffHSqr = (HPre * pInvPre - H
389 G4double delhp6Sqr = 300. * 300.;
390 #ifdef G4EVERBOSE
391 if(iverbose >= 2)
392 {
393 G4cout << " G4EP: gam " << gam << " al
394 << " diffHSqr " << diffHSqr <<
395 G4cout << " alpha= " << std::sqrt(alph
396 }
397 #endif
398 if(diffHSqr * alphaSqr > delhp6Sqr)
399 return 3;
400
401 //* *** DEFINE AVERAGE MAGNETIC FIELD AN
402 G4double pInvAver = 1. / (pInvPre + pInv
403 G4double CFACT8 = 2.997925E-4;
404 // G4double HAver
405 G4ThreeVector vHAverNorm((HPre * pInvPre
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 " << vHAve
413 << " charge " << charge << G4en
414 #endif
415
416 G4double pAver = (pPre + pPost) *
417 G4double QAver = -HAver / pAver;
418 G4double thetaAver = QAver * stepLeng
419 G4double sinThetaAver = std::sin(thetaAv
420 G4double cosThetaAver = std::cos(thetaAv
421 G4double gamma = vHAverNorm * vpP
422 G4ThreeVector AN2 = vHAverNorm.cross
423
424 #ifdef G4EVERBOSE
425 if(iverbose >= 2)
426 G4cout << " G4EP: AN2 " << AN2 << " g
427 << " theta=" << thetaAver << G
428 #endif
429 G4double AU = 1. / vpPreNorm.perp();
430 // t G4ThreeVector vU( vpPreNorm.cross(
431 G4ThreeVector vUPre(-AU * vpPreNorm.y(),
432 G4ThreeVector vVPre(-vpPreNorm.z() * vUP
433 vpPreNorm.x() * vUPr
434 vpPreNorm.y() * vU
435
436 //
437 AU = 1. / vpPostNorm.perp();
438 // t G4ThreeVector vU( vpPostNorm.cross
439 // );
440 G4ThreeVector vUPost(-AU * vpPostNorm.y(
441 G4ThreeVector vVPost(
442 -vpPostNorm.z() * vUPost.y(), vpPostNo
443 vpPostNorm.x() * vUPost.y() - vpPostNo
444 #ifdef G4EVERBOSE
445 G4cout << " vpPostNorm " << vpPostNorm <
446 if(iverbose >= 2)
447 G4cout << " G4EP: AU " << AU << " vUPr
448 << " vUPost " << vUPost << " vV
449 #endif
450 G4Point3D deltaPos(vposPre - vposPost);
451
452 // * *** COMPLETE TRANSFORMATION MATRIX
453 // * *** FIELD GRADIENT PERPENDICULAR TO
454 // * *** TAKEN INTO ACCOUNT
455
456 G4double QP = QAver * pAver; // = -HAve
457 #ifdef G4EVERBOSE
458 if(iverbose >= 2)
459 G4cout << " G4EP: QP " << QP << " QAve
460 << G4endl;
461 #endif
462 G4double ANV =
463 -(vHAverNorm.x() * vUPost.x() + vHAver
464 G4double ANU =
465 (vHAverNorm.x() * vVPost.x() + vHAverN
466 vHAverNorm.z() * vVPost.z());
467 G4double OMcosThetaAver = 1. - cosThetaA
468 #ifdef G4EVERBOSE
469 if(iverbose >= 2)
470 G4cout << "G4EP: OMcosThetaAver " << O
471 << cosThetaAver << " thetaAver
472 << QAver << " stepLengthCm " <<
473 #endif
474 G4double TMSINT = thetaAver - sinThetaAv
475 #ifdef G4EVERBOSE
476 if(iverbose >= 2)
477 G4cout << " G4EP: ANV " << ANV << " AN
478 #endif
479
480 G4ThreeVector vHUPre(
481 -vHAverNorm.z() * vUPre.y(), vHAverNor
482 vHAverNorm.x() * vUPre.y() - vHAverNor
483 #ifdef G4EVERBOSE
484 // if( iverbose >= 2 ) G4cout << "G4E
485 // << vHAverNorm.z() << " " << vUPre.
486 #endif
487 G4ThreeVector vHVPre(
488 vHAverNorm.y() * vVPre.z() - vHAverNor
489 vHAverNorm.z() * vVPre.x() - vHAverNor
490 vHAverNorm.x() * vVPre.y() - vHAverNor
491 #ifdef G4EVERBOSE
492 if(iverbose >= 2)
493 G4cout << " G4EP: HUPre " << vHUPre <<
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()
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() + vHAve
511 vHAverNorm.z() * vVPre.z()) +
512 sinThetaAver *
513 (vVPre.x() * vpPostNorm.x() + vVPre
514 vVPre.z() * vpPostNorm.z()) +
515 OMcosThetaAver *
516 (vHVPre.x() * vpPostNorm.x() + vHVP
517 vHVPre.z() * vpPostNorm.z()));
518
519 transf[0][2] =
520 -sinpPre * deltaPInv / thetaAver *
521 (TMSINT * gamma *
522 (vHAverNorm.x() * vUPre.x() + vHAve
523 sinThetaAver *
524 (vUPre.x() * vpPostNorm.x() + vUPre
525 OMcosThetaAver *
526 (vHUPre.x() * vpPostNorm.x() + vHUP
527 vHUPre.z() * vpPostNorm.z()));
528
529 transf[0][3] = -deltaPInv / stepLengthCm
530 (vUPre.x() * vpPostNorm.x
531
532 transf[0][4] = -deltaPInv / stepLengthCm
533 (vVPre.x() * vpPostNorm.x
534 vVPre.z() * vpPostNorm.z
535
536 // *** Lambda
537 transf[1][0] =
538 -QP * ANV *
539 (vpPostNorm.x() * deltaPos.x() + vpPos
540 vpPostNorm.z() * deltaPos.z()) *
541 (1. + deltaPInv * pAver);
542 #ifdef G4EVERBOSE
543 if(iverbose >= 3)
544 G4cout << "ctransf10= " << transf[1][0
545 << " " << vpPostNorm.x() << " "
546 << vpPostNorm.y() << " " << del
547 << " " << deltaPos.z() << " " <
548 << G4endl;
549 #endif
550
551 transf[1][1] =
552 cosThetaAver * (vVPre.x() * vVPost.x()
553 vVPre.z() * vVPost.z()
554 sinThetaAver * (vHVPre.x() * vVPost.x(
555 vHVPre.z() * vVPost.z(
556 OMcosThetaAver *
557 (vHAverNorm.x() * vVPre.x() + vHAver
558 vHAverNorm.z() * vVPre.z()) *
559 (vHAverNorm.x() * vVPost.x() + vHAve
560 vHAverNorm.z() * vVPost.z()) +
561 ANV * (-sinThetaAver *
562 (vVPre.x() * vpPostNorm.x() +
563 vVPre.z() * vpPostNorm.z())
564 OMcosThetaAver * (vVPre.x() * A
565 vVPre.z() * A
566 TMSINT * gamma *
567 (vHAverNorm.x() * vVPre.x() +
568 vHAverNorm.z() * vVPre.z()))
569
570 transf[1][2] =
571 cosThetaAver * (vUPre.x() * vVPost.x()
572 sinThetaAver * (vHUPre.x() * vVPost.x(
573 vHUPre.z() * vVPost.z(
574 OMcosThetaAver *
575 (vHAverNorm.x() * vUPre.x() + vHAver
576 (vHAverNorm.x() * vVPost.x() + vHAve
577 vHAverNorm.z() * vVPost.z()) +
578 ANV * (-sinThetaAver *
579 (vUPre.x() * vpPostNorm.x() +
580 OMcosThetaAver * (vUPre.x() * A
581 TMSINT * gamma *
582 (vHAverNorm.x() * vUPre.x() +
583 transf[1][2] = sinpPre * transf[1][2];
584
585 transf[1][3] = -QAver * ANV *
586 (vUPre.x() * vpPostNorm.x
587
588 transf[1][4] = -QAver * ANV *
589 (vVPre.x() * vpPostNorm.x
590 vVPre.z() * vpPostNorm.z
591
592 // *** Phi
593
594 transf[2][0] =
595 -QP * ANU *
596 (vpPostNorm.x() * deltaPos.x() + vpPos
597 vpPostNorm.z() * deltaPos.z()) *
598 sinpPostInv * (1. + deltaPInv * pAver)
599 #ifdef G4EVERBOSE
600 if(iverbose >= 3)
601 G4cout << "ctransf20= " << transf[2][0
602 << " " << vpPostNorm.x() << " "
603 << vpPostNorm.y() << " " << del
604 << " " << deltaPos.z() << " " <
605 << " " << pAver << G4endl;
606 #endif
607 transf[2][1] =
608 cosThetaAver * (vVPre.x() * vUPost.x()
609 sinThetaAver * (vHVPre.x() * vUPost.x(
610 OMcosThetaAver *
611 (vHAverNorm.x() * vVPre.x() + vHAver
612 vHAverNorm.z() * vVPre.z()) *
613 (vHAverNorm.x() * vUPost.x() + vHAve
614 ANU * (-sinThetaAver *
615 (vVPre.x() * vpPostNorm.x() +
616 vVPre.z() * vpPostNorm.z())
617 OMcosThetaAver * (vVPre.x() * A
618 vVPre.z() * A
619 TMSINT * gamma *
620 (vHAverNorm.x() * vVPre.x() +
621 vHAverNorm.z() * vVPre.z()))
622 transf[2][1] = sinpPostInv * transf[2][1
623
624 transf[2][2] =
625 cosThetaAver * (vUPre.x() * vUPost.x()
626 sinThetaAver * (vHUPre.x() * vUPost.x(
627 OMcosThetaAver *
628 (vHAverNorm.x() * vUPre.x() + vHAver
629 (vHAverNorm.x() * vUPost.x() + vHAve
630 ANU * (-sinThetaAver *
631 (vUPre.x() * vpPostNorm.x() +
632 OMcosThetaAver * (vUPre.x() * A
633 TMSINT * gamma *
634 (vHAverNorm.x() * vUPre.x() +
635 transf[2][2] = sinpPostInv * sinpPre * t
636
637 transf[2][3] = -QAver * ANU *
638 (vUPre.x() * vpPostNorm.x
639 sinpPostInv;
640 #ifdef G4EVERBOSE
641 if(iverbose >= 3)
642 G4cout << "ctransf23= " << transf[2][3
643 << " " << vUPre.x() << " " << v
644 << " " << vpPostNorm.y() << " "
645 #endif
646
647 transf[2][4] = -QAver * ANU *
648 (vVPre.x() * vpPostNorm.x
649 vVPre.z() * vpPostNorm.z
650 sinpPostInv;
651
652 // *** Yt
653
654 transf[3][0] = pAver *
655 (vUPost.x() * deltaPos.x(
656 (1. + deltaPInv * pAver);
657 #ifdef G4EVERBOSE
658 if(iverbose >= 3)
659 G4cout << "ctransf30= " << transf[3][0
660 << vUPost.x() << " " << deltaPo
661 << deltaPos.y() << " " << delta
662 #endif
663
664 transf[3][1] =
665 (sinThetaAver * (vVPre.x() * vUPost.x(
666 OMcosThetaAver * (vHVPre.x() * vUPost
667 TMSINT * (vHAverNorm.x() * vUPost.x()
668 (vHAverNorm.x() * vVPre.x() + vHAve
669 vHAverNorm.z() * vVPre.z())) /
670 QAver;
671
672 transf[3][2] =
673 (sinThetaAver * (vUPre.x() * vUPost.x(
674 OMcosThetaAver * (vHUPre.x() * vUPost
675 TMSINT * (vHAverNorm.x() * vUPost.x()
676 (vHAverNorm.x() * vUPre.x() + vHAve
677 sinpPre / QAver;
678 #ifdef G4EVERBOSE
679 if(iverbose >= 3)
680 G4cout << "ctransf32= " << transf[3][2
681 << vUPre.x() << " " << vUPost.x
682 << vUPost.y() << " " << OMcosTh
683 << " " << vUPost.x() << " " <<
684 << " " << TMSINT << " " << vHAv
685 << " " << vHAverNorm.y() << " "
686 << vHAverNorm.x() << " " << vUP
687 << " " << vUPre.y() << " " << s
688 #endif
689
690 transf[3][3] = (vUPre.x() * vUPost.x() +
691
692 transf[3][4] = (vVPre.x() * vUPost.x() +
693
694 // *** Zt
695 transf[4][0] = pAver *
696 (vVPost.x() * deltaPos.x(
697 vVPost.z() * deltaPos.z(
698 (1. + deltaPInv * pAver);
699
700 transf[4][1] =
701 (sinThetaAver * (vVPre.x() * vVPost.x(
702 vVPre.z() * vVPost.z(
703 OMcosThetaAver * (vHVPre.x() * vVPost
704 vHVPre.z() * vVPost
705 TMSINT *
706 (vHAverNorm.x() * vVPost.x() + vHAv
707 vHAverNorm.z() * vVPost.z()) *
708 (vHAverNorm.x() * vVPre.x() + vHAve
709 vHAverNorm.z() * vVPre.z())) /
710 QAver;
711 #ifdef G4EVERBOSE
712 if(iverbose >= 3)
713 G4cout << "ctransf41= " << transf[4][1
714 << OMcosThetaAver << " " << TMS
715 << vVPost << " " << vHVPre << "
716 << G4endl;
717 #endif
718
719 transf[4][2] =
720 (sinThetaAver * (vUPre.x() * vVPost.x(
721 OMcosThetaAver * (vHUPre.x() * vVPost
722 vHUPre.z() * vVPost
723 TMSINT *
724 (vHAverNorm.x() * vVPost.x() + vHAv
725 vHAverNorm.z() * vVPost.z()) *
726 (vHAverNorm.x() * vUPre.x() + vHAve
727 sinpPre / QAver;
728
729 transf[4][3] = (vUPre.x() * vVPost.x() +
730
731 transf[4][4] = (vVPre.x() * vVPost.x() +
732 vVPre.z() * vVPost.z());
733 // if(iverbose >= 3) G4cout <<"ctransf
734 // vVPre.x() <<" "<<vVPost.x() <<" "<
735 // "<< vVPre.z() <<" "<< vVPost.z() <<
736
737 #ifdef G4EVERBOSE
738 if(iverbose >= 1)
739 G4cout << "G4EP: transf matrix compute
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
750 // REGION
751 /* if( iverbose >= 1 ) G4cout << "G4EP: tra
752 << theFirstStep << G4endl; if( theFirstStep
753 theFirstStep = false;
754 }else{
755 theTransfMat = theTransfMat * transf;
756 if( iverbose >= 1 ) G4cout << "G4EP: trans
757 theTransfMat << G4endl;
758 }
759 */
760 theTransfMat = std::move(transf);
761 #ifdef G4EVERBOSE
762 if(iverbose >= 1)
763 G4cout << "G4EP: error matrix before trans
764 if(iverbose >= 2)
765 G4cout << " tf * err " << theTransfMat * f
766 << " transf matrix " << theTransfMa
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 "
774 #endif
775
776 //? S = B*S*BT S.similarity(B)
777 //? R = S
778 // not needed * *** TRANSFORM ERROR MATRIX F
779 // VARIABLES;
780
781 PropagateErrorMSC(aTrack);
782
783 PropagateErrorIoni(aTrack);
784
785 return 0;
786 }
787
788 //--------------------------------------------
789 G4int G4ErrorFreeTrajState::PropagateErrorMSC(
790 {
791 G4ThreeVector vpPre = aTrack->GetMomentum(
792 G4double pPre = vpPre.mag();
793 G4double pBeta = pPre * pPre / (aTrac
794 G4double stepLengthCm = aTrack->GetStep()->G
795
796 G4Material* mate = aTrack->GetVolume()->GetL
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() << " " <<
805 << " effZ:" << effZ << " effA:" <<
806 << " dens(g/mole):" << mate->GetDen
807 << " Radlen/cm:" << mate->GetRadlen
808 << mate->GetNuclearInterLength() /
809 #endif
810
811 G4double RI = stepLengthCm / (mate->GetRadle
812 #ifdef G4EVERBOSE
813 if(iverbose >= 4)
814 G4cout << std::setprecision(6) << std::set
815 << " stepLengthCm " << stepLengthCm
816 << (mate->GetRadlen() / cm) << " RI
817 #endif
818 G4double charge = aTrack->GetDynamicParticle
819 G4double DD = 1.8496E-4 * RI * (charge /
820 #ifdef G4EVERBOSE
821 if(iverbose >= 3)
822 G4cout << "G4EP:MSC: D*1E6= " << DD * 1.E6
823 #endif
824 G4double S1 = DD * stepLengthCm * stepLength
825 G4double S2 = DD;
826 G4double S3 = DD * stepLengthCm / 2.;
827
828 G4double CLA =
829 std::sqrt(vpPre.x() * vpPre.x() + vpPre.y(
830 #ifdef G4EVERBOSE
831 if(iverbose >= 2)
832 G4cout << std::setw(6) << "G4EP:MSC: RI "
833 << S2 << " S3 " << S3 << " CLA " <<
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 propagat
845 #endif
846
847 return 0;
848 }
849
850 //--------------------------------------------
851 void G4ErrorFreeTrajState::CalculateEffectiveZ
852
853
854 {
855 effZ = 0.;
856 effA = 0.;
857 auto nelem = mate->GetNumberOfElements();
858 const G4double* fracVec = mate->GetFractionV
859 for(G4int ii = 0; ii < (G4int)nelem; ++ii)
860 {
861 effZ += mate->GetElement(ii)->GetZ() * fra
862 effA += mate->GetElement(ii)->GetA() * fra
863 }
864 }
865
866 //--------------------------------------------
867 G4int G4ErrorFreeTrajState::PropagateErrorIoni
868 {
869 G4double stepLengthCm = aTrack->GetStep()->G
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()->GetL
879 G4double effZ, effA;
880 CalculateEffectiveZandA(mate, effZ, effA);
881
882 G4double Etot = aTrack->GetTotalEnergy() /
883 G4double beta = aTrack->GetMomentum().mag()
884 G4double mass = aTrack->GetDynamicParticle(
885 G4double gamma = Etot / mass;
886
887 // * Calculate xi factor (keV).
888 G4double XI = 153.5 * effZ * stepLengthCm *
889 (effA * beta * beta);
890
891 #ifdef G4EVERBOSE
892 if(iverbose >= 2)
893 {
894 G4cout << "G4EP:IONI: XI/keV " << XI << "
895 << gamma << G4endl;
896 G4cout << " density " << (mate->GetDensity
897 << effA << " step " << stepLengthCm
898 }
899 #endif
900 // * Maximum energy transfer to atomic e
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 * g
907 G4double Emax = 1.E+6 * F1 / F2; // no
908
909 // * *** and now sigma**2 in GeV
910 G4double dedxSq =
911 XI * Emax * (1. - (beta * beta / 2.)) * 1.
912 /*The above formula for var(1/p) good for d
913 passing through a gas it leads to overesti
914 electrons the Emax is almost equal to inci
915 k=Xi/Emax as small as e-6 and gradually
916
917 http://www2.pv.infn.it/~rotondi/kalman_1.p
918
919 Since I do not have enough info at the mom
920 sub-Landau models for k=Xi/Emax <0.01 I'll
921 */
922
923 if(XI / Emax < 0.01)
924 dedxSq *=
925 XI / Emax * 100; // Quench for low Elos
926
927 #ifdef G4EVERBOSE
928 if(iverbose >= 2)
929 G4cout << "G4EP:IONI: DEDX^2(GeV^2) " << d
930 << " Emax/keV: " << Emax << " k=Xi
931
932 #endif
933
934 G4double pPre6 =
935 (aTrack->GetStep()->GetPreStepPoint()->Get
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 <<
942 << " p^6: " << pPre6 << G4endl;
943 if(iverbose >= 2)
944 G4cout << "G4EP:IONI: error2_from_ionisati
945 << (Etot * Etot * dedxSq) / pPre6 <
946 #endif
947
948 return 0;
949 }
950