source: trunk/source/geometry/magneticfield/src/G4ChordFinder.cc@ 834

Last change on this file since 834 was 831, checked in by garnier, 17 years ago

import all except CVS
File size: 19.1 KB
Line 
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//
27// $Id: G4ChordFinder.cc,v 1.47 2006/06/29 18:23:32 gunter Exp $
28// GEANT4 tag $Name: $
29//
30//
31// 25.02.97 John Apostolakis, design and implimentation
32// 05.03.97 V. Grichine , style modification
33// -------------------------------------------------------------------
34
35#include "G4ChordFinder.hh"
36#include "G4MagneticField.hh"
37#include "G4Mag_UsualEqRhs.hh"
38#include "G4ClassicalRK4.hh"
39
40#include <iomanip>
41
42// ..........................................................................
43
44G4ChordFinder::G4ChordFinder(G4MagInt_Driver* pIntegrationDriver)
45 : fDefaultDeltaChord( 0.25 * mm ),
46 fDeltaChord( fDefaultDeltaChord ),
47 fAllocatedStepper(false),
48 fEquation(0),
49 fDriversStepper(0),
50 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
51 fMultipleRadius(15.0),
52 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
53 fStatsVerbose(0)
54{
55 // Simple constructor which does not create equation, ..
56 // fDeltaChord= fDefaultDeltaChord;
57 fIntgrDriver= pIntegrationDriver;
58 fAllocatedStepper= false;
59 fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
60
61 SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
62 // check the values and set the other parameters
63}
64
65// ..........................................................................
66
67G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField,
68 G4double stepMinimum,
69 G4MagIntegratorStepper* pItsStepper )
70 : fDefaultDeltaChord( 0.25 * mm ),
71 fDeltaChord( fDefaultDeltaChord ),
72 fAllocatedStepper(false),
73 fEquation(0),
74 fDriversStepper(0),
75 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
76 fMultipleRadius(15.0),
77 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
78 fStatsVerbose(0)
79{
80 // Construct the Chord Finder
81 // by creating in inverse order the Driver, the Stepper and EqRhs ...
82 G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
83 fEquation = pEquation;
84 fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
85 // G4FieldTrack ??
86
87 SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
88 // check the values and set the other parameters
89
90 // --->> Charge Q = 0
91 // --->> Momentum P = 1 NOMINAL VALUES !!!!!!!!!!!!!!!!!!
92
93 if( pItsStepper == 0 )
94 {
95 pItsStepper = fDriversStepper = new G4ClassicalRK4(pEquation);
96 fAllocatedStepper= true;
97 }
98 else
99 {
100 fAllocatedStepper= false;
101 }
102 fIntgrDriver = new G4MagInt_Driver(stepMinimum, pItsStepper,
103 pItsStepper->GetNumberOfVariables() );
104}
105
106// ......................................................................
107
108void
109G4ChordFinder::SetFractions_Last_Next( G4double fractLast, G4double fractNext )
110{
111 // Use -1.0 as request for Default.
112 if( fractLast == -1.0 ) fractLast = 1.0; // 0.9;
113 if( fractNext == -1.0 ) fractNext = 0.98; // 0.9;
114
115 // fFirstFraction = 0.999; // Orig 0.999 A safe value, range: ~ 0.95 - 0.999
116 // fMultipleRadius = 15.0; // For later use, range: ~ 2 - 20
117
118 if( fStatsVerbose ) {
119 G4cout << " ChordFnd> Trying to set fractions: "
120 << " first " << fFirstFraction
121 << " last " << fractLast
122 << " next " << fractNext
123 << " and multiple " << fMultipleRadius
124 << G4endl;
125 }
126
127 if( (fractLast > 0.0) && (fractLast <=1.0) )
128 { fFractionLast= fractLast; }
129 else
130 G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
131 << " fraction Last = " << fractLast
132 << " must be 0 < fractionLast <= 1 " << G4endl;
133 if( (fractNext > 0.0) && (fractNext <1.0) )
134 { fFractionNextEstimate = fractNext; }
135 else
136 G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
137 << " fraction Next = " << fractNext
138 << " must be 0 < fractionNext < 1 " << G4endl;
139}
140
141// ......................................................................
142
143G4ChordFinder::~G4ChordFinder()
144{
145 delete fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
146 if( fAllocatedStepper)
147 {
148 delete fDriversStepper;
149 } // fIntgrDriver->pIntStepper;}
150 delete fIntgrDriver;
151
152 if( fStatsVerbose ) { PrintStatistics(); }
153}
154
155void
156G4ChordFinder::PrintStatistics()
157{
158 // Print Statistics
159 G4cout << "G4ChordFinder statistics report: " << G4endl;
160 G4cout
161 << " No trials: " << fTotalNoTrials_FNC
162 << " No Calls: " << fNoCalls_FNC
163 << " Max-trial: " << fmaxTrials_FNC
164 << G4endl;
165 G4cout
166 << " Parameters: "
167 << " fFirstFraction " << fFirstFraction
168 << " fFractionLast " << fFractionLast
169 << " fFractionNextEstimate " << fFractionNextEstimate
170 << G4endl;
171}
172
173// ......................................................................
174
175G4double
176G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
177 G4double stepMax,
178 G4double epsStep,
179 const G4ThreeVector latestSafetyOrigin,
180 G4double latestSafetyRadius
181 )
182{
183 G4double stepPossible;
184 G4double dyErr;
185 G4FieldTrack yEnd( yCurrent);
186 G4double startCurveLen= yCurrent.GetCurveLength();
187
188 G4double nextStep;
189 // *************
190 stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep, &nextStep
191 , latestSafetyOrigin, latestSafetyRadius
192 );
193 // *************
194 G4bool good_advance;
195 if ( dyErr < epsStep * stepPossible )
196 {
197 // Accept this accuracy.
198 yCurrent = yEnd;
199 good_advance = true;
200 }
201 else
202 {
203 // Advance more accurately to "end of chord"
204 // ***************
205 good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible, epsStep, nextStep);
206 // ***************
207 if ( ! good_advance ){
208 // In this case the driver could not do the full distance
209 stepPossible= yCurrent.GetCurveLength()-startCurveLen;
210 }
211 }
212
213#ifdef G4DEBUG_FIELD
214 G4cout << "Exiting FindNextChord Limited with:" << G4endl
215 << " yCurrent: " << yCurrent<< G4endl;
216#endif
217
218 return stepPossible;
219}
220
221// #define TEST_CHORD_PRINT 1
222
223// ............................................................................
224
225G4double
226G4ChordFinder::FindNextChord( const G4FieldTrack yStart,
227 G4double stepMax,
228 G4FieldTrack& yEnd, // Endpoint
229 G4double& dyErrPos, // Error of endpoint
230 G4double epsStep,
231 G4double* pStepForAccuracy,
232 const G4ThreeVector, // latestSafetyOrigin,
233 G4double // latestSafetyRadius
234 )
235 // Returns Length of Step taken
236{
237 // G4int stepRKnumber=0;
238 G4FieldTrack yCurrent= yStart;
239 G4double stepTrial, stepForAccuracy;
240 G4double dydx[G4FieldTrack::ncompSVEC];
241
242 // 1.) Try to "leap" to end of interval
243 // 2.) Evaluate if resulting chord gives d_chord that is good enough.
244 // 2a.) If d_chord is not good enough, find one that is.
245
246 G4bool validEndPoint= false;
247 G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
248
249 fIntgrDriver-> GetDerivatives( yCurrent, dydx ) ;
250
251 G4int noTrials=0;
252 const G4double safetyFactor= fFirstFraction; // 0.975 or 0.99 ? was 0.999
253
254 stepTrial = std::min( stepMax,
255 safetyFactor * fLastStepEstimate_Unconstrained );
256
257 G4double newStepEst_Uncons= 0.0;
258 do
259 {
260 G4double stepForChord;
261 yCurrent = yStart; // Always start from initial point
262
263 // ************
264 fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
265 dChordStep, dyErrPos);
266 // ************
267
268 // We check whether the criterion is met here.
269 validEndPoint = AcceptableMissDist(dChordStep);
270 // && (dyErrPos < eps) ;
271
272 lastStepLength = stepTrial;
273
274 // This method estimates to step size for a good chord.
275 stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
276
277 if( ! validEndPoint ) {
278 if( stepTrial<=0.0 )
279 stepTrial = stepForChord;
280 else if (stepForChord <= stepTrial)
281 // Reduce by a fraction, possibly up to 20%
282 stepTrial = std::min( stepForChord,
283 fFractionLast * stepTrial);
284 else
285 stepTrial *= 0.1;
286
287 // if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
288 }
289 // #ifdef TEST_CHORD_PRINT
290 // TestChordPrint( noTrials, lastStepLength, dChordStep, stepTrial );
291 // #endif
292
293 noTrials++;
294 }
295 while( ! validEndPoint ); // End of do-while RKD
296
297 if( newStepEst_Uncons > 0.0 ){
298 fLastStepEstimate_Unconstrained= newStepEst_Uncons;
299 }
300
301 AccumulateStatistics( noTrials );
302
303 // stepOfLastGoodChord = stepTrial;
304
305 if( pStepForAccuracy ){
306 // Calculate the step size required for accuracy, if it is needed
307 G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
308 if( dyErr_relative > 1.0 ) {
309 stepForAccuracy =
310 fIntgrDriver->ComputeNewStepSize( dyErr_relative,
311 lastStepLength );
312 }else{
313 stepForAccuracy = 0.0; // Convention to show step was ok
314 }
315 *pStepForAccuracy = stepForAccuracy;
316 }
317
318#ifdef TEST_CHORD_PRINT
319 static int dbg=0;
320 if( dbg )
321 G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
322 << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
323#endif
324
325 yEnd= yCurrent;
326 return stepTrial;
327}
328
329// ----------------------------------------------------------------------------
330#if 0
331// #ifdef G4VERBOSE
332if( dbg ) {
333 G4cerr << "Returned from QuickAdvance with: yCur=" << yCurrent <<G4endl;
334 G4cerr << " dChordStep= "<< dChordStep <<" dyErr=" << dyErr << G4endl;
335}
336#endif
337// ----------------------------------------------------------------------------
338
339// ...........................................................................
340
341G4double G4ChordFinder::NewStep(G4double stepTrialOld,
342 G4double dChordStep, // Curr. dchord achieved
343 G4double& stepEstimate_Unconstrained )
344//
345// Is called to estimate the next step size, even for successful steps,
346// in order to predict an accurate 'chord-sensitive' first step
347// which is likely to assist in more performant 'stepping'.
348//
349
350{
351 G4double stepTrial;
352 static G4double lastStepTrial = 1., lastDchordStep= 1.;
353
354#if 1
355 // const G4double threshold = 1.21, multiplier = 0.9;
356 // 0.9 < 1 / std::sqrt(1.21)
357
358 if (dChordStep > 0.0)
359 {
360 stepEstimate_Unconstrained = stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
361 // stepTrial = 0.98 * stepEstimate_Unconstrained;
362 stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
363 }
364 else
365 {
366 // Should not update the Unconstrained Step estimate: incorrect!
367 stepTrial = stepTrialOld * 2.;
368 }
369
370 // if ( dChordStep < threshold * fDeltaChord ){
371 // stepTrial= stepTrialOld * multiplier;
372 // }
373 if( stepTrial <= 0.001 * stepTrialOld)
374 {
375 if ( dChordStep > 1000.0 * fDeltaChord ){
376 stepTrial= stepTrialOld * 0.03;
377 }else{
378 if ( dChordStep > 100. * fDeltaChord ){
379 stepTrial= stepTrialOld * 0.1;
380 }else{
381 // Try halving the length until dChordStep OK
382 stepTrial= stepTrialOld * 0.5;
383 }
384 }
385 }else if (stepTrial > 1000.0 * stepTrialOld)
386 {
387 stepTrial= 1000.0 * stepTrialOld;
388 }
389
390 if( stepTrial == 0.0 ){
391 stepTrial= 0.000001;
392 }
393
394 lastStepTrial = stepTrialOld;
395 lastDchordStep= dChordStep;
396#else
397 if ( dChordStep > 1000. * fDeltaChord ){
398 stepTrial= stepTrialOld * 0.03;
399 }else{
400 if ( dChordStep > 100. * fDeltaChord ){
401 stepTrial= stepTrialOld * 0.1;
402 }else{
403 // Keep halving the length until dChordStep OK
404 stepTrial= stepTrialOld * 0.5;
405 }
406 }
407#endif
408
409 // A more sophisticated chord-finder could figure out a better
410 // stepTrial, from dChordStep and the required d_geometry
411 // eg
412 // Calculate R, r_helix (eg at orig point)
413 // if( stepTrial < 2 pi R )
414 // stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
415 // else
416 // ??
417
418 return stepTrial;
419}
420
421//
422// Given a starting curve point A (CurveA_PointVelocity), a later
423// curve point B (CurveB_PointVelocity) and a point E which is (generally)
424// not on the curve, find and return a point F which is on the curve and
425// which is close to E. While advancing towards F utilise eps_step
426// as a measure of the relative accuracy of each Step.
427
428G4FieldTrack
429G4ChordFinder::ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
430 const G4FieldTrack& CurveB_PointVelocity,
431 const G4ThreeVector& CurrentE_Point,
432 G4double eps_step)
433{
434 // 1st implementation:
435 // if r=|AE|/|AB|, and s=true path lenght (AB)
436 // return the point that is r*s along the curve!
437
438 G4FieldTrack Current_PointVelocity= CurveA_PointVelocity;
439
440 G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
441 G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
442
443 G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
444 G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
445
446 G4double ABdist= ChordAB_Vector.mag();
447 G4double curve_length; // A curve length of AB
448 G4double AE_fraction;
449
450 curve_length= CurveB_PointVelocity.GetCurveLength()
451 - CurveA_PointVelocity.GetCurveLength();
452
453 // const
454 G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
455 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) ){
456#ifdef G4DEBUG_FIELD
457 G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
458 << G4endl
459 << " The two points are further apart than the curve length "
460 << G4endl
461 << " Dist = " << ABdist
462 << " curve length = " << curve_length
463 << " relativeDiff = " << (curve_length-ABdist)/ABdist
464 << G4endl;
465 if( curve_length < ABdist * (1. - 10*eps_step) ) {
466 G4cerr << " ERROR: the size of the above difference"
467 << " exceeds allowed limits. Aborting." << G4endl;
468 G4Exception("G4ChordFinder::ApproxCurvePointV()", "PrecisionError",
469 FatalException, "Unphysical curve length.");
470 }
471#endif
472 // Take default corrective action:
473 // --> adjust the maximum curve length.
474 // NOTE: this case only happens for relatively straight paths.
475 curve_length = ABdist;
476 }
477
478 G4double new_st_length;
479
480 if ( ABdist > 0.0 ){
481 AE_fraction = ChordAE_Vector.mag() / ABdist;
482 }else{
483 AE_fraction = 0.5; // Guess .. ?;
484#ifdef G4DEBUG_FIELD
485 G4cout << "Warning in G4ChordFinder::ApproxCurvePoint:"
486 << " A and B are the same point!" << G4endl
487 << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
488 << G4endl;
489#endif
490 }
491
492 if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) ){
493#ifdef G4DEBUG_FIELD
494 G4cerr << " G4ChordFinder::ApproxCurvePointV - Warning:"
495 << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
496 << " AE_fraction = " << AE_fraction << G4endl
497 << " Chord AE length = " << ChordAE_Vector.mag() << G4endl
498 << " Chord AB length = " << ABdist << G4endl << G4endl;
499 G4cerr << " OK if this condition occurs after a recalculation of 'B'"
500 << G4endl << " Otherwise it is an error. " << G4endl ;
501#endif
502 // This course can now result if B has been re-evaluated,
503 // without E being recomputed (1 July 99)
504 // In this case this is not a "real error" - but it undesired
505 // and we cope with it by a default corrective action ...
506 AE_fraction = 0.5; // Default value
507 }
508
509 new_st_length= AE_fraction * curve_length;
510
511 G4bool good_advance;
512 if ( AE_fraction > 0.0 ) {
513 good_advance =
514 fIntgrDriver->AccurateAdvance(Current_PointVelocity,
515 new_st_length,
516 eps_step ); // Relative accuracy
517 // In this case it does not matter if it cannot advance the full distance
518 }
519
520 // If there was a memory of the step_length actually require at the start
521 // of the integration Step, this could be re-used ...
522
523 return Current_PointVelocity;
524}
525
526void
527G4ChordFinder::TestChordPrint( G4int noTrials,
528 G4int lastStepTrial,
529 G4double dChordStep,
530 G4double nextStepTrial )
531{
532 G4int oldprec= G4cout.precision(5);
533 G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
534 << " this_step= " << std::setw(10) << lastStepTrial;
535 if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 ){
536 G4cout.precision(8);
537 }else{ G4cout.precision(6); }
538 G4cout << " dChordStep= " << std::setw(12) << dChordStep;
539 if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
540 else { G4cout << " d-"; }
541 G4cout.precision(5);
542 G4cout << " new_step= " << std::setw(10)
543 << fLastStepEstimate_Unconstrained
544 << " new_step_constr= " << std::setw(10)
545 << lastStepTrial << G4endl;
546 G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
547 G4cout.precision(oldprec);
548}
Note: See TracBrowser for help on using the repository browser.