Geant4  10.01
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 
28 #include "G4ChordFinderSaf.hh"
29 #include <iomanip>
30 
31 // ..........................................................................
32 
34  : G4ChordFinder(pIntegrationDriver)
35 {
36  // check the values and set the other parameters
37  // fNoInitialRadBig=0; fNoInitialRadSmall=0;
38  // fNoTrialsRadBig=0; fNoTrialsRadSmall=0;
39 
40 }
41 
42 // ..........................................................................
43 
45  G4double stepMinimum,
46  G4MagIntegratorStepper* pItsStepper )
47  : G4ChordFinder( theMagField, stepMinimum, pItsStepper )
48 {
49  // Let G4ChordFinder create the Driver, the Stepper and EqRhs ...
50  // ...
51 }
52 
53 // ......................................................................
54 
55 // ......................................................................
56 
58 {
59  if( SetVerbose(0) ) { PrintStatistics(); }
60  // Set verbosity 0, so that will be called in base class again
61 }
62 
63 void
65 {
66  // Print Statistics
67  G4cout << "G4ChordFinderSaf statistics report: " << G4endl;
69 
70 /*******************
71  G4cout
72  << " No radbig calls " << std::setw(10) << fNoInitialRadBig
73  << " trials " << std::setw(10) << fNoTrialsRadBig
74  << " - over " << std::setw(10) << fNoTrialsRadBig - fNoInitialRadBig
75  << G4endl
76  << " No radsm calls " << std::setw(10) << fNoInitialRadSmall
77  << " trials " << std::setw(10) << fNoTrialsRadSmall
78  << " - over " << std::setw(10) << fNoTrialsRadSmall - fNoInitialRadSmall
79  << G4endl;
80  G4cout
81  << " *** Limiting stepTrial via if Delta_chord < R_curvature "
82  << " for large to angle from Delta_chord / R_curv "
83  << " and for small with multiple " << GetMultipleRadius()
84  << G4endl;
85 ********************/
86 }
87 
88 
89 // G4SafetyDist::
90 // inline
93  G4double safetyRadius,
94  G4ThreeVector point)
95 {
96  G4double pointSafety= 0.0;
97 
98  G4ThreeVector OriginShift = point - safetyOrigin ;
99  G4double MagSqShift = OriginShift.mag2() ;
100  if( MagSqShift < sqr(safetyRadius) ){
101  pointSafety = safetyRadius - std::sqrt(MagSqShift) ;
102  }
103 
104  return pointSafety;
105 }
106 
107 // inline
108 G4bool
110  G4double safetyRadius,
111  G4ThreeVector point)
112 {
113  G4ThreeVector OriginShift = point - safetyOrigin ;
114  return ( OriginShift.mag2() < safetyRadius*safetyRadius );
115 }
116 
117 G4double
119  G4double stepMax,
120  G4FieldTrack& yEnd, // Endpoint
121  G4double& dyErrPos, // Error of endpoint
122  G4double epsStep,
123  G4double* pStepForAccuracy,
124  const G4ThreeVector latestSafetyOrigin,
125  G4double latestSafetyRadius
126  )
127  // Returns Length of Step taken
128 {
129  // G4int stepRKnumber=0;
130  G4FieldTrack yCurrent= yStart;
131  G4double stepTrial, stepForAccuracy;
133 
134  // 1.) Try to "leap" to end of interval
135  // 2.) Evaluate if resulting chord gives d_chord that is good enough.
136  // 2a.) If d_chord is not good enough, find one that is.
137 
138  G4bool validEndPoint= false;
139  G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
140 
141  GetIntegrationDriver()-> GetDerivatives( yCurrent, dydx ) ;
142 
143  G4int noTrials=0;
144  const G4double safetyFactor= GetFirstFraction(); // was 0.999
145 
146  // Figure out the starting safety
147  G4double startSafety=
148  CalculatePointSafety( latestSafetyOrigin,
149  latestSafetyRadius,
150  yCurrent.GetPosition() );
151 
152  G4double
153  likelyGood = std::max( startSafety ,
154  safetyFactor * GetLastStepEstimateUnc() );
155 
156  stepTrial = std::min( stepMax, likelyGood );
157 
159  G4double newStepEst_Uncons= 0.0;
160  do
161  {
162  G4double stepForChord;
163  yCurrent = yStart; // Always start from initial point
164 
165  // ************
166  pIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
167  dChordStep, dyErrPos);
168  // ************
169 
170  G4ThreeVector EndPointCand= yCurrent.GetPosition();
171  G4bool endPointInSafetySphere=
172  CalculatePointInside(latestSafetyOrigin, latestSafetyRadius, EndPointCand);
173 
174  // We check whether the criterion is met here.
175  validEndPoint = AcceptableMissDist(dChordStep)
176  || endPointInSafetySphere;
177  // && (dyErrPos < eps) ;
178 
179  lastStepLength = stepTrial;
180 
181  // This method estimates to step size for a good chord.
182  stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
183 
184  if( ! validEndPoint ) {
185  if( stepTrial<=0.0 )
186  stepTrial = stepForChord;
187  else if (stepForChord <= stepTrial)
188  // Reduce by a fraction, possibly up to 20%
189  stepTrial = std::min( stepForChord,
190  GetFractionLast() * stepTrial);
191  else
192  stepTrial *= 0.1;
193 
194  // if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
195  }
196  // #ifdef TEST_CHORD_PRINT
197  // TestChordPrint( noTrials, lastStepLength, dChordStep, stepTrial );
198  // #endif
199 
200  noTrials++;
201  }
202  while( ! validEndPoint ); // End of do-while RKD
203 
204  AccumulateStatistics( noTrials );
205 
206  // Should we update newStepEst_Uncons for a 'long step' via safety ??
207  if( newStepEst_Uncons > 0.0 ){
208  SetLastStepEstimateUnc( newStepEst_Uncons );
209  }
210 
211  // stepOfLastGoodChord = stepTrial;
212 
213  if( pStepForAccuracy ){
214  // Calculate the step size required for accuracy, if it is needed
215  G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
216  if( dyErr_relative > 1.0 ) {
217  stepForAccuracy =
218  pIntgrDriver->ComputeNewStepSize( dyErr_relative,
219  lastStepLength );
220  }else{
221  stepForAccuracy = 0.0; // Convention to show step was ok
222  }
223  *pStepForAccuracy = stepForAccuracy;
224  }
225 
226 #ifdef TEST_CHORD_PRINT
227  static int dbg=0;
228  if( dbg )
229  G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
230  << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
231 #endif
232 
233  yEnd= yCurrent;
234  return stepTrial;
235 }
G4double FindNextChord(const G4FieldTrack &yStart, G4double stepMax, G4FieldTrack &yEnd, G4double &dyErrPos, G4double epsStep, G4double *pStepForAccuracy, const G4ThreeVector latestSafetyOrigin, G4double latestSafetyRadius)
G4double GetFirstFraction()
CLHEP::Hep3Vector G4ThreeVector
void SetLastStepEstimateUnc(G4double stepEst)
int G4int
Definition: G4Types.hh:78
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)
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
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)