Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ChordFinderSaf.cc
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 
27 #include "G4ChordFinderSaf.hh"
28 #include <iomanip>
29 
30 // ..........................................................................
31 
33  : G4ChordFinder(pIntegrationDriver)
34 {
35  // check the values and set the other parameters
36  // fNoInitialRadBig=0; fNoInitialRadSmall=0;
37  // fNoTrialsRadBig=0; fNoTrialsRadSmall=0;
38 
39 }
40 
41 // ..........................................................................
42 
44  G4double stepMinimum,
45  G4MagIntegratorStepper* pItsStepper )
46  : G4ChordFinder( theMagField, stepMinimum, pItsStepper )
47 {
48  // Let G4ChordFinder create the Driver, the Stepper and EqRhs ...
49  // ...
50 }
51 
52 // ......................................................................
53 
54 // ......................................................................
55 
57 {
58  if( SetVerbose(0) ) { PrintStatistics(); }
59  // Set verbosity 0, so that will be called in base class again
60 }
61 
62 void
64 {
65  // Print Statistics
66  G4cout << "G4ChordFinderSaf statistics report: " << G4endl;
68 
69 /*******************
70  G4cout
71  << " No radbig calls " << std::setw(10) << fNoInitialRadBig
72  << " trials " << std::setw(10) << fNoTrialsRadBig
73  << " - over " << std::setw(10) << fNoTrialsRadBig - fNoInitialRadBig
74  << G4endl
75  << " No radsm calls " << std::setw(10) << fNoInitialRadSmall
76  << " trials " << std::setw(10) << fNoTrialsRadSmall
77  << " - over " << std::setw(10) << fNoTrialsRadSmall - fNoInitialRadSmall
78  << G4endl;
79  G4cout
80  << " *** Limiting stepTrial via if Delta_chord < R_curvature "
81  << " for large to angle from Delta_chord / R_curv "
82  << " and for small with multiple " << GetMultipleRadius()
83  << G4endl;
84 ********************/
85 }
86 
87 
88 // G4SafetyDist::
89 // inline
92  G4double safetyRadius,
93  G4ThreeVector point)
94 {
95  G4double pointSafety= 0.0;
96 
97  G4ThreeVector OriginShift = point - safetyOrigin ;
98  G4double MagSqShift = OriginShift.mag2() ;
99  if( MagSqShift < sqr(safetyRadius) ){
100  pointSafety = safetyRadius - std::sqrt(MagSqShift) ;
101  }
102 
103  return pointSafety;
104 }
105 
106 // inline
107 G4bool
109  G4double safetyRadius,
110  G4ThreeVector point)
111 {
112  G4ThreeVector OriginShift = point - safetyOrigin ;
113  return ( OriginShift.mag2() < safetyRadius*safetyRadius );
114 }
115 
116 G4double
118  G4double stepMax,
119  G4FieldTrack& yEnd, // Endpoint
120  G4double& dyErrPos, // Error of endpoint
121  G4double epsStep,
122  G4double* pStepForAccuracy,
123  const G4ThreeVector latestSafetyOrigin,
124  G4double latestSafetyRadius
125  )
126  // Returns Length of Step taken
127 {
128  // G4int stepRKnumber=0;
129  G4FieldTrack yCurrent= yStart;
130  G4double stepTrial, stepForAccuracy;
132 
133  // 1.) Try to "leap" to end of interval
134  // 2.) Evaluate if resulting chord gives d_chord that is good enough.
135  // 2a.) If d_chord is not good enough, find one that is.
136 
137  G4bool validEndPoint= false;
138  G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
139 
140  GetIntegrationDriver()-> GetDerivatives( yCurrent, dydx ) ;
141 
142  unsigned int noTrials=0;
143  const unsigned int maxTrials= 75; // Avoid endless loop for bad convergence
144 
145  const G4double safetyFactor= GetFirstFraction(); // was 0.999
146 
147  // Figure out the starting safety
148  G4double startSafety=
149  CalculatePointSafety( latestSafetyOrigin,
150  latestSafetyRadius,
151  yCurrent.GetPosition() );
152 
153  G4double
154  likelyGood = std::max( startSafety ,
155  safetyFactor * GetLastStepEstimateUnc() );
156 
157  stepTrial = std::min( stepMax, likelyGood );
158 
160  G4double newStepEst_Uncons= 0.0;
161  G4double stepForChord= -1.0;
162  do
163  {
164  yCurrent = yStart; // Always start from initial point
165 
166  // ************
167  pIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
168  dChordStep, dyErrPos);
169  // ************
170 
171  G4ThreeVector EndPointCand= yCurrent.GetPosition();
172  G4bool endPointInSafetySphere=
173  CalculatePointInside(latestSafetyOrigin, latestSafetyRadius, EndPointCand);
174 
175  // We check whether the criterion is met here.
176  validEndPoint = AcceptableMissDist(dChordStep)
177  || endPointInSafetySphere;
178  // && (dyErrPos < eps) ;
179 
180  lastStepLength = stepTrial;
181 
182  // This method estimates to step size for a good chord.
183  stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
184 
185  if( ! validEndPoint ) {
186  if( stepTrial<=0.0 )
187  stepTrial = stepForChord;
188  else if (stepForChord <= stepTrial)
189  // Reduce by a fraction, possibly up to 20%
190  stepTrial = std::min( stepForChord,
191  GetFractionLast() * stepTrial);
192  else
193  stepTrial *= 0.1;
194 
195  // if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
196  }
197  // #ifdef TEST_CHORD_PRINT
198  // TestChordPrint( noTrials, lastStepLength, dChordStep, stepTrial );
199  // #endif
200 
201  noTrials++;
202  }
203  while( (! validEndPoint) && (noTrials < maxTrials) );
204  // Loop checking, 07.10.2016, J. Apostolakis
205 
206  if( noTrials >= maxTrials )
207  {
208  std::ostringstream message;
209  message << "Exceeded maximum number of trials= " << maxTrials << G4endl
210  << "Current sagita dist= " << dChordStep << G4endl
211  << "Step sizes (actual and proposed): " << G4endl
212  << "Last trial = " << lastStepLength << G4endl
213  << "Next trial = " << stepTrial << G4endl
214  << "Proposed for chord = " << stepForChord << G4endl
215  ;
216  G4Exception("G4ChordFinderSaf::FindNextChord()", "GeomField0003",
217  JustWarning, message);
218  }
219  AccumulateStatistics( noTrials );
220 
221  // Should we update newStepEst_Uncons for a 'long step' via safety ??
222  if( newStepEst_Uncons > 0.0 ){
223  SetLastStepEstimateUnc( newStepEst_Uncons );
224  }
225 
226  // stepOfLastGoodChord = stepTrial;
227 
228  if( pStepForAccuracy ){
229  // Calculate the step size required for accuracy, if it is needed
230  G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
231  if( dyErr_relative > 1.0 ) {
232  stepForAccuracy =
233  pIntgrDriver->ComputeNewStepSize( dyErr_relative,
234  lastStepLength );
235  }else{
236  stepForAccuracy = 0.0; // Convention to show step was ok
237  }
238  *pStepForAccuracy = stepForAccuracy;
239  }
240 
241 #ifdef TEST_CHORD_PRINT
242  static int dbg=0;
243  if( dbg )
244  G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
245  << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
246 #endif
247 
248  yEnd= yCurrent;
249  return stepTrial;
250 }
G4double FindNextChord(const G4FieldTrack &yStart, G4double stepMax, G4FieldTrack &yEnd, G4double &dyErrPos, G4double epsStep, G4double *pStepForAccuracy, const G4ThreeVector latestSafetyOrigin, G4double latestSafetyRadius)
G4double GetFirstFraction()
void SetLastStepEstimateUnc(G4double stepEst)
virtual void PrintStatistics()
G4int SetVerbose(G4int newvalue=1)
G4bool CalculatePointInside(G4ThreeVector safetyOrigin, G4double safetyRadius, G4ThreeVector point)
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4bool AcceptableMissDist(G4double dChordStep) const
void AccumulateStatistics(G4int noTrials)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double NewStep(G4double stepTrialOld, G4double dChordStep, G4double &stepEstimate_Unconstrained)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double mag2() const
G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)
G4double GetLastStepEstimateUnc()
#define G4endl
Definition: G4ios.hh:61
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4MagInt_Driver * GetIntegrationDriver()
G4double GetFractionLast()
G4double CalculatePointSafety(G4ThreeVector safetyOrigin, G4double safetyRadius, G4ThreeVector point)
G4ChordFinderSaf(G4MagInt_Driver *pIntegrationDriver)