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