Geant4-11
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4ChordFinderDelegate.icc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4ChordFinderDelegate inline methods implementation
27//
28// Created: D.Sorokin
29// --------------------------------------------------------------------
30
31#include <iomanip>
32
33template <class Driver>
34G4ChordFinderDelegate<Driver>::~G4ChordFinderDelegate()
35{
36#ifdef G4VERBOSE
37 if (GetDriver().GetVerboseLevel() > 0)
38 {
39 PrintStatistics();
40 }
41#endif
42}
43
44template <class Driver>
45void G4ChordFinderDelegate<Driver>::ResetStepEstimate()
46{
47 fLastStepEstimate_Unconstrained = DBL_MAX;
48}
49
50template <class Driver>
51Driver& G4ChordFinderDelegate<Driver>::GetDriver()
52{
53 return static_cast<Driver&>(*this);
54}
55
56template <class Driver>
57G4double G4ChordFinderDelegate<Driver>::
58AdvanceChordLimitedImpl(G4FieldTrack& yCurrent, G4double stepMax,
59 G4double epsStep, G4double chordDistance)
60{
61 G4double dyErr;
62 G4FieldTrack yEnd = yCurrent;
63 G4double nextStep;
64
65 const G4double stepPossible = FindNextChord(yCurrent, stepMax,
66 epsStep, chordDistance,
67 yEnd, dyErr, nextStep);
68 if (dyErr < epsStep * stepPossible)
69 {
70 // Accept this accuracy.
71 //
72 yCurrent = yEnd;
73 return stepPossible;
74 }
75
76 // Advance more accurately to "end of chord"
77 //
78 const G4double startCurveLen = yCurrent.GetCurveLength();
79 const G4bool goodAdvance =
80 GetDriver().AccurateAdvance(yCurrent,stepPossible,epsStep,nextStep);
81
82 return goodAdvance ? stepPossible
83 : yCurrent.GetCurveLength() - startCurveLen;
84}
85
86// Returns Length of Step taken
87//
88template <class T>
89G4double G4ChordFinderDelegate<T>::
90FindNextChord(const G4FieldTrack& yStart,
91 G4double stepMax,
92 G4double epsStep,
93 G4double chordDistance,
94 G4FieldTrack& yEnd, // Endpoint
95 G4double& dyErrPos, // Error of endpoint
96 G4double& stepForAccuracy)
97{
98 // 1.) Try to "leap" to end of interval
99 // 2.) Evaluate if resulting chord gives d_chord that is good enough.
100 // 2a.) If d_chord is not good enough, find one that is.
101
102 G4double dydx[G4FieldTrack::ncompSVEC];
103
104 G4bool validEndPoint = false;
105 G4double dChordStep, lastStepLength;
106
107 GetDriver().GetDerivatives(yStart, dydx);
108
109 const G4double safetyFactor = fFirstFraction; // 0.975 or 0.99 ? was 0.999
110
111 G4double stepTrial = std::min(stepMax,
112 safetyFactor*fLastStepEstimate_Unconstrained);
113
114 G4double newStepEst_Uncons = 0.0;
115 G4double stepForChord;
116
117 G4int noTrials = 1;
118 constexpr G4int maxTrials = 75; // Avoid endless loop for bad convergence
119 for (; noTrials < maxTrials; ++noTrials)
120 {
121 yEnd = yStart; // Always start from initial point
122 GetDriver().QuickAdvance(yEnd, dydx, stepTrial, dChordStep, dyErrPos);
123 lastStepLength = stepTrial;
124
125 validEndPoint = dChordStep < chordDistance;
126 stepForChord = NewStep(stepTrial, dChordStep,
127 chordDistance, newStepEst_Uncons);
128 if (validEndPoint)
129 {
130 break;
131 }
132
133 if (stepTrial <= 0.0)
134 {
135 stepTrial = stepForChord;
136 }
137 else if (stepForChord <= stepTrial)
138 {
139 // Reduce by a fraction, possibly up to 20%
140 stepTrial = std::min( stepForChord, fFractionLast * stepTrial);
141 }
142 else
143 {
144 stepTrial *= 0.1;
145 }
146 }
147
148 if (noTrials >= maxTrials)
149 {
150 std::ostringstream message;
151 message << "Exceeded maximum number of trials= " << maxTrials << G4endl
152 << "Current sagita dist= " << dChordStep << G4endl
153 << "Max sagita dist= " << chordDistance << G4endl
154 << "Step sizes (actual and proposed): " << G4endl
155 << "Last trial = " << lastStepLength << G4endl
156 << "Next trial = " << stepTrial << G4endl
157 << "Proposed for chord = " << stepForChord << G4endl;
158 G4Exception("G4ChordFinder::FindNextChord()", "GeomField0003",
159 JustWarning, message);
160 }
161
162 if (newStepEst_Uncons > 0.0)
163 {
164 fLastStepEstimate_Unconstrained = newStepEst_Uncons;
165 }
166
167 AccumulateStatistics(noTrials);
168
169
170 // Calculate the step size required for accuracy, if it is needed
171 G4double dyErr_relative = dyErrPos / (epsStep * lastStepLength);
172 stepForAccuracy = dyErr_relative > 1 ?
173 GetDriver().ComputeNewStepSize(dyErr_relative, lastStepLength) : 0;
174
175 return stepTrial;
176}
177
178// Is called to estimate the next step size, even for successful steps,
179// in order to predict an accurate 'chord-sensitive' first step
180// which is likely to assist in more performant 'stepping'.
181//
182template <class T>
183G4double G4ChordFinderDelegate<T>::
184NewStep(G4double stepTrialOld,
185 G4double dChordStep, // Curr. dchord achieved
186 G4double fDeltaChord,
187 G4double& stepEstimate_Unconstrained)
188{
189 G4double stepTrial;
190
191 if (dChordStep > 0.0)
192 {
193 stepEstimate_Unconstrained =
194 stepTrialOld * std::sqrt(fDeltaChord / dChordStep);
195 stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
196 }
197 else
198 {
199 // Should not update the Unconstrained Step estimate: incorrect!
200 stepTrial = stepTrialOld * 2.;
201 }
202
203 if (stepTrial <= 0.001 * stepTrialOld)
204 {
205 if (dChordStep > 1000.0 * fDeltaChord)
206 {
207 stepTrial = stepTrialOld * 0.03;
208 }
209 else
210 {
211 if (dChordStep > 100. * fDeltaChord)
212 {
213 stepTrial = stepTrialOld * 0.1;
214 }
215 else // Try halving the length until dChordStep OK
216 {
217 stepTrial = stepTrialOld * 0.5;
218 }
219 }
220 }
221 else if (stepTrial > 1000.0 * stepTrialOld)
222 {
223 stepTrial = 1000.0 * stepTrialOld;
224 }
225
226 if (stepTrial == 0.0)
227 {
228 stepTrial= 0.000001;
229 }
230
231 // A more sophisticated chord-finder could figure out a better
232 // stepTrial, from dChordStep and the required d_geometry
233 // e.g.
234 // Calculate R, r_helix (eg at orig point)
235 // if( stepTrial < 2 pi R )
236 // stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
237 // else
238 // ??
239
240 return stepTrial;
241}
242
243template <class T>
244void G4ChordFinderDelegate<T>::AccumulateStatistics(G4int noTrials)
245{
246 fTotalNoTrials += noTrials;
247 ++fNoCalls;
248
249 if (noTrials > fmaxTrials)
250 {
251 fmaxTrials = noTrials;
252 }
253}
254
255template <class T>
256void G4ChordFinderDelegate<T>::PrintStatistics()
257{
258 // Print Statistics
259 G4cout << "G4ChordFinder statistics report: \n"
260 << " No trials: " << fTotalNoTrials
261 << " No Calls: " << fNoCalls
262 << " Max-trial: " << fmaxTrials << "\n"
263 << " Parameters: "
264 << " fFirstFraction " << fFirstFraction
265 << " fFractionLast " << fFractionLast
266 << " fFractionNextEstimate " << fFractionNextEstimate
267 << G4endl;
268}
269
270template <class T>
271G4int G4ChordFinderDelegate<T>::GetNoCalls()
272{
273 return fNoCalls;
274}
275
276template <class T>
277G4int G4ChordFinderDelegate<T>::GetNoTrials()
278{
279 return fTotalNoTrials;
280}
281
282template <class T>
283G4int G4ChordFinderDelegate<T>::GetNoMaxTrials()
284{
285 return fmaxTrials;
286}
287
288template <class T>
289void G4ChordFinderDelegate<T>::SetFractions_Last_Next(G4double fractLast,
290 G4double fractNext)
291{
292 // Use -1.0 as request for Default.
293 if (fractLast == -1.0) fractLast = 1.0; // 0.9;
294 if (fractNext == -1.0) fractNext = 0.98; // 0.9;
295
296 // fFirstFraction = 0.999; // Safe value, range: ~ 0.95 - 0.999
297 if (GetDriver().GetVerboseLevel() > 0)
298 {
299 G4cout << " ChordFnd> Trying to set fractions: "
300 << " first " << fFirstFraction
301 << " last " << fractLast
302 << " next " << fractNext
303 << G4endl;
304 }
305
306 if (fractLast > 0 && fractLast <= 1)
307 {
308 fFractionLast = fractLast;
309 } else
310 {
311 std::ostringstream message;
312 message << "Invalid fraction Last = " << fractLast
313 << "; must be 0 < fractionLast <= 1 ";
314 G4Exception("G4ChordFinderDelegate::SetFractions_Last_Next()",
315 "GeomField1001", JustWarning, message);
316 }
317 if (fractNext > 0. && fractNext < 1)
318 {
319 fFractionNextEstimate = fractNext;
320 } else
321 {
322 std::ostringstream message;
323 message << "Invalid fraction Next = " << fractNext
324 << "; must be 0 < fractionNext < 1 ";
325 G4Exception("G4ChordFinderDelegate::SetFractions_Last_Next()",
326 "GeomField1001", JustWarning, message);
327 }
328}
329
330template <class T>
331void G4ChordFinderDelegate<T>::SetFirstFraction(G4double fractFirst)
332{
333 fFirstFraction = fractFirst;
334}
335
336template <class T>
337G4double G4ChordFinderDelegate<T>::GetFirstFraction()
338{
339 return fFirstFraction;
340}
341
342template <class T>
343G4double G4ChordFinderDelegate<T>::GetFractionLast()
344{
345 return fFractionLast;
346}
347
348template <class T>
349G4double G4ChordFinderDelegate<T>::GetFractionNextEstimate()
350{
351 return fFractionNextEstimate;
352}
353
354template <class T>
355G4double G4ChordFinderDelegate<T>::GetLastStepEstimateUnc()
356{
357 return fLastStepEstimate_Unconstrained;
358}
359
360template <class T>
361void G4ChordFinderDelegate<T>::SetLastStepEstimateUnc(G4double stepEst)
362{
363 fLastStepEstimate_Unconstrained = stepEst;
364}
365
366template <class T>
367void G4ChordFinderDelegate<T>::TestChordPrint(G4int noTrials,
368 G4int lastStepTrial,
369 G4double dChordStep,
370 G4double fDeltaChord,
371 G4double nextStepTrial)
372{
373 G4int oldprec = G4cout.precision(5);
374 G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
375 << " this_step= " << std::setw(10) << lastStepTrial;
376 if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
377 {
378 G4cout.precision(8);
379 }
380 else
381 {
382 G4cout.precision(6);
383 }
384 G4cout << " dChordStep= " << std::setw(12) << dChordStep;
385 if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
386 else { G4cout << " d-"; }
387 G4cout.precision(5);
388 G4cout << " new_step= " << std::setw(10)
389 << fLastStepEstimate_Unconstrained
390 << " new_step_constr= " << std::setw(10)
391 << lastStepTrial << G4endl;
392 G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
393 G4cout.precision(oldprec);
394}
395
396template <class T>
397void G4ChordFinderDelegate<T>::StreamDelegateInfo( std::ostream& os ) const
398{
399// Write out the parameters / state of the driver
400 os << "State of G4ChordFinderDelegate: " << std::endl;
401 os << "--Parameters: " << std::endl;
402 os << " First Fraction = " << fFirstFraction << std::endl;
403 os << " Last Fraction = " << fFractionLast << std::endl;
404 os << " Fract Next est = " << fFractionNextEstimate << std::endl;
405
406 os << "--State (fungible): " << std::endl;
407 os << " Maximum No Trials (seen) = " << fmaxTrials << std::endl;
408 os << " LastStepEstimate (Unconstrained) = " << fLastStepEstimate_Unconstrained
409 << std::endl;
410 // os << " Statistics NOT printed. " << std::endl;
411 os << "--Statistics: trials= " << fTotalNoTrials
412 << " calls= " << fNoCalls << std::endl;
413}
414
Generated by doxygen 1.9.3