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