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