| 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 // $Id: G4ChordFinder.cc 107508 2017-11-20 08:23:14Z gcosmo $
>> 28 //
>> 29 //
>> 30 // 25.02.97 - John Apostolakis - Design and implementation
29 // ------------------------------------------- 31 // -------------------------------------------------------------------
30 32
31 #include <iomanip> 33 #include <iomanip>
32 34
33 #include "G4ChordFinder.hh" 35 #include "G4ChordFinder.hh"
34 #include "G4SystemOfUnits.hh" 36 #include "G4SystemOfUnits.hh"
35 #include "G4MagneticField.hh" 37 #include "G4MagneticField.hh"
36 #include "G4Mag_UsualEqRhs.hh" 38 #include "G4Mag_UsualEqRhs.hh"
37 #include "G4MagIntegratorDriver.hh" 39 #include "G4MagIntegratorDriver.hh"
38 // #include "G4ClassicalRK4.hh" 40 // #include "G4ClassicalRK4.hh"
39 // #include "G4CashKarpRKF45.hh" 41 // #include "G4CashKarpRKF45.hh"
40 // #include "G4NystromRK4.hh" <<
41 // #include "G4BogackiShampine23.hh" 42 // #include "G4BogackiShampine23.hh"
42 // #include "G4BogackiShampine45.hh" 43 // #include "G4BogackiShampine45.hh"
43 <<
44 #include "G4DormandPrince745.hh" 44 #include "G4DormandPrince745.hh"
45 45
46 // New templated stepper(s) -- avoid virtual c << 46 // New FSAL type driver / steppers -----
47 #include "G4TDormandPrince45.hh" <<
48 <<
49 // FSAL type driver / steppers ----- <<
50 #include "G4FSALIntegrationDriver.hh" 47 #include "G4FSALIntegrationDriver.hh"
51 #include "G4VFSALIntegrationStepper.hh" 48 #include "G4VFSALIntegrationStepper.hh"
52 #include "G4RK547FEq1.hh" 49 #include "G4RK547FEq1.hh"
53 // #include "G4RK547FEq2.hh" 50 // #include "G4RK547FEq2.hh"
54 // #include "G4RK547FEq3.hh" 51 // #include "G4RK547FEq3.hh"
55 // #include "G4FSALBogackiShampine45.hh" << 52 // #include "G4NystromRK4.hh"
56 // #include "G4FSALDormandPrince745.hh" <<
57 53
58 // Templated type drivers ----- << 54 // New FSAL type driver / steppers -----
59 #include "G4IntegrationDriver.hh" 55 #include "G4IntegrationDriver.hh"
60 #include "G4InterpolationDriver.hh" << 56 #include "G4FSALBogackiShampine45.hh"
61 << 57 // #include "G4FSALDormandPrince745.hh"
62 #include "G4HelixHeum.hh" <<
63 #include "G4BFieldIntegrationDriver.hh" <<
64 <<
65 #include "G4QSSDriverCreator.hh" <<
66 <<
67 #include "G4CachedMagneticField.hh" <<
68 58
69 #include <cassert> <<
70 #include <memory> <<
71 59
72 G4bool G4ChordFinder::gVerboseCtor = false; <<
73 // ........................................... 60 // ..........................................................................
74 61
75 G4ChordFinder::G4ChordFinder(G4VIntegrationDri 62 G4ChordFinder::G4ChordFinder(G4VIntegrationDriver* pIntegrationDriver)
76 : fDefaultDeltaChord(0.25 * mm), fIntgrDrive << 63 : fDefaultDeltaChord( 0.25 * mm ), // Parameters
>> 64 fDeltaChord( fDefaultDeltaChord ), // Internal parameters
>> 65 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
>> 66 fMultipleRadius(15.0),
>> 67 fStatsVerbose(0),
>> 68 fRegularStepperOwned(nullptr), // Dependent objects
>> 69 fEquation(0),
>> 70 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0)
77 { 71 {
78 // Simple constructor -- it does not create 72 // Simple constructor -- it does not create equation
79 if( gVerboseCtor ) << 73 fIntgrDriver= pIntegrationDriver;
80 { << 74
81 G4cout << "G4ChordFinder: Simple construct << 75 fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
82 } << 76
83 << 77 SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
84 fDeltaChord = fDefaultDeltaChord; // P << 78 // check the values and set the other parameters
>> 79
>> 80 // G4cout << "G4ChordFinder 1st Constructor called - (driver given). " << G4endl;
85 } 81 }
86 82
>> 83
87 // ........................................... 84 // ..........................................................................
88 85
89 G4ChordFinder::G4ChordFinder( G4MagneticField* 86 G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField,
90 G4double << 87 G4double stepMinimum,
91 G4MagIntegratorS << 88 G4MagIntegratorStepper* pItsStepper, // nullptr is default
92 G4int << 89 // G4bool useHigherEfficiencyStepper, // false by default
93 : fDefaultDeltaChord(0.25 * mm) << 90 G4bool useFSALstepper ) // false by default
>> 91 : fDefaultDeltaChord( 0.25 * mm ), // Constants
>> 92 fDeltaChord( fDefaultDeltaChord ), // Parameters
>> 93 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
>> 94 fMultipleRadius(15.0),
>> 95 fStatsVerbose(0),
>> 96 // fRegularStepperOwned(nullptr), // Dependent objects
>> 97 fEquation(0),
>> 98 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0) // State - stats
94 { 99 {
95 // Construct the Chord Finder << 100 // Construct the Chord Finder
96 // by creating in inverse order the Driver, << 101 // 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 <<
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 102
115 using EquationType = G4Mag_UsualEqRhs; << 103 using NewFsalStepperType = G4RK547FEq1; // or 2 or 3
116 << 104 const char* NewFSALStepperName = "G4RK574FEq1> FSAL 4th/5th order 7-stage 'Equilibrium-type' #1.";
117 using TemplatedStepperType = << 105 // using OldFsalStepperType = G4FSALBogackiShampine45;
118 G4TDormandPrince45<EquationType,nVar6 << 106 // const char* OldFSALStepperName = "FSAL BogackiShampine 45 (Embedded 5th/4th Order, 7-stage)";
119 const char* TemplatedStepperName = << 107 // = G4FSALDormandPrince745; // = "FSAL Dormand Prince 745 stepper";
120 "G4TDormandPrince745 (templated Dormand- <<
121 <<
122 using RegularStepperType = 108 using RegularStepperType =
123 G4DormandPrince745; // 5th order embe << 109 G4DormandPrince745; // Famous DOPRI5 (MatLab) 5th order embedded method. High efficiency.
124 // G4ClassicalRK4; // The old 110 // G4ClassicalRK4; // The old default
125 // G4CashKarpRKF45; // First em 111 // G4CashKarpRKF45; // First embedded method in G4
126 // G4BogackiShampine45; // High eff 112 // G4BogackiShampine45; // High efficiency 5th order embedded method
127 // G4NystromRK4; // Nystrom << 113 // G4NystromRK4(pEquation, 0.1*millimeter ); // *clhep::millimeter );
128 // G4RK547FEq1; // or 2 or 3 114 // G4RK547FEq1; // or 2 or 3
129 const char* RegularStepperName = << 115 const char* RegularStepperName = "G4DormandPrince745 (aka DOPRI5): 5th/4th Order 7-stage embedded stepper";
130 "G4DormandPrince745 (aka DOPRI5): 5th/4t <<
131 // "BogackiShampine 45 (Embedded 5th/4th 116 // "BogackiShampine 45 (Embedded 5th/4th Order, 7-stage)";
132 // "Nystrom stepper 4th order"; <<
133 117
134 using NewFsalStepperType = G4DormandPrince74 << 118 // Configurable
135 << 119 G4bool forceFSALstepper= false; // Choice - true to enable !!
136 const char* NewFSALStepperName = << 120 // G4bool useNewFSALtype= true;
137 "G4RK574FEq1> FSAL 4th/5th order 7-stage << 121 // G4bool forceHigherEffiencyStepper = false;
138 << 122 G4bool report = false; // Report type of stepper used
139 #ifdef G4DEBUG_FIELD << 123
140 static G4bool verboseDebug = true; << 124 bool recallFSALflag = useFSALstepper;
141 if( verboseDebug ) << 125 useFSALstepper = forceFSALstepper || useFSALstepper;
142 { << 126
>> 127 if( report ) {
143 G4cout << "G4ChordFinder 2nd Constructor 128 G4cout << "G4ChordFinder 2nd Constructor called. " << G4endl;
144 G4cout << " Arguments: " << G4endl << 129 G4cout << " Parameters: " << G4endl;
145 << " - min step = " << stepMinimum << 130 G4cout << " useFSAL stepper= " << useFSALstepper
146 << " - stepper ptr provided : " << 131 << " (request = " << recallFSALflag
147 << ( pItsStepper==nullptr ? " no << 132 << " 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 } 133 }
156 #endif <<
157 134
158 // useHigherStepper = forceHigherEffiencySte 135 // useHigherStepper = forceHigherEffiencyStepper || useHigherStepper;
159 << 136
160 auto pEquation = new G4Mag_UsualEqRhs(theMa << 137 G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
161 fEquation = pEquation; 138 fEquation = pEquation;
162 << 139 fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
163 // G4MagIntegratorStepper* regularStepper = << 140 // G4FieldTrack ??
164 // G4VFSALIntegrationStepper* fsalStepper = << 141 SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
165 // G4MagIntegratorStepper* oldFSALStepper = << 142 // check the values and set the other parameters
>> 143
>> 144 // G4MagIntegratorStepper* regularStepper = nullptr;
>> 145 // G4VFSALIntegrationStepper* fsalSepper = nullptr; // for new-type FSAL steppers only
>> 146 // NewFsalStepperType* fsalStepper =nullptr;
>> 147 // G4MagIntegratorStepper* oldFSALStepper =nullptr;
166 148
167 G4bool errorInStepperCreation = false; 149 G4bool errorInStepperCreation = false;
168 150
169 std::ostringstream message; // In case of f 151 std::ostringstream message; // In case of failure, load with description !
>> 152 message << "G4ChordFinder 2nd Constructor called. " << G4endl;
170 153
171 if( pItsStepper != nullptr ) 154 if( pItsStepper != nullptr )
172 { 155 {
173 if( gVerboseCtor ) << 156 // Type is not known - so must use old class
174 { << 157 fIntgrDriver = new G4MagInt_Driver(stepMinimum, pItsStepper,
175 G4cout << " G4ChordFinder: Creating G4I << 158 pItsStepper->GetNumberOfVariables() );
176 << " stepMinimum = " << stepMini << 159 }
177 << " numVar= " << pItsStepper->G << 160 else if ( !useFSALstepper )
178 } << 161 {
179 << 162 // RegularStepperType* regularStepper =nullptr; // To check the exception
180 // Stepper type is not known - so must us << 163 auto regularStepper =new RegularStepperType(pEquation);
181 if(pItsStepper->isQSS()) << 164 // *** ******************
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 // 165 //
213 // Alternative - for G4NystromRK4: 166 // Alternative - for G4NystromRK4:
214 // = new G4NystromRK4(pEquation, 0.1*mm ) << 167 // = 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; 168 fRegularStepperOwned = regularStepper;
239 169
240 if( gVerboseCtor ) <<
241 { <<
242 G4cout << " G4ChordFinder: Creating Dr <<
243 } <<
244 <<
245 if( regularStepper == nullptr ) 170 if( regularStepper == nullptr )
246 { 171 {
247 message << "Regular Stepper instantiat << 172 message << " ERROR> 'Regular' RK Stepper instantiation FAILED." << G4endl;
248 message << "G4ChordFinder: Attempted t 173 message << "G4ChordFinder: Attempted to instantiate "
249 << RegularStepperName << " typ 174 << RegularStepperName << " type stepper " << G4endl;
250 errorInStepperCreation = true; 175 errorInStepperCreation = true;
251 } 176 }
252 else 177 else
253 { 178 {
254 auto dp5= dynamic_cast<G4DormandPrince << 179 fIntgrDriver =
255 if( dp5 != nullptr ) << 180 new G4MagInt_Driver(stepMinimum,
256 { << 181 regularStepper,
257 fIntgrDriver = new G4InterpolationD << 182 regularStepper->GetNumberOfVariables() );
258 stepMinimum, << 183 // ==== Create the old type of driver
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 184
286 fLongStepper = std::make_unique<G4Heli << 185 // Alternative:
287 << 186 // new G4IntegrationDriver<RegularStepperType>(stepMinimum,
288 fIntgrDriver = new G4BFieldIntegration << 187 // ==== Create the driver which knows the class type
289 std::make_unique<SmallStepDriver>(st <<
290 regularStepper, regularStepper-> <<
291 std::make_unique<LargeStepDriver>(st <<
292 fLongStepper.get(), regularStepp <<
293 188
294 if( fIntgrDriver == nullptr) << 189 if( (fIntgrDriver==nullptr) || report ) {
295 { << 190 message << "G4ChordFinder: Using G4IntegrationDriver with "
296 message << "Using G4BFieldIntegrati <<
297 << RegularStepperName << " 191 << RegularStepperName << " type stepper " << G4endl;
298 message << "Driver instantiation FA <<
299 G4Exception("G4ChordFinder::G4Chord <<
300 "GeomField1001", JustWa <<
301 } 192 }
302 } << 193 if(fIntgrDriver==nullptr) {
303 } << 194 message << " ERROR> 'Regular' RK Driver instantiation FAILED." << G4endl;
304 else if( useG4QSSDriver ) << 195 }
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 } 196 }
328 } 197 }
329 else 198 else
330 { 199 {
331 auto fsalStepper= new NewFsalStepperType 200 auto fsalStepper= new NewFsalStepperType(pEquation);
332 // *** ****************** << 201 // ******************
333 fNewFSALStepperOwned = fsalStepper; 202 fNewFSALStepperOwned = fsalStepper;
>> 203 // delete fsalStepper;
>> 204 // /*NewFsalStepperType* */ fsalStepper =nullptr; // To check the exception
334 205
335 if( fsalStepper == nullptr ) 206 if( fsalStepper == nullptr )
336 { 207 {
337 message << "Stepper instantiation FAIL << 208 message << " ERROR> 'FSAL' RK Stepper instantiation FAILED." << G4endl;
338 message << "Attempted to instantiate " << 209 message << "G4ChordFinder: Attempted to instantiate "
339 << NewFSALStepperName << " typ 210 << NewFSALStepperName << " type stepper " << G4endl;
340 G4Exception("G4ChordFinder::G4ChordFin <<
341 "GeomField1001", JustWarni <<
342 errorInStepperCreation = true; 211 errorInStepperCreation = true;
343 } 212 }
344 else 213 else
345 { 214 {
346 fIntgrDriver = new 215 fIntgrDriver = new
347 G4FSALIntegrationDriver<NewFsalStep << 216 G4FSALIntegrationDriver<NewFsalStepperType>(stepMinimum,
348 fsal << 217 fsalStepper,
>> 218 fsalStepper->GetNumberOfVariables() );
349 // ==== Create the driver which k 219 // ==== Create the driver which knows the class type
350 220
351 if( fIntgrDriver == nullptr ) << 221 if( (fIntgrDriver==nullptr) || report ) {
352 { << 222 message << "G4ChordFinder: Using G4FSALIntegrationDriver with stepper type: " << G4endl
353 message << "Using G4FSALIntegration << 223 << NewFSALStepperName << " (new-FSAL type stepper.) " << G4endl;
354 << NewFSALStepperName << G4 << 224 }
355 message << "Integration Driver inst << 225 if(fIntgrDriver==nullptr) {
356 G4Exception("G4ChordFinder::G4Chord << 226 message << " ERROR> FSAL Integration Driver instantiation FAILED." << G4endl;
357 "GeomField1001", JustWa <<
358 } 227 }
359 } 228 }
360 } 229 }
361 230
362 // -- Main work is now done 231 // -- Main work is now done
363 232
364 // Now check that no error occured, and r 233 // Now check that no error occured, and report it if one did.
365 234
366 // To test failure to create driver 235 // To test failure to create driver
367 // delete fIntgrDriver; 236 // delete fIntgrDriver;
368 // fIntgrDriver = nullptr; << 237 // fIntgrDriver= nullptr;
369 238
370 // Detect and report Error conditions 239 // Detect and report Error conditions
371 // <<
372 if( errorInStepperCreation || (fIntgrDriver 240 if( errorInStepperCreation || (fIntgrDriver == nullptr ))
373 { 241 {
374 std::ostringstream errmsg; 242 std::ostringstream errmsg;
375 243
376 if( errorInStepperCreation ) 244 if( errorInStepperCreation )
377 { 245 {
378 errmsg << "ERROR> Failure to create S 246 errmsg << "ERROR> Failure to create Stepper object." << G4endl
379 << " ------------------- 247 << " --------------------------------" << G4endl;
380 } 248 }
381 if (fIntgrDriver == nullptr ) 249 if (fIntgrDriver == nullptr )
382 { 250 {
383 errmsg << "ERROR> Failure to create I << 251 errmsg << "ERROR> Failure to create Integration-Driver object." << G4endl
384 << G4endl << 252 << " -------------------------------------------" << G4endl;
385 << " ------------------- <<
386 << G4endl; <<
387 } 253 }
388 const std::string BoolName[2]= { "False", 254 const std::string BoolName[2]= { "False", "True" };
389 errmsg << " Configuration: (constructor 255 errmsg << " Configuration: (constructor arguments) " << G4endl
390 << " provided Stepper = " << pI 256 << " provided Stepper = " << pItsStepper << G4endl
391 << " stepper/driver Id = " << step << 257 << " use FSAL stepper = " << BoolName[useFSALstepper]
392 << " useTemplated = " << BoolNam << 258 // ( useFSALstepper ? "True" : "False" )
393 << " useRegular = " << BoolName[ << 259 << " (request = " << BoolName[recallFSALflag]
394 << " useFSAL = " << BoolName[use << 260 << " force FSAL = " << BoolName[forceFSALstepper] << " )" << G4endl;
395 << " using combo BField Driver = << 261 // << " use new FSAL stp = " << ( useNewFSALstepper ? "True" : "False" ) << G4endl;
396 BoolName[ ! (useFSALstepper <<
397 || useRegularSt <<
398 << G4endl; <<
399 errmsg << message.str(); 262 errmsg << message.str();
400 errmsg << "Aborting."; 263 errmsg << "Aborting.";
401 G4Exception("G4ChordFinder::G4ChordFinder << 264 G4Exception("G4ChordFinder::G4ChordFinder() - constructor 2", "GeomField0003",
402 "GeomField0003", FatalExcepti << 265 FatalException, errmsg);
>> 266 }
>> 267 else if ( report )
>> 268 {
>> 269 G4cout << message.str();
403 } 270 }
404 271
405 assert( ( pItsStepper != nullptr ) << 272 assert( ( pItsStepper != nullptr )
406 || ( fRegularStepperOwned != nullptr 273 || ( fRegularStepperOwned != nullptr )
407 || ( fNewFSALStepperOwned != nullptr << 274 || ( fNewFSALStepperOwned != nullptr )
408 || useG4QSSDriver << 275 // || ( fOldFSALStepperOwned != nullptr )
409 ); 276 );
410 assert( fIntgrDriver != nullptr ); 277 assert( fIntgrDriver != nullptr );
>> 278
411 } 279 }
412 280
>> 281
413 // ........................................... 282 // ......................................................................
414 283
415 G4ChordFinder::~G4ChordFinder() 284 G4ChordFinder::~G4ChordFinder()
416 { 285 {
417 delete fEquation; << 286 delete fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
418 delete fRegularStepperOwned; << 287 delete fRegularStepperOwned;
419 delete fNewFSALStepperOwned; << 288 delete fNewFSALStepperOwned;
420 delete fCachedField; << 289 // delete fOldFSALStepperOwned;
421 delete fIntgrDriver; << 290 delete fIntgrDriver;
>> 291
>> 292 if( fStatsVerbose ) { PrintStatistics(); }
>> 293 }
>> 294
>> 295
>> 296 // ......................................................................
>> 297
>> 298 void
>> 299 G4ChordFinder::SetFractions_Last_Next( G4double fractLast, G4double fractNext )
>> 300 {
>> 301 // Use -1.0 as request for Default.
>> 302 if( fractLast == -1.0 ) fractLast = 1.0; // 0.9;
>> 303 if( fractNext == -1.0 ) fractNext = 0.98; // 0.9;
>> 304
>> 305 // fFirstFraction = 0.999; // Orig 0.999 A safe value, range: ~ 0.95 - 0.999
>> 306 // fMultipleRadius = 15.0; // For later use, range: ~ 2 - 20
>> 307
>> 308 if( fStatsVerbose )
>> 309 {
>> 310 G4cout << " ChordFnd> Trying to set fractions: "
>> 311 << " first " << fFirstFraction
>> 312 << " last " << fractLast
>> 313 << " next " << fractNext
>> 314 << " and multiple " << fMultipleRadius
>> 315 << G4endl;
>> 316 }
>> 317
>> 318 if( (fractLast > 0.0) && (fractLast <=1.0) )
>> 319 {
>> 320 fFractionLast= fractLast;
>> 321 }
>> 322 else
>> 323 {
>> 324 G4cerr << "G4ChordFinder::SetFractions_Last_Next: Invalid "
>> 325 << " fraction Last = " << fractLast
>> 326 << " must be 0 < fractionLast <= 1 " << G4endl;
>> 327 }
>> 328 if( (fractNext > 0.0) && (fractNext <1.0) )
>> 329 {
>> 330 fFractionNextEstimate = fractNext;
>> 331 }
>> 332 else
>> 333 {
>> 334 G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
>> 335 << " fraction Next = " << fractNext
>> 336 << " must be 0 < fractionNext < 1 " << G4endl;
>> 337 }
422 } 338 }
423 339
>> 340
>> 341 // ......................................................................
>> 342
>> 343 G4double
>> 344 G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
>> 345 G4double stepMax,
>> 346 G4double epsStep,
>> 347 const G4ThreeVector& latestSafetyOrigin,
>> 348 G4double latestSafetyRadius )
>> 349 {
>> 350 G4double stepPossible;
>> 351 G4double dyErr;
>> 352 G4FieldTrack yEnd( yCurrent);
>> 353 G4double startCurveLen= yCurrent.GetCurveLength();
>> 354 G4double nextStep;
>> 355 // *************
>> 356 stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep,
>> 357 &nextStep, latestSafetyOrigin, latestSafetyRadius
>> 358 );
>> 359 // *************
>> 360
>> 361 G4bool good_advance;
>> 362
>> 363 if ( dyErr < epsStep * stepPossible )
>> 364 {
>> 365 // Accept this accuracy.
>> 366
>> 367 yCurrent = yEnd;
>> 368 good_advance = true;
>> 369 }
>> 370 else
>> 371 {
>> 372 // Advance more accurately to "end of chord"
>> 373 // ***************
>> 374 good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible,
>> 375 epsStep, nextStep);
>> 376 if ( ! good_advance )
>> 377 {
>> 378 // In this case the driver could not do the full distance
>> 379 stepPossible= yCurrent.GetCurveLength()-startCurveLen;
>> 380 }
>> 381 }
>> 382 return stepPossible;
>> 383 }
>> 384
>> 385
>> 386 // ............................................................................
>> 387
>> 388 G4double
>> 389 G4ChordFinder::FindNextChord( const G4FieldTrack& yStart,
>> 390 G4double stepMax,
>> 391 G4FieldTrack& yEnd, // Endpoint
>> 392 G4double& dyErrPos, // Error of endpoint
>> 393 G4double epsStep,
>> 394 G4double* pStepForAccuracy,
>> 395 const G4ThreeVector, // latestSafetyOrigin,
>> 396 G4double // latestSafetyRadius
>> 397 )
>> 398 {
>> 399 // Returns Length of Step taken
>> 400
>> 401 // G4cout << ">G4ChordFinder::FindNextChord called." << G4endl;
>> 402
>> 403 G4FieldTrack yCurrent= yStart;
>> 404 G4double stepTrial, stepForAccuracy;
>> 405 G4double dydx[G4FieldTrack::ncompSVEC];
>> 406
>> 407 // 1.) Try to "leap" to end of interval
>> 408 // 2.) Evaluate if resulting chord gives d_chord that is good enough.
>> 409 // 2a.) If d_chord is not good enough, find one that is.
>> 410
>> 411 G4bool validEndPoint= false;
>> 412 G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
>> 413
>> 414 fIntgrDriver-> GetDerivatives( yCurrent, dydx );
>> 415
>> 416 unsigned int noTrials=0;
>> 417 const unsigned int maxTrials= 75; // Avoid endless loop for bad convergence
>> 418
>> 419 const G4double safetyFactor= fFirstFraction; // 0.975 or 0.99 ? was 0.999
>> 420
>> 421 stepTrial = std::min( stepMax, safetyFactor*fLastStepEstimate_Unconstrained );
>> 422
>> 423 G4double newStepEst_Uncons= 0.0;
>> 424 G4double stepForChord;
>> 425 do
>> 426 {
>> 427 yCurrent = yStart; // Always start from initial point
>> 428 // ************
>> 429 fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
>> 430 dChordStep, dyErrPos);
>> 431 // ************
>> 432
>> 433 // We check whether the criterion is met here.
>> 434 validEndPoint = AcceptableMissDist(dChordStep);
>> 435
>> 436 lastStepLength = stepTrial;
>> 437
>> 438 // This method estimates to step size for a good chord.
>> 439 stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
>> 440
>> 441 if( ! validEndPoint )
>> 442 {
>> 443 if( stepTrial<=0.0 )
>> 444 {
>> 445 stepTrial = stepForChord;
>> 446 }
>> 447 else if (stepForChord <= stepTrial)
>> 448 {
>> 449 // Reduce by a fraction, possibly up to 20%
>> 450 stepTrial = std::min( stepForChord, fFractionLast * stepTrial);
>> 451 }
>> 452 else
>> 453 {
>> 454 stepTrial *= 0.1;
>> 455 }
>> 456 }
>> 457 noTrials++;
>> 458 }
>> 459 while( (! validEndPoint) && (noTrials < maxTrials) );
>> 460 // Loop checking, 07.10.2016, J. Apostolakis
>> 461
>> 462 if( noTrials >= maxTrials )
>> 463 {
>> 464 std::ostringstream message;
>> 465 message << "Exceeded maximum number of trials= " << maxTrials << G4endl
>> 466 << "Current sagita dist= " << dChordStep << G4endl
>> 467 << "Step sizes (actual and proposed): " << G4endl
>> 468 << "Last trial = " << lastStepLength << G4endl
>> 469 << "Next trial = " << stepTrial << G4endl
>> 470 << "Proposed for chord = " << stepForChord << G4endl
>> 471 ;
>> 472 G4Exception("G4ChordFinder::FindNextChord()", "GeomField0003",
>> 473 JustWarning, message);
>> 474 }
>> 475
>> 476 if( newStepEst_Uncons > 0.0 )
>> 477 {
>> 478 fLastStepEstimate_Unconstrained= newStepEst_Uncons;
>> 479 }
>> 480
>> 481 AccumulateStatistics( noTrials );
>> 482
>> 483 if( pStepForAccuracy )
>> 484 {
>> 485 // Calculate the step size required for accuracy, if it is needed
>> 486 //
>> 487 G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
>> 488 if( dyErr_relative > 1.0 )
>> 489 {
>> 490 stepForAccuracy = fIntgrDriver->ComputeNewStepSize( dyErr_relative,
>> 491 lastStepLength );
>> 492 }
>> 493 else
>> 494 {
>> 495 stepForAccuracy = 0.0; // Convention to show step was ok
>> 496 }
>> 497 *pStepForAccuracy = stepForAccuracy;
>> 498 }
>> 499
>> 500 #ifdef TEST_CHORD_PRINT
>> 501 static int dbg=0;
>> 502 if( dbg )
>> 503 {
>> 504 G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
>> 505 << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
>> 506 }
>> 507 #endif
>> 508 yEnd= yCurrent;
>> 509 return stepTrial;
>> 510 }
>> 511
>> 512
>> 513 // ...........................................................................
>> 514
>> 515 G4double G4ChordFinder::NewStep(G4double stepTrialOld,
>> 516 G4double dChordStep, // Curr. dchord achieved
>> 517 G4double& stepEstimate_Unconstrained )
>> 518 {
>> 519 // Is called to estimate the next step size, even for successful steps,
>> 520 // in order to predict an accurate 'chord-sensitive' first step
>> 521 // which is likely to assist in more performant 'stepping'.
>> 522
>> 523 G4double stepTrial;
>> 524
>> 525 #if 1
>> 526
>> 527 if (dChordStep > 0.0)
>> 528 {
>> 529 stepEstimate_Unconstrained =
>> 530 stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
>> 531 stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
>> 532 }
>> 533 else
>> 534 {
>> 535 // Should not update the Unconstrained Step estimate: incorrect!
>> 536 stepTrial = stepTrialOld * 2.;
>> 537 }
>> 538
>> 539 if( stepTrial <= 0.001 * stepTrialOld)
>> 540 {
>> 541 if ( dChordStep > 1000.0 * fDeltaChord )
>> 542 {
>> 543 stepTrial= stepTrialOld * 0.03;
>> 544 }
>> 545 else
>> 546 {
>> 547 if ( dChordStep > 100. * fDeltaChord )
>> 548 {
>> 549 stepTrial= stepTrialOld * 0.1;
>> 550 }
>> 551 else // Try halving the length until dChordStep OK
>> 552 {
>> 553 stepTrial= stepTrialOld * 0.5;
>> 554 }
>> 555 }
>> 556 }
>> 557 else if (stepTrial > 1000.0 * stepTrialOld)
>> 558 {
>> 559 stepTrial= 1000.0 * stepTrialOld;
>> 560 }
>> 561
>> 562 if( stepTrial == 0.0 )
>> 563 {
>> 564 stepTrial= 0.000001;
>> 565 }
>> 566
>> 567 #else
>> 568
>> 569 if ( dChordStep > 1000. * fDeltaChord )
>> 570 {
>> 571 stepTrial= stepTrialOld * 0.03;
>> 572 }
>> 573 else
>> 574 {
>> 575 if ( dChordStep > 100. * fDeltaChord )
>> 576 {
>> 577 stepTrial= stepTrialOld * 0.1;
>> 578 }
>> 579 else // Keep halving the length until dChordStep OK
>> 580 {
>> 581 stepTrial= stepTrialOld * 0.5;
>> 582 }
>> 583 }
>> 584
>> 585 #endif
>> 586
>> 587 // A more sophisticated chord-finder could figure out a better
>> 588 // stepTrial, from dChordStep and the required d_geometry
>> 589 // e.g.
>> 590 // Calculate R, r_helix (eg at orig point)
>> 591 // if( stepTrial < 2 pi R )
>> 592 // stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
>> 593 // else
>> 594 // ??
>> 595
>> 596 return stepTrial;
>> 597 }
>> 598
>> 599
424 // ........................................... 600 // ...........................................................................
425 601
426 G4FieldTrack 602 G4FieldTrack
427 G4ChordFinder::ApproxCurvePointS( const G4Fiel 603 G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
428 const G4Fiel 604 const G4FieldTrack& CurveB_PointVelocity,
429 const G4Fiel 605 const G4FieldTrack& ApproxCurveV,
430 const G4Thre 606 const G4ThreeVector& CurrentE_Point,
431 const G4Thre 607 const G4ThreeVector& CurrentF_Point,
432 const G4Thre 608 const G4ThreeVector& PointG,
433 G4bool << 609 G4bool first, G4double eps_step)
434 { 610 {
435 // ApproxCurvePointS is 2nd implementation o 611 // ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
436 // Use Brent Algorithm (or InvParabolic) whe 612 // Use Brent Algorithm (or InvParabolic) when possible.
437 // Given a starting curve point A (CurveA_Po 613 // Given a starting curve point A (CurveA_PointVelocity), curve point B
438 // (CurveB_PointVelocity), a point E which i 614 // (CurveB_PointVelocity), a point E which is (generally) not on the curve
439 // and a point F which is on the curve (fir 615 // and a point F which is on the curve (first approximation), find new
440 // point S on the curve closer to point E. 616 // point S on the curve closer to point E.
441 // While advancing towards S utilise 'eps_st 617 // While advancing towards S utilise 'eps_step' as a measure of the
442 // relative accuracy of each Step. 618 // relative accuracy of each Step.
443 619
444 G4FieldTrack EndPoint(CurveA_PointVelocity); 620 G4FieldTrack EndPoint(CurveA_PointVelocity);
445 if(!first) { EndPoint = ApproxCurveV; } << 621 if(!first){EndPoint= ApproxCurveV;}
446 622
447 G4ThreeVector Point_A,Point_B; 623 G4ThreeVector Point_A,Point_B;
448 Point_A=CurveA_PointVelocity.GetPosition(); 624 Point_A=CurveA_PointVelocity.GetPosition();
449 Point_B=CurveB_PointVelocity.GetPosition(); 625 Point_B=CurveB_PointVelocity.GetPosition();
450 626
451 G4double xa,xb,xc,ya,yb,yc; 627 G4double xa,xb,xc,ya,yb,yc;
452 628
453 // InverseParabolic. AF Intersects (First Pa 629 // InverseParabolic. AF Intersects (First Part of Curve)
454 630
455 if(first) 631 if(first)
456 { 632 {
457 xa=0.; 633 xa=0.;
458 ya=(PointG-Point_A).mag(); 634 ya=(PointG-Point_A).mag();
459 xb=(Point_A-CurrentF_Point).mag(); 635 xb=(Point_A-CurrentF_Point).mag();
460 yb=-(PointG-CurrentF_Point).mag(); 636 yb=-(PointG-CurrentF_Point).mag();
461 xc=(Point_A-Point_B).mag(); 637 xc=(Point_A-Point_B).mag();
462 yc=-(CurrentE_Point-Point_B).mag(); 638 yc=-(CurrentE_Point-Point_B).mag();
463 } 639 }
464 else 640 else
465 { 641 {
466 xa=0.; 642 xa=0.;
467 ya=(Point_A-CurrentE_Point).mag(); 643 ya=(Point_A-CurrentE_Point).mag();
468 xb=(Point_A-CurrentF_Point).mag(); 644 xb=(Point_A-CurrentF_Point).mag();
469 yb=(PointG-CurrentF_Point).mag(); 645 yb=(PointG-CurrentF_Point).mag();
470 xc=(Point_A-Point_B).mag(); 646 xc=(Point_A-Point_B).mag();
471 yc=-(Point_B-PointG).mag(); 647 yc=-(Point_B-PointG).mag();
472 if(xb==0.) 648 if(xb==0.)
473 { 649 {
474 EndPoint = ApproxCurvePointV(CurveA_Poin << 650 EndPoint=
475 CurrentE_Po << 651 ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity,
>> 652 CurrentE_Point, eps_step);
476 return EndPoint; 653 return EndPoint;
477 } 654 }
478 } 655 }
479 656
480 const G4double tolerance = 1.e-12; << 657 const G4double tolerance= 1.e-12;
481 if(std::abs(ya)<=tolerance||std::abs(yc)<=to 658 if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
482 { 659 {
483 ; // What to do for the moment: return the 660 ; // What to do for the moment: return the same point as at start
484 // then PropagatorInField will take care 661 // then PropagatorInField will take care
485 } 662 }
486 else 663 else
487 { 664 {
488 G4double test_step = InvParabolic(xa,ya,xb 665 G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
489 G4double curve; 666 G4double curve;
490 if(first) 667 if(first)
491 { 668 {
492 curve=std::abs(EndPoint.GetCurveLength() 669 curve=std::abs(EndPoint.GetCurveLength()
493 -ApproxCurveV.GetCurveLeng 670 -ApproxCurveV.GetCurveLength());
494 } 671 }
495 else 672 else
496 { 673 {
497 test_step = test_step - xb; << 674 test_step=(test_step-xb);
498 curve=std::abs(EndPoint.GetCurveLength() 675 curve=std::abs(EndPoint.GetCurveLength()
499 -CurveB_PointVelocity.GetC 676 -CurveB_PointVelocity.GetCurveLength());
500 xb = (CurrentF_Point-Point_B).mag(); << 677 xb=(CurrentF_Point-Point_B).mag();
501 } 678 }
502 679
503 if(test_step<=0) { test_step=0.1*xb; } 680 if(test_step<=0) { test_step=0.1*xb; }
504 if(test_step>=xb) { test_step=0.5*xb; } 681 if(test_step>=xb) { test_step=0.5*xb; }
505 if(test_step>=curve){ test_step=0.5*curve; 682 if(test_step>=curve){ test_step=0.5*curve; }
506 683
507 if(curve*(1.+eps_step)<xb) // Similar to R 684 if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from
508 { // G4VIntersect 685 { // G4VIntersectionLocator
509 test_step=0.5*curve; 686 test_step=0.5*curve;
510 } 687 }
511 688
512 fIntgrDriver->AccurateAdvance(EndPoint,tes 689 fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
513 690
514 #ifdef G4DEBUG_FIELD 691 #ifdef G4DEBUG_FIELD
515 // Printing Brent and Linear Approximation 692 // Printing Brent and Linear Approximation
516 // 693 //
517 G4cout << "G4ChordFinder::ApproxCurvePoint 694 G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
518 << test_step << " EndPoint = " << 695 << test_step << " EndPoint = " << EndPoint << G4endl;
519 696
520 // Test Track 697 // Test Track
521 // 698 //
522 G4FieldTrack TestTrack( CurveA_PointVeloci 699 G4FieldTrack TestTrack( CurveA_PointVelocity);
523 TestTrack = ApproxCurvePointV( CurveA_Poin 700 TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
524 CurveB_Poin 701 CurveB_PointVelocity,
525 CurrentE_Po 702 CurrentE_Point, eps_step );
526 G4cout.precision(14); 703 G4cout.precision(14);
527 G4cout << "G4ChordFinder::BrentApprox = " 704 G4cout << "G4ChordFinder::BrentApprox = " << EndPoint << G4endl;
528 G4cout << "G4ChordFinder::LinearApprox= " 705 G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl;
529 #endif 706 #endif
530 } 707 }
531 return EndPoint; 708 return EndPoint;
532 } 709 }
533 710
534 711
535 // ........................................... 712 // ...........................................................................
536 713
537 G4FieldTrack G4ChordFinder:: 714 G4FieldTrack G4ChordFinder::
538 ApproxCurvePointV( const G4FieldTrack& CurveA_ 715 ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
539 const G4FieldTrack& CurveB_ 716 const G4FieldTrack& CurveB_PointVelocity,
540 const G4ThreeVector& Curren 717 const G4ThreeVector& CurrentE_Point,
541 G4double eps_step) 718 G4double eps_step)
542 { 719 {
543 // If r=|AE|/|AB|, and s=true path lenght (A 720 // If r=|AE|/|AB|, and s=true path lenght (AB)
544 // return the point that is r*s along the cu 721 // return the point that is r*s along the curve!
545 722
546 G4FieldTrack Current_PointVelocity = Curve 723 G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
547 724
548 G4ThreeVector CurveA_Point= CurveA_PointVel 725 G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
549 G4ThreeVector CurveB_Point= CurveB_PointVel 726 G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
550 727
551 G4ThreeVector ChordAB_Vector= CurveB_Point 728 G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
552 G4ThreeVector ChordAE_Vector= CurrentE_Poin 729 G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
553 730
554 G4double ABdist= ChordAB_Vector.mag(); 731 G4double ABdist= ChordAB_Vector.mag();
555 G4double curve_length; // A curve length 732 G4double curve_length; // A curve length of AB
556 G4double AE_fraction; 733 G4double AE_fraction;
557 734
558 curve_length= CurveB_PointVelocity.GetCurveL 735 curve_length= CurveB_PointVelocity.GetCurveLength()
559 - CurveA_PointVelocity.GetCurveL 736 - CurveA_PointVelocity.GetCurveLength();
560 737
561 G4double integrationInaccuracyLimit= std::ma << 738 G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
562 if( curve_length < ABdist * (1. - integratio 739 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
563 { 740 {
564 #ifdef G4DEBUG_FIELD 741 #ifdef G4DEBUG_FIELD
565 G4cerr << " Warning in G4ChordFinder::Appr 742 G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
566 << G4endl 743 << G4endl
567 << " The two points are further apa 744 << " The two points are further apart than the curve length "
568 << G4endl 745 << G4endl
569 << " Dist = " << ABdist 746 << " Dist = " << ABdist
570 << " curve length = " << curve_leng 747 << " curve length = " << curve_length
571 << " relativeDiff = " << (curve_len 748 << " relativeDiff = " << (curve_length-ABdist)/ABdist
572 << G4endl; 749 << G4endl;
573 if( curve_length < ABdist * (1. - 10*eps_s 750 if( curve_length < ABdist * (1. - 10*eps_step) )
574 { 751 {
575 std::ostringstream message; 752 std::ostringstream message;
576 message << "Unphysical curve length." << 753 message << "Unphysical curve length." << G4endl
577 << "The size of the above differ 754 << "The size of the above difference exceeds allowed limits."
578 << G4endl 755 << G4endl
579 << "Aborting."; 756 << "Aborting.";
580 G4Exception("G4ChordFinder::ApproxCurveP 757 G4Exception("G4ChordFinder::ApproxCurvePointV()", "GeomField0003",
581 FatalException, message); 758 FatalException, message);
582 } 759 }
583 #endif 760 #endif
584 // Take default corrective action: adjust 761 // Take default corrective action: adjust the maximum curve length.
585 // NOTE: this case only happens for relati 762 // NOTE: this case only happens for relatively straight paths.
586 // curve_length = ABdist; 763 // curve_length = ABdist;
587 } 764 }
588 765
589 G4double new_st_length; << 766 G4double new_st_length;
590 767
591 if ( ABdist > 0.0 ) 768 if ( ABdist > 0.0 )
592 { 769 {
593 AE_fraction = ChordAE_Vector.mag() / ABdi 770 AE_fraction = ChordAE_Vector.mag() / ABdist;
594 } 771 }
595 else 772 else
596 { 773 {
597 AE_fraction = 0.5; 774 AE_fraction = 0.5; // Guess .. ?;
598 #ifdef G4DEBUG_FIELD 775 #ifdef G4DEBUG_FIELD
599 G4cout << "Warning in G4ChordFinder::Appr 776 G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
600 << " A and B are the same point!" 777 << " A and B are the same point!" << G4endl
601 << " Chord AB length = " << ChordA 778 << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
602 << G4endl; 779 << G4endl;
603 #endif 780 #endif
604 } 781 }
605 782
606 if( (AE_fraction> 1.0 + perMillion) || (AE_f 783 if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) )
607 { 784 {
608 #ifdef G4DEBUG_FIELD 785 #ifdef G4DEBUG_FIELD
609 G4cerr << " G4ChordFinder::ApproxCurvePoin 786 G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
610 << " Anomalous condition:AE > AB or 787 << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
611 << " AE_fraction = " << AE_fract 788 << " AE_fraction = " << AE_fraction << G4endl
612 << " Chord AE length = " << Chord 789 << " Chord AE length = " << ChordAE_Vector.mag() << G4endl
613 << " Chord AB length = " << ABdis 790 << " Chord AB length = " << ABdist << G4endl << G4endl;
614 G4cerr << " OK if this condition occurs af 791 G4cerr << " OK if this condition occurs after a recalculation of 'B'"
615 << G4endl << " Otherwise it is an e 792 << G4endl << " Otherwise it is an error. " << G4endl ;
616 #endif 793 #endif
617 // This course can now result if B has be 794 // This course can now result if B has been re-evaluated,
618 // without E being recomputed (1 July 99) 795 // without E being recomputed (1 July 99).
619 // In this case this is not a "real error 796 // In this case this is not a "real error" - but it is undesired
620 // and we cope with it by a default corre 797 // and we cope with it by a default corrective action ...
621 // 798 //
622 AE_fraction = 0.5; 799 AE_fraction = 0.5; // Default value
623 } 800 }
624 801
625 new_st_length = AE_fraction * curve_length; << 802 new_st_length= AE_fraction * curve_length;
626 803
627 if ( AE_fraction > 0.0 ) 804 if ( AE_fraction > 0.0 )
628 { 805 {
629 fIntgrDriver->AccurateAdvance(Current_Poi 806 fIntgrDriver->AccurateAdvance(Current_PointVelocity,
630 new_st_leng 807 new_st_length, eps_step );
631 // 808 //
632 // In this case it does not matter if it 809 // In this case it does not matter if it cannot advance the full distance
633 } 810 }
634 811
635 // If there was a memory of the step_length 812 // If there was a memory of the step_length actually required at the start
636 // of the integration Step, this could be re 813 // of the integration Step, this could be re-used ...
637 814
638 G4cout.precision(14); 815 G4cout.precision(14);
639 816
640 return Current_PointVelocity; 817 return Current_PointVelocity;
641 } 818 }
642 819
643 // ........................................... <<
644 820
645 std::ostream& operator<<( std::ostream& os, co << 821 // ......................................................................
>> 822
>> 823 void
>> 824 G4ChordFinder::PrintStatistics()
646 { 825 {
647 // Dumping the state of G4ChordFinder << 826 // Print Statistics
648 os << "State of G4ChordFinder : " << std::e << 827
649 os << " delta_chord = " << cf.fDeltaCh << 828 G4cout << "G4ChordFinder statistics report: " << G4endl;
650 os << " Default d_c = " << cf.fDefault << 829 G4cout
>> 830 << " No trials: " << fTotalNoTrials_FNC
>> 831 << " No Calls: " << fNoCalls_FNC
>> 832 << " Max-trial: " << fmaxTrials_FNC
>> 833 << G4endl;
>> 834 G4cout
>> 835 << " Parameters: "
>> 836 << " fFirstFraction " << fFirstFraction
>> 837 << " fFractionLast " << fFractionLast
>> 838 << " fFractionNextEstimate " << fFractionNextEstimate
>> 839 << G4endl;
>> 840 }
>> 841
651 842
652 os << " stats-verbose = " << cf.fStatsVe << 843 // ...........................................................................
653 844
654 return os; << 845 void G4ChordFinder::TestChordPrint( G4int noTrials,
>> 846 G4int lastStepTrial,
>> 847 G4double dChordStep,
>> 848 G4double nextStepTrial )
>> 849 {
>> 850 G4int oldprec= G4cout.precision(5);
>> 851 G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
>> 852 << " this_step= " << std::setw(10) << lastStepTrial;
>> 853 if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
>> 854 {
>> 855 G4cout.precision(8);
>> 856 }
>> 857 else
>> 858 {
>> 859 G4cout.precision(6);
>> 860 }
>> 861 G4cout << " dChordStep= " << std::setw(12) << dChordStep;
>> 862 if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
>> 863 else { G4cout << " d-"; }
>> 864 G4cout.precision(5);
>> 865 G4cout << " new_step= " << std::setw(10)
>> 866 << fLastStepEstimate_Unconstrained
>> 867 << " new_step_constr= " << std::setw(10)
>> 868 << lastStepTrial << G4endl;
>> 869 G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
>> 870 G4cout.precision(oldprec);
655 } 871 }
656 872