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