| 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 // G4ChordFinder implementation
27 //
28 // Author: J.Apostolakis - Design and implemen
29 // -------------------------------------------
30
31 #include <iomanip>
32
33 #include "G4ChordFinder.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4MagneticField.hh"
36 #include "G4Mag_UsualEqRhs.hh"
37 #include "G4MagIntegratorDriver.hh"
38 // #include "G4ClassicalRK4.hh"
39 // #include "G4CashKarpRKF45.hh"
40 // #include "G4NystromRK4.hh"
41 // #include "G4BogackiShampine23.hh"
42 // #include "G4BogackiShampine45.hh"
43
44 #include "G4DormandPrince745.hh"
45
46 // New templated stepper(s) -- avoid virtual c
47 #include "G4TDormandPrince45.hh"
48
49 // FSAL type driver / steppers -----
50 #include "G4FSALIntegrationDriver.hh"
51 #include "G4VFSALIntegrationStepper.hh"
52 #include "G4RK547FEq1.hh"
53 // #include "G4RK547FEq2.hh"
54 // #include "G4RK547FEq3.hh"
55 // #include "G4FSALBogackiShampine45.hh"
56 // #include "G4FSALDormandPrince745.hh"
57
58 // Templated type drivers -----
59 #include "G4IntegrationDriver.hh"
60 #include "G4InterpolationDriver.hh"
61
62 #include "G4HelixHeum.hh"
63 #include "G4BFieldIntegrationDriver.hh"
64
65 #include "G4QSSDriverCreator.hh"
66
67 #include "G4CachedMagneticField.hh"
68
69 #include <cassert>
70 #include <memory>
71
72 G4bool G4ChordFinder::gVerboseCtor = false;
73 // ...........................................
74
75 G4ChordFinder::G4ChordFinder(G4VIntegrationDri
76 : fDefaultDeltaChord(0.25 * mm), fIntgrDrive
77 {
78 // Simple constructor -- it does not create
79 if( gVerboseCtor )
80 {
81 G4cout << "G4ChordFinder: Simple construct
82 }
83
84 fDeltaChord = fDefaultDeltaChord; // P
85 }
86
87 // ...........................................
88
89 G4ChordFinder::G4ChordFinder( G4MagneticField*
90 G4double
91 G4MagIntegratorS
92 G4int
93 : fDefaultDeltaChord(0.25 * mm)
94 {
95 // Construct the Chord Finder
96 // by creating in inverse order the Driver,
97 constexpr G4int nVar6 = 6; // Components i
98
99 fDeltaChord = fDefaultDeltaChord; // P
100
101 G4cout << " G4ChordFinder: stepperDriverId:
102
103 G4bool useFSALstepper = (stepperDriverId
104 G4bool useTemplatedStepper= (stepperDriverId
105 G4bool useRegularStepper = (stepperDriverId
106 G4bool useBfieldDriver = (stepperDriverId
107 G4bool useG4QSSDriver = (stepperDriverId
108
109 if( stepperDriverId == kQss3DriverType)
110 {
111 stepperDriverId = kQss2DriverType;
112 G4cout << " G4ChordFinder: QSS 3 is curren
113 }
114
115 using EquationType = G4Mag_UsualEqRhs;
116
117 using TemplatedStepperType =
118 G4TDormandPrince45<EquationType,nVar6
119 const char* TemplatedStepperName =
120 "G4TDormandPrince745 (templated Dormand-
121
122 using RegularStepperType =
123 G4DormandPrince745; // 5th order embe
124 // G4ClassicalRK4; // The old
125 // G4CashKarpRKF45; // First em
126 // G4BogackiShampine45; // High eff
127 // G4NystromRK4; // Nystrom
128 // G4RK547FEq1; // or 2 or 3
129 const char* RegularStepperName =
130 "G4DormandPrince745 (aka DOPRI5): 5th/4t
131 // "BogackiShampine 45 (Embedded 5th/4th
132 // "Nystrom stepper 4th order";
133
134 using NewFsalStepperType = G4DormandPrince74
135
136 const char* NewFSALStepperName =
137 "G4RK574FEq1> FSAL 4th/5th order 7-stage
138
139 #ifdef G4DEBUG_FIELD
140 static G4bool verboseDebug = true;
141 if( verboseDebug )
142 {
143 G4cout << "G4ChordFinder 2nd Constructor
144 G4cout << " Arguments: " << G4endl
145 << " - min step = " << stepMinimum
146 << " - stepper ptr provided : "
147 << ( pItsStepper==nullptr ? " no
148 if( pItsStepper==nullptr )
149 G4cout << " - stepper/driver Id = " <<
150 << " useFSAL = " << useFSALste
151 << " , useTemplated = " << use
152 << " , useRegular = " << useRe
153 << " , useFSAL = " << useFSALs
154 << G4endl;
155 }
156 #endif
157
158 // useHigherStepper = forceHigherEffiencySte
159
160 auto pEquation = new G4Mag_UsualEqRhs(theMa
161 fEquation = pEquation;
162
163 // G4MagIntegratorStepper* regularStepper =
164 // G4VFSALIntegrationStepper* fsalStepper =
165 // G4MagIntegratorStepper* oldFSALStepper =
166
167 G4bool errorInStepperCreation = false;
168
169 std::ostringstream message; // In case of f
170
171 if( pItsStepper != nullptr )
172 {
173 if( gVerboseCtor )
174 {
175 G4cout << " G4ChordFinder: Creating G4I
176 << " stepMinimum = " << stepMini
177 << " numVar= " << pItsStepper->G
178 }
179
180 // Stepper type is not known - so must us
181 if(pItsStepper->isQSS())
182 {
183 // fIntgrDriver = pItsStepper->build_
184 G4Exception("G4ChordFinder::G4ChordFi
185 "GeomField1001", FatalEx
186 "Cannot provide QSS ste
187 }
188 else
189 {
190 fIntgrDriver = new G4IntegrationDrive
191 pItsStepper,
192 // Stepper type is not known - so mus
193 // Non-interpolating driver used by d
194 // WAS: fIntgrDriver = pItsStepper->
195 }
196 // -- Older:
197 // G4cout << " G4ChordFinder: Creating G4
198 // Type is not known - so must use old cl
199 // fIntgrDriver = new G4MagInt_Driver( st
200 // pItsSt
201 }
202 else if ( useTemplatedStepper )
203 {
204 if( gVerboseCtor )
205 {
206 G4cout << " G4ChordFinder: Creating Te
207 << TemplatedStepperName << G4en
208 }
209 // RegularStepperType* regularStepper = n
210 auto templatedStepper = new TemplatedStep
211 // *** ***************
212 //
213 // Alternative - for G4NystromRK4:
214 // = new G4NystromRK4(pEquation, 0.1*mm )
215 fRegularStepperOwned = templatedStepper;
216 if( templatedStepper == nullptr )
217 {
218 message << "Templated Stepper instanti
219 message << "G4ChordFinder: Attempted t
220 << TemplatedStepperName << " t
221 errorInStepperCreation = true;
222 }
223 else
224 {
225 fIntgrDriver = new G4IntegrationDriver
226 stepMinimum, templatedStepper, nVar
227 if( gVerboseCtor )
228 {
229 G4cout << " G4ChordFinder: Using G4
230 }
231 }
232
233 }
234 else if ( useRegularStepper ) // Plain st
235 {
236 auto regularStepper = new RegularStepperT
237 // *** ***************
238 fRegularStepperOwned = regularStepper;
239
240 if( gVerboseCtor )
241 {
242 G4cout << " G4ChordFinder: Creating Dr
243 }
244
245 if( regularStepper == nullptr )
246 {
247 message << "Regular Stepper instantiat
248 message << "G4ChordFinder: Attempted t
249 << RegularStepperName << " typ
250 errorInStepperCreation = true;
251 }
252 else
253 {
254 auto dp5= dynamic_cast<G4DormandPrince
255 if( dp5 != nullptr )
256 {
257 fIntgrDriver = new G4InterpolationD
258 stepMinimum,
259 if( gVerboseCtor )
260 {
261 G4cout << " Using InterpolationD
262 }
263 }
264 else
265 {
266 fIntgrDriver = new G4IntegrationDri
267 stepMinimum,
268 if( gVerboseCtor )
269 {
270 G4cout << " Using IntegrationDri
271 }
272 }
273 }
274 }
275 else if ( useBfieldDriver )
276 {
277 auto regularStepper = new G4DormandPrince
278 // *** ***************
279 //
280 fRegularStepperOwned = regularStepper;
281
282 {
283 using SmallStepDriver = G4Interpolatio
284 using LargeStepDriver = G4IntegrationD
285
286 fLongStepper = std::make_unique<G4Heli
287
288 fIntgrDriver = new G4BFieldIntegration
289 std::make_unique<SmallStepDriver>(st
290 regularStepper, regularStepper->
291 std::make_unique<LargeStepDriver>(st
292 fLongStepper.get(), regularStepp
293
294 if( fIntgrDriver == nullptr)
295 {
296 message << "Using G4BFieldIntegrati
297 << RegularStepperName << "
298 message << "Driver instantiation FA
299 G4Exception("G4ChordFinder::G4Chord
300 "GeomField1001", JustWa
301 }
302 }
303 }
304 else if( useG4QSSDriver )
305 {
306 if( stepperDriverId == kQss2DriverType )
307 {
308 auto qssStepper2 = G4QSSDriverCreator::
309 if( gVerboseCtor )
310 {
311 G4cout << "-- Created QSS-2 stepper"
312 }
313 fIntgrDriver = G4QSSDriverCreator::Crea
314 }
315 else
316 {
317 auto qssStepper3 = G4QSSDriverCreator::
318 if( gVerboseCtor )
319 {
320 G4cout << "-- Created QSS-3 stepper"
321 }
322 fIntgrDriver = G4QSSDriverCreator::Crea
323 }
324 if( gVerboseCtor )
325 {
326 G4cout << "-- G4ChordFinder: Using QSS
327 }
328 }
329 else
330 {
331 auto fsalStepper= new NewFsalStepperType
332 // *** ******************
333 fNewFSALStepperOwned = fsalStepper;
334
335 if( fsalStepper == nullptr )
336 {
337 message << "Stepper instantiation FAIL
338 message << "Attempted to instantiate "
339 << NewFSALStepperName << " typ
340 G4Exception("G4ChordFinder::G4ChordFin
341 "GeomField1001", JustWarni
342 errorInStepperCreation = true;
343 }
344 else
345 {
346 fIntgrDriver = new
347 G4FSALIntegrationDriver<NewFsalStep
348 fsal
349 // ==== Create the driver which k
350
351 if( fIntgrDriver == nullptr )
352 {
353 message << "Using G4FSALIntegration
354 << NewFSALStepperName << G4
355 message << "Integration Driver inst
356 G4Exception("G4ChordFinder::G4Chord
357 "GeomField1001", JustWa
358 }
359 }
360 }
361
362 // -- Main work is now done
363
364 // Now check that no error occured, and r
365
366 // To test failure to create driver
367 // delete fIntgrDriver;
368 // fIntgrDriver = nullptr;
369
370 // Detect and report Error conditions
371 //
372 if( errorInStepperCreation || (fIntgrDriver
373 {
374 std::ostringstream errmsg;
375
376 if( errorInStepperCreation )
377 {
378 errmsg << "ERROR> Failure to create S
379 << " -------------------
380 }
381 if (fIntgrDriver == nullptr )
382 {
383 errmsg << "ERROR> Failure to create I
384 << G4endl
385 << " -------------------
386 << G4endl;
387 }
388 const std::string BoolName[2]= { "False",
389 errmsg << " Configuration: (constructor
390 << " provided Stepper = " << pI
391 << " stepper/driver Id = " << step
392 << " useTemplated = " << BoolNam
393 << " useRegular = " << BoolName[
394 << " useFSAL = " << BoolName[use
395 << " using combo BField Driver =
396 BoolName[ ! (useFSALstepper
397 || useRegularSt
398 << G4endl;
399 errmsg << message.str();
400 errmsg << "Aborting.";
401 G4Exception("G4ChordFinder::G4ChordFinder
402 "GeomField0003", FatalExcepti
403 }
404
405 assert( ( pItsStepper != nullptr )
406 || ( fRegularStepperOwned != nullptr
407 || ( fNewFSALStepperOwned != nullptr
408 || useG4QSSDriver
409 );
410 assert( fIntgrDriver != nullptr );
411 }
412
413 // ...........................................
414
415 G4ChordFinder::~G4ChordFinder()
416 {
417 delete fEquation;
418 delete fRegularStepperOwned;
419 delete fNewFSALStepperOwned;
420 delete fCachedField;
421 delete fIntgrDriver;
422 }
423
424 // ...........................................
425
426 G4FieldTrack
427 G4ChordFinder::ApproxCurvePointS( const G4Fiel
428 const G4Fiel
429 const G4Fiel
430 const G4Thre
431 const G4Thre
432 const G4Thre
433 G4bool
434 {
435 // ApproxCurvePointS is 2nd implementation o
436 // Use Brent Algorithm (or InvParabolic) whe
437 // Given a starting curve point A (CurveA_Po
438 // (CurveB_PointVelocity), a point E which i
439 // and a point F which is on the curve (fir
440 // point S on the curve closer to point E.
441 // While advancing towards S utilise 'eps_st
442 // relative accuracy of each Step.
443
444 G4FieldTrack EndPoint(CurveA_PointVelocity);
445 if(!first) { EndPoint = ApproxCurveV; }
446
447 G4ThreeVector Point_A,Point_B;
448 Point_A=CurveA_PointVelocity.GetPosition();
449 Point_B=CurveB_PointVelocity.GetPosition();
450
451 G4double xa,xb,xc,ya,yb,yc;
452
453 // InverseParabolic. AF Intersects (First Pa
454
455 if(first)
456 {
457 xa=0.;
458 ya=(PointG-Point_A).mag();
459 xb=(Point_A-CurrentF_Point).mag();
460 yb=-(PointG-CurrentF_Point).mag();
461 xc=(Point_A-Point_B).mag();
462 yc=-(CurrentE_Point-Point_B).mag();
463 }
464 else
465 {
466 xa=0.;
467 ya=(Point_A-CurrentE_Point).mag();
468 xb=(Point_A-CurrentF_Point).mag();
469 yb=(PointG-CurrentF_Point).mag();
470 xc=(Point_A-Point_B).mag();
471 yc=-(Point_B-PointG).mag();
472 if(xb==0.)
473 {
474 EndPoint = ApproxCurvePointV(CurveA_Poin
475 CurrentE_Po
476 return EndPoint;
477 }
478 }
479
480 const G4double tolerance = 1.e-12;
481 if(std::abs(ya)<=tolerance||std::abs(yc)<=to
482 {
483 ; // What to do for the moment: return the
484 // then PropagatorInField will take care
485 }
486 else
487 {
488 G4double test_step = InvParabolic(xa,ya,xb
489 G4double curve;
490 if(first)
491 {
492 curve=std::abs(EndPoint.GetCurveLength()
493 -ApproxCurveV.GetCurveLeng
494 }
495 else
496 {
497 test_step = test_step - xb;
498 curve=std::abs(EndPoint.GetCurveLength()
499 -CurveB_PointVelocity.GetC
500 xb = (CurrentF_Point-Point_B).mag();
501 }
502
503 if(test_step<=0) { test_step=0.1*xb; }
504 if(test_step>=xb) { test_step=0.5*xb; }
505 if(test_step>=curve){ test_step=0.5*curve;
506
507 if(curve*(1.+eps_step)<xb) // Similar to R
508 { // G4VIntersect
509 test_step=0.5*curve;
510 }
511
512 fIntgrDriver->AccurateAdvance(EndPoint,tes
513
514 #ifdef G4DEBUG_FIELD
515 // Printing Brent and Linear Approximation
516 //
517 G4cout << "G4ChordFinder::ApproxCurvePoint
518 << test_step << " EndPoint = " <<
519
520 // Test Track
521 //
522 G4FieldTrack TestTrack( CurveA_PointVeloci
523 TestTrack = ApproxCurvePointV( CurveA_Poin
524 CurveB_Poin
525 CurrentE_Po
526 G4cout.precision(14);
527 G4cout << "G4ChordFinder::BrentApprox = "
528 G4cout << "G4ChordFinder::LinearApprox= "
529 #endif
530 }
531 return EndPoint;
532 }
533
534
535 // ...........................................
536
537 G4FieldTrack G4ChordFinder::
538 ApproxCurvePointV( const G4FieldTrack& CurveA_
539 const G4FieldTrack& CurveB_
540 const G4ThreeVector& Curren
541 G4double eps_step)
542 {
543 // If r=|AE|/|AB|, and s=true path lenght (A
544 // return the point that is r*s along the cu
545
546 G4FieldTrack Current_PointVelocity = Curve
547
548 G4ThreeVector CurveA_Point= CurveA_PointVel
549 G4ThreeVector CurveB_Point= CurveB_PointVel
550
551 G4ThreeVector ChordAB_Vector= CurveB_Point
552 G4ThreeVector ChordAE_Vector= CurrentE_Poin
553
554 G4double ABdist= ChordAB_Vector.mag();
555 G4double curve_length; // A curve length
556 G4double AE_fraction;
557
558 curve_length= CurveB_PointVelocity.GetCurveL
559 - CurveA_PointVelocity.GetCurveL
560
561 G4double integrationInaccuracyLimit= std::ma
562 if( curve_length < ABdist * (1. - integratio
563 {
564 #ifdef G4DEBUG_FIELD
565 G4cerr << " Warning in G4ChordFinder::Appr
566 << G4endl
567 << " The two points are further apa
568 << G4endl
569 << " Dist = " << ABdist
570 << " curve length = " << curve_leng
571 << " relativeDiff = " << (curve_len
572 << G4endl;
573 if( curve_length < ABdist * (1. - 10*eps_s
574 {
575 std::ostringstream message;
576 message << "Unphysical curve length." <<
577 << "The size of the above differ
578 << G4endl
579 << "Aborting.";
580 G4Exception("G4ChordFinder::ApproxCurveP
581 FatalException, message);
582 }
583 #endif
584 // Take default corrective action: adjust
585 // NOTE: this case only happens for relati
586 // curve_length = ABdist;
587 }
588
589 G4double new_st_length;
590
591 if ( ABdist > 0.0 )
592 {
593 AE_fraction = ChordAE_Vector.mag() / ABdi
594 }
595 else
596 {
597 AE_fraction = 0.5;
598 #ifdef G4DEBUG_FIELD
599 G4cout << "Warning in G4ChordFinder::Appr
600 << " A and B are the same point!"
601 << " Chord AB length = " << ChordA
602 << G4endl;
603 #endif
604 }
605
606 if( (AE_fraction> 1.0 + perMillion) || (AE_f
607 {
608 #ifdef G4DEBUG_FIELD
609 G4cerr << " G4ChordFinder::ApproxCurvePoin
610 << " Anomalous condition:AE > AB or
611 << " AE_fraction = " << AE_fract
612 << " Chord AE length = " << Chord
613 << " Chord AB length = " << ABdis
614 G4cerr << " OK if this condition occurs af
615 << G4endl << " Otherwise it is an e
616 #endif
617 // This course can now result if B has be
618 // without E being recomputed (1 July 99)
619 // In this case this is not a "real error
620 // and we cope with it by a default corre
621 //
622 AE_fraction = 0.5;
623 }
624
625 new_st_length = AE_fraction * curve_length;
626
627 if ( AE_fraction > 0.0 )
628 {
629 fIntgrDriver->AccurateAdvance(Current_Poi
630 new_st_leng
631 //
632 // In this case it does not matter if it
633 }
634
635 // If there was a memory of the step_length
636 // of the integration Step, this could be re
637
638 G4cout.precision(14);
639
640 return Current_PointVelocity;
641 }
642
643 // ...........................................
644
645 std::ostream& operator<<( std::ostream& os, co
646 {
647 // Dumping the state of G4ChordFinder
648 os << "State of G4ChordFinder : " << std::e
649 os << " delta_chord = " << cf.fDeltaCh
650 os << " Default d_c = " << cf.fDefault
651
652 os << " stats-verbose = " << cf.fStatsVe
653
654 return os;
655 }
656