| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4MagInt_Driver implementation <<
27 // 26 //
28 // V.Grichine, 07.10.1996 - Created << 27 // $Id: G4MagIntegratorDriver.cc,v 1.57 2010-07-14 10:00:36 gcosmo Exp $
29 // W.Wander, 28.01.1998 - Added ability for lo << 28 // GEANT4 tag $Name: not supported by cvs2svn $
30 // J.Apostolakis, 08.11.2001 - Respect minimum << 29 //
>> 30 //
>> 31 //
>> 32 // Implementation for class G4MagInt_Driver
>> 33 // Tracking in space dependent magnetic field
>> 34 //
>> 35 // History of major changes:
>> 36 // 8 Nov 01 J. Apostolakis: Respect minimum step in AccurateAdvance
>> 37 // 27 Jul 99 J. Apostolakis: Ensured that AccurateAdvance does not loop
>> 38 // due to very small eps & step size (precision)
>> 39 // 28 Jan 98 W. Wander: Added ability for low order integrators
>> 40 // 7 Oct 96 V. Grichine First version
31 // ------------------------------------------- 41 // --------------------------------------------------------------------
32 42
33 #include <iomanip> <<
34 <<
35 #include "globals.hh" 43 #include "globals.hh"
36 #include "G4SystemOfUnits.hh" <<
37 #include "G4GeometryTolerance.hh" 44 #include "G4GeometryTolerance.hh"
>> 45 #include <iomanip>
38 #include "G4MagIntegratorDriver.hh" 46 #include "G4MagIntegratorDriver.hh"
39 #include "G4FieldTrack.hh" 47 #include "G4FieldTrack.hh"
40 48
41 #ifdef G4DEBUG_FIELD << 49 // Stepsize can increase by no more than 5.0
42 #include "G4DriverReporter.hh" << 50 // and decrease by no more than 1/10. = 0.1
>> 51 //
>> 52 const G4double G4MagInt_Driver::max_stepping_increase = 5.0;
>> 53 const G4double G4MagInt_Driver::max_stepping_decrease = 0.1;
>> 54
>> 55 // The (default) maximum number of steps is Base
>> 56 // divided by the order of Stepper
>> 57 //
>> 58 const G4int G4MagInt_Driver::fMaxStepBase = 250; // Was 5000
>> 59
>> 60 #ifndef G4NO_FIELD_STATISTICS
>> 61 #define G4FLD_STATS 1
43 #endif 62 #endif
44 63
45 // ------------------------------------------- 64 // ---------------------------------------------------------
46 65
47 // Constructor 66 // Constructor
48 // 67 //
49 G4MagInt_Driver::G4MagInt_Driver( G4double 68 G4MagInt_Driver::G4MagInt_Driver( G4double hminimum,
50 G4MagIntegra << 69 G4MagIntegratorStepper *pStepper,
51 G4int 70 G4int numComponents,
52 G4int 71 G4int statisticsVerbose)
53 : fNoIntegrationVariables(numComponents), << 72 : fSmallestFraction( 1.0e-12 ),
>> 73 fNoIntegrationVariables(numComponents),
>> 74 fMinNoVars(12),
54 fNoVars( std::max( fNoIntegrationVariables 75 fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
55 fStatisticsVerboseLevel(statisticsVerbose) << 76 fStatisticsVerboseLevel(statisticsVerbose),
>> 77 fNoTotalSteps(0), fNoBadSteps(0), fNoSmallSteps(0),
>> 78 fNoInitialSmallSteps(0),
>> 79 fDyerr_max(0.0), fDyerr_mx2(0.0),
>> 80 fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
>> 81 fSumH_sm(0.0), fSumH_lg(0.0),
>> 82 fVerboseLevel(0)
56 { 83 {
57 // In order to accomodate "Laboratory Time", 84 // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
58 // is required. For proper time of flight an 85 // is required. For proper time of flight and spin, fMinNoVars must be 12
59 86
60 RenewStepperAndAdjust( pStepper ); 87 RenewStepperAndAdjust( pStepper );
61 fMinimumStep = hminimum; << 88 fMinimumStep= hminimum;
62 <<
63 fMaxNoSteps = fMaxStepBase / pIntStepper->In 89 fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
64 #ifdef G4DEBUG_FIELD 90 #ifdef G4DEBUG_FIELD
65 fVerboseLevel=2; 91 fVerboseLevel=2;
66 #endif 92 #endif
67 93
68 if( (fVerboseLevel > 0) || (fStatisticsVerbo 94 if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
69 { 95 {
70 G4cout << "MagIntDriver version: Accur-Adv 96 G4cout << "MagIntDriver version: Accur-Adv: "
71 << "invE_nS, QuickAdv-2sqrt with St 97 << "invE_nS, QuickAdv-2sqrt with Statistics "
72 #ifdef G4FLD_STATS 98 #ifdef G4FLD_STATS
73 << " enabled " 99 << " enabled "
74 #else 100 #else
75 << " disabled " 101 << " disabled "
76 #endif 102 #endif
77 << G4endl; 103 << G4endl;
78 } 104 }
79 } 105 }
80 106
81 // ------------------------------------------- 107 // ---------------------------------------------------------
82 108
83 // Destructor 109 // Destructor
84 // 110 //
85 G4MagInt_Driver::~G4MagInt_Driver() 111 G4MagInt_Driver::~G4MagInt_Driver()
86 { 112 {
87 if( fStatisticsVerboseLevel > 1 ) 113 if( fStatisticsVerboseLevel > 1 )
88 { 114 {
89 PrintStatisticsReport(); 115 PrintStatisticsReport();
90 } 116 }
91 } 117 }
92 118
>> 119 // To add much printing for debugging purposes, uncomment the following
>> 120 // and set verbose level to 1 or higher value !
>> 121 // #define G4DEBUG_FIELD 1
>> 122
93 // ------------------------------------------- 123 // ---------------------------------------------------------
94 124
95 G4bool 125 G4bool
96 G4MagInt_Driver::AccurateAdvance(G4FieldTrack& 126 G4MagInt_Driver::AccurateAdvance(G4FieldTrack& y_current,
97 G4double << 127 G4double hstep,
98 G4double << 128 G4double eps,
99 G4double << 129 G4double hinitial )
100 { 130 {
101 // Runge-Kutta driver with adaptive stepsize 131 // Runge-Kutta driver with adaptive stepsize control. Integrate starting
102 // values at y_current over hstep x2 with ac 132 // values at y_current over hstep x2 with accuracy eps.
103 // On output ystart is replaced by values at 133 // On output ystart is replaced by values at the end of the integration
104 // interval. RightHandSide is the right-hand 134 // interval. RightHandSide is the right-hand side of ODE system.
105 // The source is similar to odeint routine f 135 // The source is similar to odeint routine from NRC p.721-722 .
106 136
107 G4int nstp, i; << 137 G4int nstp, i, no_warnings=0;
108 G4double x, hnext, hdid, h; 138 G4double x, hnext, hdid, h;
109 139
110 #ifdef G4DEBUG_FIELD 140 #ifdef G4DEBUG_FIELD
111 G4int no_warnings = 0; << 141 static G4int dbg=1;
112 static G4int dbg = 1; << 142 static G4int nStpPr=50; // For debug printing of long integrations
113 static G4int nStpPr = 50; // For debug pri <<
114 G4double ySubStepStart[G4FieldTrack::ncompSV 143 G4double ySubStepStart[G4FieldTrack::ncompSVEC];
115 G4FieldTrack yFldTrkStart(y_current); 144 G4FieldTrack yFldTrkStart(y_current);
116 #endif 145 #endif
117 146
118 G4double y[G4FieldTrack::ncompSVEC] = {0., 0 << 147 G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
119 G4double dydx[G4FieldTrack::ncompSVEC] = {0. << 148 G4double ystart[G4FieldTrack::ncompSVEC], yEnd[G4FieldTrack::ncompSVEC];
120 G4double ystart[G4FieldTrack::ncompSVEC] = { <<
121 G4double yEnd[G4FieldTrack::ncompSVEC] = {0. <<
122 G4double x1, x2; 149 G4double x1, x2;
123 G4bool succeeded = true; << 150 G4bool succeeded = true, lastStepSucceeded;
124 151
125 G4double startCurveLength; 152 G4double startCurveLength;
126 153
127 const G4int nvar = fNoVars; << 154 G4int noFullIntegr=0, noSmallIntegr = 0 ;
>> 155 static G4int noGoodSteps =0 ; // Bad = chord > curve-len
>> 156 const G4int nvar= fNoVars;
128 157
129 G4FieldTrack yStartFT(y_current); 158 G4FieldTrack yStartFT(y_current);
130 159
131 // Ensure that hstep > 0 160 // Ensure that hstep > 0
132 // 161 //
133 if( hstep <= 0.0 ) 162 if( hstep <= 0.0 )
134 { 163 {
135 if( hstep == 0.0 ) << 164 if(hstep==0.0)
136 { 165 {
137 std::ostringstream message; << 166 G4cerr << "WARNING - G4MagIntegratorDriver::AccurateAdvance()" << G4endl
138 message << "Proposed step is zero; hstep << 167 << " Proposed step is zero; hstep = " << hstep
139 G4Exception("G4MagInt_Driver::AccurateAd << 168 << " !" << G4endl;
140 "GeomField1001", JustWarning <<
141 return succeeded; 169 return succeeded;
142 } 170 }
143 << 171 else
144 std::ostringstream message; << 172 {
145 message << "Invalid run condition." << G4e << 173 G4cerr << "ERROR - G4MagIntegratorDriver::AccurateAdvance()" << G4endl
146 << "Proposed step is negative; hst << 174 << " Proposed step is negative; hstep = " << hstep
147 << "Requested step cannot be negat << 175 << " !" << G4endl;
148 G4Exception("G4MagInt_Driver::AccurateAdva << 176 G4Exception("G4MagInt_Driver::AccurateAdvance()",
149 "GeomField0003", EventMustBeAb << 177 "InvalidCall", EventMustBeAborted,
150 return false; << 178 "Requested step cannot be negative! Aborting event.");
>> 179 return false;
>> 180 }
151 } 181 }
152 182
153 y_current.DumpToArray( ystart ); 183 y_current.DumpToArray( ystart );
154 184
155 startCurveLength= y_current.GetCurveLength() 185 startCurveLength= y_current.GetCurveLength();
156 x1= startCurveLength; 186 x1= startCurveLength;
157 x2= x1 + hstep; 187 x2= x1 + hstep;
158 188
159 if ( (hinitial > 0.0) && (hinitial < hstep) 189 if ( (hinitial > 0.0) && (hinitial < hstep)
160 && (hinitial > perMillion * hstep) ) 190 && (hinitial > perMillion * hstep) )
161 { 191 {
162 h = hinitial; 192 h = hinitial;
163 } 193 }
164 else // Initial Step size "h" defaults to 194 else // Initial Step size "h" defaults to the full interval
165 { 195 {
166 h = hstep; 196 h = hstep;
167 } 197 }
168 198
169 x = x1; 199 x = x1;
170 200
171 for ( i=0; i<nvar; ++i) { y[i] = ystart[i]; << 201 for (i=0;i<nvar;i++) { y[i] = ystart[i]; }
172 202
173 G4bool lastStep= false; 203 G4bool lastStep= false;
174 nstp = 1; << 204 nstp=1;
175 205
176 do 206 do
177 { 207 {
178 G4ThreeVector StartPos( y[0], y[1], y[2] ) 208 G4ThreeVector StartPos( y[0], y[1], y[2] );
179 209
180 #ifdef G4DEBUG_FIELD 210 #ifdef G4DEBUG_FIELD
181 G4double xSubStepStart= x; 211 G4double xSubStepStart= x;
182 for (i=0; i<nvar; ++i) { ySubStepStart[i] << 212 for (i=0;i<nvar;i++) { ySubStepStart[i] = y[i]; }
183 yFldTrkStart.LoadFromArray(y, fNoIntegrati 213 yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
184 yFldTrkStart.SetCurveLength(x); 214 yFldTrkStart.SetCurveLength(x);
185 #endif 215 #endif
186 216
187 pIntStepper->RightHandSide( y, dydx ); << 217 // Old method - inline call to Equation of Motion
188 ++fNoTotalSteps; << 218 // pIntStepper->RightHandSide( y, dydx );
>> 219 // New method allows to cache field, or state (eg momentum magnitude)
>> 220 pIntStepper->ComputeRightHandSide( y, dydx );
>> 221 fNoTotalSteps++;
189 222
190 // Perform the Integration 223 // Perform the Integration
191 // 224 //
192 if( h > fMinimumStep ) 225 if( h > fMinimumStep )
193 { 226 {
194 OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ; 227 OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
195 //-------------------------------------- 228 //--------------------------------------
>> 229 lastStepSucceeded= (hdid == h);
196 #ifdef G4DEBUG_FIELD 230 #ifdef G4DEBUG_FIELD
197 if (dbg>2) << 231 if (dbg>2) {
198 { << 232 PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp); // Only
199 // PrintStatus( ySubStepStart, xSubSt <<
200 G4DriverReporter::PrintStatus( ySubSte <<
201 } 233 }
202 #endif 234 #endif
203 } 235 }
204 else 236 else
205 { 237 {
206 G4FieldTrack yFldTrk( G4ThreeVector(0,0, 238 G4FieldTrack yFldTrk( G4ThreeVector(0,0,0),
207 G4ThreeVector(0,0, 239 G4ThreeVector(0,0,0), 0., 0., 0., 0. );
208 G4double dchord_step, dyerr, dyerr_len; 240 G4double dchord_step, dyerr, dyerr_len; // What to do with these ?
209 yFldTrk.LoadFromArray(y, fNoIntegrationV 241 yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
210 yFldTrk.SetCurveLength( x ); 242 yFldTrk.SetCurveLength( x );
211 243
212 QuickAdvance( yFldTrk, dydx, h, dchord_s << 244 QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
213 //-------------------------------------- 245 //-----------------------------------------------------
214 246
215 yFldTrk.DumpToArray(y); 247 yFldTrk.DumpToArray(y);
216 248
217 #ifdef G4FLD_STATS 249 #ifdef G4FLD_STATS
218 ++fNoSmallSteps; << 250 fNoSmallSteps++;
219 if ( dyerr_len > fDyerr_max ) { fDyerr_ << 251 if ( dyerr_len > fDyerr_max) { fDyerr_max= dyerr_len; }
220 fDyerrPos_smTot += dyerr_len; 252 fDyerrPos_smTot += dyerr_len;
221 fSumH_sm += h; // Length total for 'sma 253 fSumH_sm += h; // Length total for 'small' steps
222 if (nstp==1) { ++fNoInitialSmallSteps; << 254 if (nstp<=1) { fNoInitialSmallSteps++; }
223 #endif 255 #endif
224 #ifdef G4DEBUG_FIELD 256 #ifdef G4DEBUG_FIELD
225 if (dbg>1) 257 if (dbg>1)
226 { 258 {
227 if(fNoSmallSteps<2) { PrintStatus(ySub 259 if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
228 G4cout << "Another sub-min step, no " 260 G4cout << "Another sub-min step, no " << fNoSmallSteps
229 << " of " << fNoTotalSteps << " 261 << " of " << fNoTotalSteps << " this time " << nstp << G4endl;
230 PrintStatus( ySubStepStart, x1, y, x, 262 PrintStatus( ySubStepStart, x1, y, x, h, nstp); // Only this
231 G4cout << " dyerr= " << dyerr_len << " 263 G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h
232 << " epsilon= " << eps << " hst 264 << " epsilon= " << eps << " hstep= " << hstep
233 << " h= " << h << " hmin= " << 265 << " h= " << h << " hmin= " << fMinimumStep << G4endl;
234 } 266 }
235 #endif 267 #endif
236 if( h == 0.0 ) 268 if( h == 0.0 )
237 { 269 {
238 G4Exception("G4MagInt_Driver::Accurate 270 G4Exception("G4MagInt_Driver::AccurateAdvance()",
239 "GeomField0003", FatalExce << 271 "Integration Step became Zero",
240 "Integration Step became Z << 272 FatalException, "IntegrationStepUnderflow.");
241 } 273 }
242 dyerr = dyerr_len / h; 274 dyerr = dyerr_len / h;
243 hdid = h; << 275 hdid= h;
244 x += hdid; 276 x += hdid;
245 277
246 // Compute suggested new step 278 // Compute suggested new step
247 hnext = ComputeNewStepSize( dyerr/eps, h << 279 hnext= ComputeNewStepSize( dyerr/eps, h);
>> 280
>> 281 // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
>> 282 lastStepSucceeded= (dyerr<= eps);
248 } 283 }
249 284
>> 285 if (lastStepSucceeded) { noFullIntegr++; }
>> 286 else { noSmallIntegr++; }
>> 287
250 G4ThreeVector EndPos( y[0], y[1], y[2] ); 288 G4ThreeVector EndPos( y[0], y[1], y[2] );
251 289
252 #ifdef G4DEBUG_FIELD 290 #ifdef G4DEBUG_FIELD
253 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr)) 291 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
254 { 292 {
255 if( nstp==nStpPr ) { G4cout << "***** M 293 if( nstp==nStpPr ) { G4cout << "***** Many steps ****" << G4endl; }
256 G4cout << "MagIntDrv: " ; 294 G4cout << "MagIntDrv: " ;
257 G4cout << "hdid=" << std::setw(12) << h 295 G4cout << "hdid=" << std::setw(12) << hdid << " "
258 << "hnext=" << std::setw(12) << h 296 << "hnext=" << std::setw(12) << hnext << " "
259 << "hstep=" << std::setw(12) << h << 297 << "hstep=" << std::setw(12) << hstep << " (requested) "
260 << G4endl; << 298 << G4endl;
261 PrintStatus( ystart, x1, y, x, h, (nstp= 299 PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
262 } 300 }
263 #endif 301 #endif
264 302
265 // Check the endpoint 303 // Check the endpoint
266 G4double endPointDist= (EndPos-StartPos).m 304 G4double endPointDist= (EndPos-StartPos).mag();
267 if ( endPointDist >= hdid*(1.+perMillion) 305 if ( endPointDist >= hdid*(1.+perMillion) )
268 { 306 {
269 ++fNoBadSteps; << 307 fNoBadSteps++;
270 308
271 // Issue a warning only for gross differ 309 // Issue a warning only for gross differences -
272 // we understand how small difference oc 310 // we understand how small difference occur.
273 if ( endPointDist >= hdid*(1.+perThousan 311 if ( endPointDist >= hdid*(1.+perThousand) )
274 { 312 {
275 #ifdef G4DEBUG_FIELD 313 #ifdef G4DEBUG_FIELD
276 if (dbg) 314 if (dbg)
277 { 315 {
278 WarnEndPointTooFar ( endPointDist, h 316 WarnEndPointTooFar ( endPointDist, hdid, eps, dbg );
279 G4cerr << " Total steps: bad " << 317 G4cerr << " Total steps: bad " << fNoBadSteps
280 << " current h= " << hdid << << 318 << " good " << noGoodSteps << " current h= " << hdid
>> 319 << G4endl;
281 PrintStatus( ystart, x1, y, x, hstep 320 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
282 } 321 }
283 ++no_warnings; <<
284 #endif 322 #endif
>> 323 no_warnings++;
285 } 324 }
286 } 325 }
>> 326 else
>> 327 {
>> 328 noGoodSteps ++;
>> 329 }
>> 330 // #endif
287 331
288 // Avoid numerous small last steps 332 // Avoid numerous small last steps
289 if( (h < eps * hstep) || (h < fSmallestFra 333 if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
290 { 334 {
291 // No more integration -- the next step 335 // No more integration -- the next step will not happen
292 lastStep = true; 336 lastStep = true;
293 } 337 }
294 else 338 else
295 { 339 {
296 // Check the proposed next stepsize 340 // Check the proposed next stepsize
297 if(std::fabs(hnext) <= Hmin()) 341 if(std::fabs(hnext) <= Hmin())
298 { 342 {
299 #ifdef G4DEBUG_FIELD 343 #ifdef G4DEBUG_FIELD
300 // If simply a very small interval is 344 // If simply a very small interval is being integrated, do not warn
301 if( (x < x2 * (1-eps) ) && // T 345 if( (x < x2 * (1-eps) ) && // The last step can be small: OK
302 (std::fabs(hstep) > Hmin()) ) // a 346 (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
303 { 347 {
304 if(dbg>0) 348 if(dbg>0)
305 { 349 {
306 WarnSmallStepSize( hnext, hstep, h 350 WarnSmallStepSize( hnext, hstep, h, x-x1, nstp );
307 PrintStatus( ystart, x1, y, x, hst 351 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
308 } 352 }
309 ++no_warnings; << 353 no_warnings++;
310 } 354 }
311 #endif 355 #endif
312 // Make sure that the next step is at 356 // Make sure that the next step is at least Hmin.
313 h = Hmin(); 357 h = Hmin();
314 } 358 }
315 else 359 else
316 { 360 {
317 h = hnext; 361 h = hnext;
318 } 362 }
319 363
320 // Ensure that the next step does not o 364 // Ensure that the next step does not overshoot
321 if ( x+h > x2 ) 365 if ( x+h > x2 )
322 { // When stepsize oversh 366 { // When stepsize overshoots, decrease it!
323 h = x2 - x ; // Must cope with diffi 367 h = x2 - x ; // Must cope with difficult rounding-error
324 } // issues if hstep << x 368 } // issues if hstep << x2
325 369
326 if ( h == 0.0 ) 370 if ( h == 0.0 )
327 { 371 {
328 // Cannot progress - accept this as la 372 // Cannot progress - accept this as last step - by default
329 lastStep = true; 373 lastStep = true;
330 #ifdef G4DEBUG_FIELD 374 #ifdef G4DEBUG_FIELD
331 if (dbg>2) 375 if (dbg>2)
332 { 376 {
333 int prec= G4cout.precision(12); 377 int prec= G4cout.precision(12);
334 G4cout << "Warning: G4MagIntegratorD 378 G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
335 << G4endl 379 << G4endl
336 << " Integration step 'h' be 380 << " Integration step 'h' became "
337 << h << " due to roundoff. " 381 << h << " due to roundoff. " << G4endl
338 << " Calculated as difference << 382 << " Calculated as difference of x2= "<< x2 << " and x=" << x
339 << " Forcing termination of 383 << " Forcing termination of advance." << G4endl;
340 G4cout.precision(prec); 384 G4cout.precision(prec);
341 } 385 }
342 #endif 386 #endif
343 } 387 }
344 } 388 }
345 } while ( ((++nstp)<=fMaxNoSteps) && (x < x2 << 389 } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
346 // Loop checking, 07.10.2016, J. Apostolakis <<
347 <<
348 // Have we reached the end ? 390 // Have we reached the end ?
349 // --> a better test might be x-x2 > an_e 391 // --> a better test might be x-x2 > an_epsilon
350 392
351 succeeded = (x>=x2); // If it was a "forced << 393 succeeded= (x>=x2); // If it was a "forced" last step
352 394
353 for (i=0; i<nvar; ++i) { yEnd[i] = y[i]; } << 395 for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
354 396
355 // Put back the values. 397 // Put back the values.
356 y_current.LoadFromArray( yEnd, fNoIntegratio 398 y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
357 y_current.SetCurveLength( x ); 399 y_current.SetCurveLength( x );
358 400
359 if(nstp > fMaxNoSteps) 401 if(nstp > fMaxNoSteps)
360 { 402 {
>> 403 no_warnings++;
361 succeeded = false; 404 succeeded = false;
362 #ifdef G4DEBUG_FIELD 405 #ifdef G4DEBUG_FIELD
363 ++no_warnings; <<
364 if (dbg) 406 if (dbg)
365 { 407 {
366 WarnTooManyStep( x1, x2, x ); // Issue 408 WarnTooManyStep( x1, x2, x ); // Issue WARNING
367 PrintStatus( yEnd, x1, y, x, hstep, -nst 409 PrintStatus( yEnd, x1, y, x, hstep, -nstp);
368 } 410 }
369 #endif 411 #endif
370 } 412 }
371 413
372 #ifdef G4DEBUG_FIELD 414 #ifdef G4DEBUG_FIELD
373 if( dbg && no_warnings ) 415 if( dbg && no_warnings )
374 { 416 {
375 G4cerr << "G4MagIntegratorDriver exit stat << 417 G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
376 PrintStatus( yEnd, x1, y, x, hstep, nstp); 418 PrintStatus( yEnd, x1, y, x, hstep, nstp);
377 } 419 }
378 #endif 420 #endif
379 421
380 return succeeded; 422 return succeeded;
381 } // end of AccurateAdvance ................. 423 } // end of AccurateAdvance ...........................
382 424
383 // ------------------------------------------- 425 // ---------------------------------------------------------
384 426
385 void 427 void
386 G4MagInt_Driver::WarnSmallStepSize( G4double h 428 G4MagInt_Driver::WarnSmallStepSize( G4double hnext, G4double hstep,
387 G4double h 429 G4double h, G4double xDone,
388 G4int nstp 430 G4int nstp)
389 { 431 {
390 static G4ThreadLocal G4int noWarningsIssued << 432 static G4int noWarningsIssued =0;
391 const G4int maxNoWarnings = 10; // Number << 433 const G4int maxNoWarnings = 10; // Number of verbose warnings
392 std::ostringstream message; <<
393 if( (noWarningsIssued < maxNoWarnings) || fV 434 if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
394 { 435 {
395 message << "The stepsize for the next iter << 436 G4cerr << " Warning (G4MagIntegratorDriver::AccurateAdvance):" << G4endl
396 << ", is too small - in Step numbe << 437 << " The stepsize for the next iteration, " << hnext
397 << "The minimum for the driver is << 438 << ", is too small - in Step number " << nstp << "." << G4endl;
398 << "Requested integr. length was " << 439 G4cerr << " The minimum for the driver is " << Hmin() << G4endl;
399 << "The size of this sub-step was << 440 G4cerr << " Requested integr. length was " << hstep << " ." << G4endl;
400 << "The integrations has already g << 441 G4cerr << " The size of this sub-step was " << h << " ." << G4endl;
>> 442 G4cerr << " The integrations has already gone " << xDone << G4endl;
401 } 443 }
402 else 444 else
403 { 445 {
404 message << "Too small 'next' step " << hne << 446 G4cerr << " G4MagInt_Driver: Too small 'next' step " << hnext
405 << ", step-no: " << nstp << G4endl << 447 << " step-no " << nstp << G4endl;
406 << ", this sub-step: " << h << 448 G4cerr << " this sub-step " << h
407 << ", req_tot_len: " << hstep << 449 << " req_tot_len " << hstep
408 << ", done: " << xDone << ", min: << 450 << " done " << xDone
409 } << 451 << " min " << Hmin()
410 G4Exception("G4MagInt_Driver::WarnSmallStepS << 452 << G4endl;
411 JustWarning, message); << 453 }
412 ++noWarningsIssued; << 454 noWarningsIssued++;
413 } 455 }
414 456
415 // ------------------------------------------- 457 // ---------------------------------------------------------
416 458
417 void 459 void
418 G4MagInt_Driver::WarnTooManyStep( G4double x1s 460 G4MagInt_Driver::WarnTooManyStep( G4double x1start,
419 G4double x2e 461 G4double x2end,
420 G4double xCu << 462 G4double xCurrent)
421 { 463 {
422 std::ostringstream message; << 464 G4cerr << " Warning (G4MagIntegratorDriver):" << G4endl
423 message << "The number of steps used in the << 465 << " The number of steps used in the Integration driver"
424 << " (Runge-Kutta) is too many." << << 466 << " (Runge-Kutta) is too many." << G4endl;
425 << "Integration of the interval was << 467 G4cerr << " Integration of the interval was not completed !" << G4endl
426 << "Only a " << (xCurrent-x1start)* << 468 << " Only a " << (xCurrent-x1start)*100/(x2end-x1start)
427 << " % fraction of it was done."; << 469 <<" % fraction of it was done." << G4endl;
428 G4Exception("G4MagInt_Driver::WarnTooManySt <<
429 JustWarning, message); <<
430 } 470 }
431 471
432 // ------------------------------------------- 472 // ---------------------------------------------------------
433 473
434 void 474 void
435 G4MagInt_Driver::WarnEndPointTooFar (G4double 475 G4MagInt_Driver::WarnEndPointTooFar (G4double endPointDist,
436 G4double << 476 G4double h ,
437 G4double << 477 G4double eps,
438 G4int << 478 G4int dbg)
439 { 479 {
440 static G4ThreadLocal G4double maxRelError = << 480 static G4double maxRelError=0.0;
441 G4bool isNewMax, prNewMax; 481 G4bool isNewMax, prNewMax;
442 482
443 isNewMax = endPointDist > (1.0 + maxRelError 483 isNewMax = endPointDist > (1.0 + maxRelError) * h;
444 prNewMax = endPointDist > (1.0 + 1.05 * maxR 484 prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
445 if( isNewMax ) { maxRelError= endPointDist / 485 if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
446 486
447 if( (dbg != 0) && (h > G4GeometryTolerance:: << 487 if( dbg && (h > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
448 && ( (dbg>1) || prNewMax || (endPoin 488 && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
449 { 489 {
450 static G4ThreadLocal G4int noWarnings = 0; << 490 static G4int noWarnings = 0;
451 std::ostringstream message; << 491 if( (noWarnings ++ < 10) || (dbg>2) )
452 if( (noWarnings++ < 10) || (dbg>2) ) <<
453 { 492 {
454 message << "The integration produced an << 493 G4cerr << " Warning (G4MagIntegratorDriver): " << G4endl
455 << "is further from the start-po << 494 << " The integration produced an endpoint which " << G4endl
456 << G4endl; << 495 << " is further from the startpoint than the curve length."
>> 496 << G4endl;
>> 497
>> 498 G4cerr << " Distance of endpoints = " << endPointDist
>> 499 << " curve length = " << h
>> 500 << " Difference (curveLen-endpDist)= " << (h - endPointDist)
>> 501 << " relative = " << (h-endPointDist) / h
>> 502 << " epsilon = " << eps
>> 503 << G4endl;
>> 504 }
>> 505 else
>> 506 {
>> 507 G4cerr << " Warning:"
>> 508 << " dist_e= " << endPointDist
>> 509 << " h_step = " << h
>> 510 << " Diff (hs-ed)= " << (h - endPointDist)
>> 511 << " rel = " << (h-endPointDist) / h
>> 512 << " eps = " << eps
>> 513 << " (from G4MagInt_Driver)" << G4endl;
457 } 514 }
458 message << " Distance of endpoints = " << <<
459 << ", curve length = " << h << G4e <<
460 << " Difference (curveLen-endpDis <<
461 << ", relative = " << (h-endPointD <<
462 << ", epsilon = " << eps; <<
463 G4Exception("G4MagInt_Driver::WarnEndPoint <<
464 JustWarning, message); <<
465 } 515 }
466 } 516 }
467 517
468 // ------------------------------------------- 518 // ---------------------------------------------------------
469 519
470 void 520 void
471 G4MagInt_Driver::OneGoodStep( G4double y[ 521 G4MagInt_Driver::OneGoodStep( G4double y[], // InOut
472 const G4double dy 522 const G4double dydx[],
473 G4double& x 523 G4double& x, // InOut
474 G4double ht 524 G4double htry,
475 G4double ep 525 G4double eps_rel_max,
476 G4double& h 526 G4double& hdid, // Out
477 G4double& h 527 G4double& hnext ) // Out
478 528
479 // Driver for one Runge-Kutta Step with monito 529 // Driver for one Runge-Kutta Step with monitoring of local truncation error
480 // to ensure accuracy and adjust stepsize. Inp 530 // to ensure accuracy and adjust stepsize. Input are dependent variable
481 // array y[0,...,5] and its derivative dydx[0, 531 // array y[0,...,5] and its derivative dydx[0,...,5] at the
482 // starting value of the independent variable 532 // starting value of the independent variable x . Also input are stepsize
483 // to be attempted htry, and the required accu 533 // to be attempted htry, and the required accuracy eps. On output y and x
484 // are replaced by their new values, hdid is t 534 // are replaced by their new values, hdid is the stepsize that was actually
485 // accomplished, and hnext is the estimated ne 535 // accomplished, and hnext is the estimated next stepsize.
486 // This is similar to the function rkqs from t 536 // This is similar to the function rkqs from the book:
487 // Numerical Recipes in C: The Art of Scientif 537 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
488 // Edition, by William H. Press, Saul A. Teuko 538 // Edition, by William H. Press, Saul A. Teukolsky, William T.
489 // Vetterling, and Brian P. Flannery (Cambridg 539 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
490 // 16.2 Adaptive StepSize Control for Runge-Ku 540 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
491 541
492 { 542 {
493 G4double errmax_sq; 543 G4double errmax_sq;
494 G4double h, htemp, xnew ; 544 G4double h, htemp, xnew ;
495 545
496 G4double yerr[G4FieldTrack::ncompSVEC], ytem 546 G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
497 547
498 h = htry ; // Set stepsize to the initial tr 548 h = htry ; // Set stepsize to the initial trial value
499 549
500 G4double inv_eps_vel_sq = 1.0 / (eps_rel_max 550 G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
501 551
502 G4double errpos_sq = 0.0; // square of di << 552 G4double errpos_sq=0.0; // square of displacement error
503 G4double errvel_sq = 0.0; // square of mo << 553 G4double errvel_sq=0.0; // square of momentum vector difference
504 G4double errspin_sq = 0.0; // square of sp << 554 G4double errspin_sq=0.0; // square of spin vector difference
>> 555
>> 556 G4int iter;
505 557
>> 558 static G4int tot_no_trials=0;
506 const G4int max_trials=100; 559 const G4int max_trials=100;
507 560
508 G4ThreeVector Spin(y[9],y[10],y[11]); 561 G4ThreeVector Spin(y[9],y[10],y[11]);
509 G4double spin_mag2 = Spin.mag2(); << 562 G4bool hasSpin= (Spin.mag2() > 0.0);
510 G4bool hasSpin = (spin_mag2 > 0.0); <<
511 563
512 for (G4int iter=0; iter<max_trials; ++iter) << 564 for (iter=0; iter<max_trials ;iter++)
513 { 565 {
>> 566 tot_no_trials++;
514 pIntStepper-> Stepper(y,dydx,h,ytemp,yerr) 567 pIntStepper-> Stepper(y,dydx,h,ytemp,yerr);
515 // ******* 568 // *******
516 G4double eps_pos = eps_rel_max * std::max( 569 G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
517 G4double inv_eps_pos_sq = 1.0 / (eps_pos*e 570 G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
518 571
519 // Evaluate accuracy 572 // Evaluate accuracy
520 // 573 //
521 errpos_sq = sqr(yerr[0]) + sqr(yerr[1]) + 574 errpos_sq = sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
522 errpos_sq *= inv_eps_pos_sq; // Scale rela 575 errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
523 576
524 // Accuracy for momentum 577 // Accuracy for momentum
525 G4double magvel_sq= sqr(y[3]) + sqr(y[4]) << 578 errvel_sq = (sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) )
526 G4double sumerr_sq = sqr(yerr[3]) + sqr(y << 579 / (sqr(y[3]) + sqr(y[4]) + sqr(y[5]) );
527 if( magvel_sq > 0.0 ) <<
528 { <<
529 errvel_sq = sumerr_sq / magvel_sq; <<
530 } <<
531 else <<
532 { <<
533 std::ostringstream message; <<
534 message << "Found case of zero momentum <<
535 << "- iteration= " << iter << " <<
536 G4Exception("G4MagInt_Driver::OneGoodSt <<
537 "GeomField1001", JustWarnin <<
538 errvel_sq = sumerr_sq; <<
539 } <<
540 errvel_sq *= inv_eps_vel_sq; 580 errvel_sq *= inv_eps_vel_sq;
541 errmax_sq = std::max( errpos_sq, errvel_sq 581 errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
542 582
543 if( hasSpin ) 583 if( hasSpin )
544 { 584 {
545 // Accuracy for spin 585 // Accuracy for spin
546 errspin_sq = ( sqr(yerr[9]) + sqr(yerr[ 586 errspin_sq = ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
547 / spin_mag2; // ( sqr(y[9 << 587 / ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
548 errspin_sq *= inv_eps_vel_sq; 588 errspin_sq *= inv_eps_vel_sq;
549 errmax_sq = std::max( errmax_sq, errspin 589 errmax_sq = std::max( errmax_sq, errspin_sq );
550 } << 590 }
551 591
552 if ( errmax_sq <= 1.0 ) { break; } // Ste 592 if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.
553 593
554 // Step failed; compute the size of retria 594 // Step failed; compute the size of retrial Step.
555 htemp = GetSafety() * h * std::pow( errmax << 595 htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
556 596
557 if (htemp >= 0.1*h) { h = htemp; } // Tr 597 if (htemp >= 0.1*h) { h = htemp; } // Truncation error too large,
558 else { h = 0.1*h; } // re 598 else { h = 0.1*h; } // reduce stepsize, but no more
559 // th 599 // than a factor of 10
560 xnew = x + h; 600 xnew = x + h;
561 if(xnew == x) 601 if(xnew == x)
562 { 602 {
563 std::ostringstream message; << 603 G4cerr << "G4MagIntegratorDriver::OneGoodStep:" << G4endl
564 message << "Stepsize underflow in Steppe << 604 << " Stepsize underflow in Stepper " << G4endl ;
565 << "- Step's start x=" << x << " << 605 G4cerr << " Step's start x=" << x << " and end x= " << xnew
566 << " are equal !! " << G4endl << 606 << " are equal !! " << G4endl
567 << " Due to step-size= " << h << 607 <<" Due to step-size= " << h
568 << ". Note that input step was " << 608 << " . Note that input step was " << htry << G4endl;
569 G4Exception("G4MagInt_Driver::OneGoodSte <<
570 "GeomField1001", JustWarning <<
571 break; 609 break;
572 } 610 }
573 } 611 }
574 612
>> 613 #ifdef G4FLD_STATS
>> 614 // Sum of squares of position error // and momentum dir (underestimated)
>> 615 fSumH_lg += h;
>> 616 fDyerrPos_lgTot += errpos_sq;
>> 617 fDyerrVel_lgTot += errvel_sq * h * h;
>> 618 #endif
>> 619
575 // Compute size of next Step 620 // Compute size of next Step
576 if (errmax_sq > errcon*errcon) 621 if (errmax_sq > errcon*errcon)
577 { 622 {
578 hnext = GetSafety()*h*std::pow(errmax_sq, 623 hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
579 } 624 }
580 else 625 else
581 { 626 {
582 hnext = max_stepping_increase*h ; // No mo 627 hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
583 } 628 }
584 x += (hdid = h); 629 x += (hdid = h);
585 630
586 for(G4int k=0; k<fNoIntegrationVariables; ++ << 631 for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
587 632
588 return; 633 return;
589 } << 634 } // end of OneGoodStep .............................
590 635
591 //-------------------------------------------- 636 //----------------------------------------------------------------------
592 637
593 // QuickAdvance just tries one Step - it does 638 // QuickAdvance just tries one Step - it does not ensure accuracy
594 // 639 //
595 G4bool G4MagInt_Driver::QuickAdvance(G4FieldTr << 640 G4bool G4MagInt_Driver::QuickAdvance(
596 const G4double << 641 G4FieldTrack& y_posvel, // INOUT
597 G4double << 642 const G4double dydx[],
598 G4double& << 643 G4double hstep, // In
599 G4double& << 644 G4double& dchord_step,
600 G4double& << 645 G4double& dyerr_pos_sq,
>> 646 G4double& dyerr_mom_rel_sq )
601 { 647 {
602 G4Exception("G4MagInt_Driver::QuickAdvance() << 648 G4Exception("G4MagInt_Driver::QuickAdvance()", "NotImplemented",
603 FatalException, "Not yet impleme 649 FatalException, "Not yet implemented.");
604 650
605 // Use the parameters of this method, to ple 651 // Use the parameters of this method, to please compiler
606 // <<
607 dchord_step = dyerr_pos_sq = hstep * hstep * 652 dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
608 dyerr_mom_rel_sq = y_posvel.GetPosition().ma 653 dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
609 return true; 654 return true;
610 } 655 }
611 656
612 //-------------------------------------------- 657 //----------------------------------------------------------------------
613 658
614 G4bool G4MagInt_Driver::QuickAdvance(G4FieldTr << 659 G4bool G4MagInt_Driver::QuickAdvance(
615 const G4double << 660 G4FieldTrack& y_posvel, // INOUT
616 G4double << 661 const G4double dydx[],
617 G4double& << 662 G4double hstep, // In
618 G4double& << 663 G4double& dchord_step,
>> 664 G4double& dyerr )
619 { 665 {
620 G4double dyerr_pos_sq, dyerr_mom_rel_sq; 666 G4double dyerr_pos_sq, dyerr_mom_rel_sq;
621 G4double yerr_vec[G4FieldTrack::ncompSVEC], 667 G4double yerr_vec[G4FieldTrack::ncompSVEC],
622 yarrin[G4FieldTrack::ncompSVEC], ya 668 yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC];
623 G4double s_start; 669 G4double s_start;
624 G4double dyerr_mom_sq, vel_mag_sq, inv_vel_m 670 G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
625 671
>> 672 static G4int no_call=0;
>> 673 no_call ++;
>> 674
626 // Move data into array 675 // Move data into array
627 y_posvel.DumpToArray( yarrin ); // yar 676 y_posvel.DumpToArray( yarrin ); // yarrin <== y_posvel
628 s_start = y_posvel.GetCurveLength(); 677 s_start = y_posvel.GetCurveLength();
629 678
630 // Do an Integration Step 679 // Do an Integration Step
631 pIntStepper-> Stepper(yarrin, dydx, hstep, y 680 pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
>> 681 // *******
632 682
633 // Estimate curve-chord distance 683 // Estimate curve-chord distance
634 dchord_step= pIntStepper-> DistChord(); 684 dchord_step= pIntStepper-> DistChord();
>> 685 // *********
635 686
636 // Put back the values. yarrout ==> y_posve 687 // Put back the values. yarrout ==> y_posvel
637 y_posvel.LoadFromArray( yarrout, fNoIntegrat 688 y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
638 y_posvel.SetCurveLength( s_start + hstep ); 689 y_posvel.SetCurveLength( s_start + hstep );
639 690
640 #ifdef G4DEBUG_FIELD 691 #ifdef G4DEBUG_FIELD
641 if(fVerboseLevel>2) 692 if(fVerboseLevel>2)
642 { 693 {
643 G4cout << "G4MagIntDrv: Quick Advance" << 694 G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
644 PrintStatus( yarrin, s_start, yarrout, s_s 695 PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
645 } 696 }
646 #endif 697 #endif
647 698
648 // A single measure of the error 699 // A single measure of the error
649 // TO-DO : account for energy, spin, 700 // TO-DO : account for energy, spin, ... ?
650 vel_mag_sq = ( sqr(yarrout[3])+sqr(yarrout 701 vel_mag_sq = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
651 inv_vel_mag_sq = 1.0 / vel_mag_sq; 702 inv_vel_mag_sq = 1.0 / vel_mag_sq;
652 dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_v 703 dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
653 dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_v 704 dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
654 dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_ma 705 dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
655 706
656 // Calculate also the change in the momentum 707 // Calculate also the change in the momentum squared also ???
657 // G4double veloc_square = y_posvel.GetVeloc 708 // G4double veloc_square = y_posvel.GetVelocity().mag2();
658 // ... 709 // ...
659 710
660 #ifdef RETURN_A_NEW_STEP_LENGTH 711 #ifdef RETURN_A_NEW_STEP_LENGTH
661 // The following step cannot be done here be 712 // The following step cannot be done here because "eps" is not known.
662 dyerr_len = std::sqrt( dyerr_len_sq ); 713 dyerr_len = std::sqrt( dyerr_len_sq );
663 dyerr_len_sq /= eps ; 714 dyerr_len_sq /= eps ;
664 715
665 // Look at the velocity deviation ? 716 // Look at the velocity deviation ?
666 // sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(ye 717 // sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
667 718
668 // Set suggested new step 719 // Set suggested new step
669 hstep = ComputeNewStepSize( dyerr_len, hstep << 720 hstep= ComputeNewStepSize( dyerr_len, hstep);
670 #endif 721 #endif
671 722
672 if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr( 723 if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
673 { 724 {
674 dyerr = std::sqrt(dyerr_pos_sq); 725 dyerr = std::sqrt(dyerr_pos_sq);
675 } 726 }
676 else 727 else
677 { 728 {
678 // Scale it to the current step size - for 729 // Scale it to the current step size - for now
679 dyerr = std::sqrt(dyerr_mom_rel_sq) * hste 730 dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
680 } 731 }
681 732
682 return true; 733 return true;
683 } 734 }
684 735
685 // ------------------------------------------- 736 // --------------------------------------------------------------------------
686 737
687 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT 738 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
688 G4bool G4MagInt_Driver::QuickAdvance(G4double << 739 G4bool G4MagInt_Driver::QuickAdvance(
689 const G4double << 740 G4double yarrin[], // In
690 G4double << 741 const G4double dydx[],
691 G4double << 742 G4double hstep, // In
692 G4double << 743 G4double yarrout[],
693 G4double << 744 G4double& dchord_step,
>> 745 G4double& dyerr ) // In length
694 { 746 {
695 G4Exception("G4MagInt_Driver::QuickAdvance() << 747 G4Exception("G4MagInt_Driver::QuickAdvance()", "NotImplemented",
696 FatalException, "Not yet impleme 748 FatalException, "Not yet implemented.");
697 dyerr = dchord_step = hstep * yarrin[0] * dy 749 dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
698 yarrout[0]= yarrin[0]; 750 yarrout[0]= yarrin[0];
699 } 751 }
700 #endif 752 #endif
701 753
702 // ------------------------------------------- 754 // --------------------------------------------------------------------------
703 755
704 // This method computes new step sizes - but d << 756 // This method computes new step sizes - but does not limit changes to
705 // within certain factors << 757 // within certain factors
706 // 758 //
707 G4double G4MagInt_Driver:: << 759 G4double
708 ComputeNewStepSize_WithoutReductionLimit(G4dou << 760 G4MagInt_Driver::ComputeNewStepSize(
709 G4double hstepCurrent) // << 761 G4double errMaxNorm, // max error (normalised)
>> 762 G4double hstepCurrent) // current step size
710 { 763 {
711 G4double hnew; 764 G4double hnew;
712 765
713 // Compute size of next Step for a failed st 766 // Compute size of next Step for a failed step
714 if(errMaxNorm > 1.0 ) 767 if(errMaxNorm > 1.0 )
715 { 768 {
716 // Step failed; compute the size of retria 769 // Step failed; compute the size of retrial Step.
717 hnew = GetSafety()*hstepCurrent*std::pow(e 770 hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
718 } << 771 } else if(errMaxNorm > 0.0 ) {
719 else if(errMaxNorm > 0.0 ) <<
720 { <<
721 // Compute size of next Step for a success 772 // Compute size of next Step for a successful step
722 hnew = GetSafety()*hstepCurrent*std::pow(e 773 hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
723 } << 774 } else {
724 else <<
725 { <<
726 // if error estimate is zero (possible) or 775 // if error estimate is zero (possible) or negative (dubious)
727 hnew = max_stepping_increase * hstepCurren 776 hnew = max_stepping_increase * hstepCurrent;
728 } 777 }
729 778
730 return hnew; 779 return hnew;
731 } 780 }
732 781
733 // ------------------------------------------- 782 // ---------------------------------------------------------------------------
734 783
735 G4double <<
736 G4MagInt_Driver::ComputeNewStepSize( <<
737 G4double errMaxNorm <<
738 G4double hstepCurre <<
739 { <<
740 // Legacy behaviour: <<
741 return ComputeNewStepSize_WithoutReductionL <<
742 // 'Improved' behaviour - at least more con <<
743 // return ComputeNewStepSize_WithinLimits( <<
744 } <<
745 <<
746 // This method computes new step sizes limitin 784 // This method computes new step sizes limiting changes within certain factors
747 // 785 //
748 // It shares its logic with AccurateAdvance. 786 // It shares its logic with AccurateAdvance.
749 // They are kept separate currently for optimi 787 // They are kept separate currently for optimisation.
750 // 788 //
751 G4double 789 G4double
752 G4MagInt_Driver::ComputeNewStepSize_WithinLimi 790 G4MagInt_Driver::ComputeNewStepSize_WithinLimits(
753 G4double errMaxNorm 791 G4double errMaxNorm, // max error (normalised)
754 G4double hstepCurre 792 G4double hstepCurrent) // current step size
755 { 793 {
756 G4double hnew; 794 G4double hnew;
757 795
758 // Compute size of next Step for a failed st 796 // Compute size of next Step for a failed step
759 if (errMaxNorm > 1.0 ) 797 if (errMaxNorm > 1.0 )
760 { 798 {
761 // Step failed; compute the size of retria 799 // Step failed; compute the size of retrial Step.
762 hnew = GetSafety()*hstepCurrent*std::pow(e 800 hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
763 801
764 if (hnew < max_stepping_decrease*hstepCurr 802 if (hnew < max_stepping_decrease*hstepCurrent)
765 { 803 {
766 hnew = max_stepping_decrease*hstepCurren 804 hnew = max_stepping_decrease*hstepCurrent ;
767 // reduce stepsize, b 805 // reduce stepsize, but no more
768 // than this factor ( 806 // than this factor (value= 1/10)
769 } 807 }
770 } 808 }
771 else 809 else
772 { 810 {
773 // Compute size of next Step for a success 811 // Compute size of next Step for a successful step
774 if (errMaxNorm > errcon) 812 if (errMaxNorm > errcon)
775 { hnew = GetSafety()*hstepCurrent*std::po 813 { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
776 else // No more than a factor of 5 increa 814 else // No more than a factor of 5 increase
777 { hnew = max_stepping_increase * hstepCur 815 { hnew = max_stepping_increase * hstepCurrent; }
778 } 816 }
779 return hnew; 817 return hnew;
780 } 818 }
781 819
782 // ------------------------------------------- 820 // ---------------------------------------------------------------------------
783 821
784 void G4MagInt_Driver::PrintStatus( const G4dou << 822 void G4MagInt_Driver::PrintStatus( const G4double* StartArr,
785 G4dou << 823 G4double xstart,
786 const G4dou << 824 const G4double* CurrentArr,
787 G4dou << 825 G4double xcurrent,
788 G4dou << 826 G4double requestStep,
789 G4int << 827 G4int subStepNo)
790 // Potentially add as arguments: 828 // Potentially add as arguments:
791 // <dydx> 829 // <dydx> - as Initial Force
792 // stepTaken 830 // stepTaken(hdid) - last step taken
793 // nextStep 831 // nextStep (hnext) - proposal for size
794 { 832 {
795 G4FieldTrack StartFT(G4ThreeVector(0,0,0), 833 G4FieldTrack StartFT(G4ThreeVector(0,0,0),
796 G4ThreeVector(0,0,0), 0., 0., 834 G4ThreeVector(0,0,0), 0., 0., 0., 0. );
797 G4FieldTrack CurrentFT (StartFT); 835 G4FieldTrack CurrentFT (StartFT);
798 836
799 StartFT.LoadFromArray( StartArr, fNoIntegra 837 StartFT.LoadFromArray( StartArr, fNoIntegrationVariables);
800 StartFT.SetCurveLength( xstart); 838 StartFT.SetCurveLength( xstart);
801 CurrentFT.LoadFromArray( CurrentArr, fNoInt 839 CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables);
802 CurrentFT.SetCurveLength( xcurrent ); 840 CurrentFT.SetCurveLength( xcurrent );
803 841
804 PrintStatus(StartFT, CurrentFT, requestStep 842 PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
805 } 843 }
806 844
807 // ------------------------------------------- 845 // ---------------------------------------------------------------------------
808 846
809 void G4MagInt_Driver::PrintStatus(const G4Fiel << 847 void G4MagInt_Driver::PrintStatus(
810 const G4Fiel << 848 const G4FieldTrack& StartFT,
811 G4doub << 849 const G4FieldTrack& CurrentFT,
812 G4int << 850 G4double requestStep,
>> 851 G4int subStepNo)
813 { 852 {
814 G4int verboseLevel= fVerboseLevel; 853 G4int verboseLevel= fVerboseLevel;
815 const G4int noPrecision = 5; << 854 static G4int noPrecision= 5;
816 G4long oldPrec= G4cout.precision(noPrecisi << 855 G4int oldPrec= G4cout.precision(noPrecision);
817 // G4cout.setf(ios_base::fixed,ios_base::f 856 // G4cout.setf(ios_base::fixed,ios_base::floatfield);
818 857
819 const G4ThreeVector StartPosition= S 858 const G4ThreeVector StartPosition= StartFT.GetPosition();
820 const G4ThreeVector StartUnitVelocity= S 859 const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir();
821 const G4ThreeVector CurrentPosition= C 860 const G4ThreeVector CurrentPosition= CurrentFT.GetPosition();
822 const G4ThreeVector CurrentUnitVelocity= C 861 const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
823 862
824 G4double DotStartCurrentVeloc= StartUnitV 863 G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
825 864
826 G4double step_len= CurrentFT.GetCurveLengt 865 G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
827 G4double subStepSize = step_len; 866 G4double subStepSize = step_len;
828 867
829 if( (subStepNo <= 1) || (verboseLevel > 3) 868 if( (subStepNo <= 1) || (verboseLevel > 3) )
830 { 869 {
831 subStepNo = - subStepNo; // To a 870 subStepNo = - subStepNo; // To allow printing banner
832 871
833 G4cout << std::setw( 6) << " " << std: 872 G4cout << std::setw( 6) << " " << std::setw( 25)
834 << " G4MagInt_Driver: Current Po 873 << " G4MagInt_Driver: Current Position and Direction" << " "
835 << G4endl; 874 << G4endl;
836 G4cout << std::setw( 5) << "Step#" << " 875 G4cout << std::setw( 5) << "Step#" << " "
837 << std::setw( 7) << "s-curve" << 876 << std::setw( 7) << "s-curve" << " "
838 << std::setw( 9) << "X(mm)" << " 877 << std::setw( 9) << "X(mm)" << " "
839 << std::setw( 9) << "Y(mm)" << " 878 << std::setw( 9) << "Y(mm)" << " "
840 << std::setw( 9) << "Z(mm)" << " 879 << std::setw( 9) << "Z(mm)" << " "
841 << std::setw( 8) << " N_x " << " 880 << std::setw( 8) << " N_x " << " "
842 << std::setw( 8) << " N_y " << " 881 << std::setw( 8) << " N_y " << " "
843 << std::setw( 8) << " N_z " << " 882 << std::setw( 8) << " N_z " << " "
844 << std::setw( 8) << " N^2-1 " << 883 << std::setw( 8) << " N^2-1 " << " "
845 << std::setw(10) << " N(0).N " < 884 << std::setw(10) << " N(0).N " << " "
846 << std::setw( 7) << "KinEner " < 885 << std::setw( 7) << "KinEner " << " "
847 << std::setw(12) << "Track-l" << 886 << std::setw(12) << "Track-l" << " " // Add the Sub-step ??
848 << std::setw(12) << "Step-len" < 887 << std::setw(12) << "Step-len" << " "
849 << std::setw(12) << "Step-len" < 888 << std::setw(12) << "Step-len" << " "
850 << std::setw( 9) << "ReqStep" << 889 << std::setw( 9) << "ReqStep" << " "
851 << G4endl; 890 << G4endl;
852 } 891 }
853 892
854 if( (subStepNo <= 0) ) 893 if( (subStepNo <= 0) )
855 { 894 {
856 PrintStat_Aux( StartFT, requestStep, 0. 895 PrintStat_Aux( StartFT, requestStep, 0.,
857 0, 0.0, 896 0, 0.0, 1.0);
>> 897 //*************
858 } 898 }
859 899
860 if( verboseLevel <= 3 ) 900 if( verboseLevel <= 3 )
861 { 901 {
862 G4cout.precision(noPrecision); 902 G4cout.precision(noPrecision);
863 PrintStat_Aux( CurrentFT, requestStep, s 903 PrintStat_Aux( CurrentFT, requestStep, step_len,
864 subStepNo, subStepSize, D 904 subStepNo, subStepSize, DotStartCurrentVeloc );
>> 905 //*************
865 } 906 }
866 907
>> 908 else // if( verboseLevel > 3 )
>> 909 {
>> 910 // Multi-line output
>> 911
>> 912 // G4cout << "Current Position is " << CurrentPosition << G4endl
>> 913 // << " and UnitVelocity is " << CurrentUnitVelocity << G4endl;
>> 914 // G4cout << "Step taken was " << step_len
>> 915 // << " out of PhysicalStep= " << requestStep << G4endl;
>> 916 // G4cout << "Final safety is: " << safety << G4endl;
>> 917 // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
>> 918 // << G4endl << G4endl;
>> 919 }
867 G4cout.precision(oldPrec); 920 G4cout.precision(oldPrec);
868 } 921 }
869 922
870 // ------------------------------------------- 923 // ---------------------------------------------------------------------------
871 924
872 void G4MagInt_Driver::PrintStat_Aux(const G4Fi << 925 void G4MagInt_Driver::PrintStat_Aux(
873 G4do << 926 const G4FieldTrack& aFieldTrack,
874 G4do << 927 G4double requestStep,
875 G4in << 928 G4double step_len,
876 G4do << 929 G4int subStepNo,
877 G4do << 930 G4double subStepSize,
>> 931 G4double dotVeloc_StartCurr)
878 { 932 {
879 const G4ThreeVector Position = aFieldTrack << 933 const G4ThreeVector Position= aFieldTrack.GetPosition();
880 const G4ThreeVector UnitVelocity = aFieldT << 934 const G4ThreeVector UnitVelocity= aFieldTrack.GetMomentumDir();
881 935
882 if( subStepNo >= 0) 936 if( subStepNo >= 0)
883 { 937 {
884 G4cout << std::setw( 5) << subStepNo << 938 G4cout << std::setw( 5) << subStepNo << " ";
885 } 939 }
886 else 940 else
887 { 941 {
888 G4cout << std::setw( 5) << "Start" << " 942 G4cout << std::setw( 5) << "Start" << " ";
889 } 943 }
890 G4double curveLen= aFieldTrack.GetCurveLen 944 G4double curveLen= aFieldTrack.GetCurveLength();
891 G4cout << std::setw( 7) << curveLen; 945 G4cout << std::setw( 7) << curveLen;
892 G4cout << std::setw( 9) << Position.x() << 946 G4cout << std::setw( 9) << Position.x() << " "
893 << std::setw( 9) << Position.y() << 947 << std::setw( 9) << Position.y() << " "
894 << std::setw( 9) << Position.z() << 948 << std::setw( 9) << Position.z() << " "
895 << std::setw( 8) << UnitVelocity.x( 949 << std::setw( 8) << UnitVelocity.x() << " "
896 << std::setw( 8) << UnitVelocity.y( 950 << std::setw( 8) << UnitVelocity.y() << " "
897 << std::setw( 8) << UnitVelocity.z( 951 << std::setw( 8) << UnitVelocity.z() << " ";
898 G4long oldprec= G4cout.precision(3); << 952 G4int oldprec= G4cout.precision(3);
899 G4cout << std::setw( 8) << UnitVelocity.ma 953 G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
900 G4cout.precision(6); 954 G4cout.precision(6);
901 G4cout << std::setw(10) << dotVeloc_StartC 955 G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
902 G4cout.precision(oldprec); 956 G4cout.precision(oldprec);
903 G4cout << std::setw( 7) << aFieldTrack.Get 957 G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
904 G4cout << std::setw(12) << step_len << " " 958 G4cout << std::setw(12) << step_len << " ";
905 959
906 static G4ThreadLocal G4double oldCurveLeng << 960 static G4double oldCurveLength= 0.0;
907 static G4ThreadLocal G4double oldSubStepLe << 961 static G4double oldSubStepLength= 0.0;
908 static G4ThreadLocal G4int oldSubStepNo = << 962 static G4int oldSubStepNo= -1;
909 963
910 G4double subStep_len = 0.0; << 964 G4double subStep_len=0.0;
911 if( curveLen > oldCurveLength ) 965 if( curveLen > oldCurveLength )
912 { 966 {
913 subStep_len= curveLen - oldCurveLength; 967 subStep_len= curveLen - oldCurveLength;
914 } 968 }
915 else if (subStepNo == oldSubStepNo) 969 else if (subStepNo == oldSubStepNo)
916 { 970 {
917 subStep_len= oldSubStepLength; 971 subStep_len= oldSubStepLength;
918 } 972 }
919 oldCurveLength= curveLen; 973 oldCurveLength= curveLen;
920 oldSubStepLength= subStep_len; 974 oldSubStepLength= subStep_len;
921 975
922 G4cout << std::setw(12) << subStep_len << 976 G4cout << std::setw(12) << subStep_len << " ";
923 G4cout << std::setw(12) << subStepSize << 977 G4cout << std::setw(12) << subStepSize << " ";
924 if( requestStep != -1.0 ) 978 if( requestStep != -1.0 )
925 { 979 {
926 G4cout << std::setw( 9) << requestStep < 980 G4cout << std::setw( 9) << requestStep << " ";
927 } 981 }
928 else 982 else
929 { 983 {
930 G4cout << std::setw( 9) << " InitialSte 984 G4cout << std::setw( 9) << " InitialStep " << " ";
931 } 985 }
932 G4cout << G4endl; 986 G4cout << G4endl;
933 } 987 }
934 988
935 // ------------------------------------------- 989 // ---------------------------------------------------------------------------
936 990
937 void G4MagInt_Driver::PrintStatisticsReport() 991 void G4MagInt_Driver::PrintStatisticsReport()
938 { 992 {
939 G4int noPrecBig = 6; << 993 G4int noPrecBig= 6;
940 G4long oldPrec = G4cout.precision(noPrecBig) << 994 G4int oldPrec= G4cout.precision(noPrecBig);
941 995
942 G4cout << "G4MagInt_Driver Statistics of ste 996 G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
943 G4cout << "G4MagInt_Driver: Number of Steps: 997 G4cout << "G4MagInt_Driver: Number of Steps: "
944 << " Total= " << fNoTotalSteps 998 << " Total= " << fNoTotalSteps
945 << " Bad= " << fNoBadSteps 999 << " Bad= " << fNoBadSteps
946 << " Small= " << fNoSmallSteps 1000 << " Small= " << fNoSmallSteps
947 << " Non-initial small= " << (fNoSmal 1001 << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
948 << G4endl; 1002 << G4endl;
>> 1003
>> 1004 #ifdef G4FLD_STATS
>> 1005 G4cout << "MID dyerr: "
>> 1006 << " maximum= " << fDyerr_max
>> 1007 << " Sum small= " << fDyerrPos_smTot
>> 1008 << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
>> 1009 << " vel= " << std::sqrt( fDyerrVel_lgTot )
>> 1010 << " Total h-distance: small= " << fSumH_sm
>> 1011 << " large= " << fSumH_lg
>> 1012 << G4endl;
>> 1013
>> 1014 #if 0
>> 1015 G4int noPrecSmall=4;
>> 1016 // Single line precis of statistics ... optional
>> 1017 G4cout.precision(noPrecSmall);
>> 1018 G4cout << "MIDnums: " << fMinimumStep
>> 1019 << " " << fNoTotalSteps
>> 1020 << " " << fNoSmallSteps
>> 1021 << " " << fNoSmallSteps-fNoInitialSmallSteps
>> 1022 << " " << fNoBadSteps
>> 1023 << " " << fDyerr_max
>> 1024 << " " << fDyerr_mx2
>> 1025 << " " << fDyerrPos_smTot
>> 1026 << " " << fSumH_sm
>> 1027 << " " << fDyerrPos_lgTot
>> 1028 << " " << fDyerrVel_lgTot
>> 1029 << " " << fSumH_lg
>> 1030 << G4endl;
>> 1031 #endif
>> 1032 #endif
>> 1033
949 G4cout.precision(oldPrec); 1034 G4cout.precision(oldPrec);
950 } 1035 }
951 1036
952 // ------------------------------------------- 1037 // ---------------------------------------------------------------------------
953 1038
954 void G4MagInt_Driver::SetSmallestFraction(G4do 1039 void G4MagInt_Driver::SetSmallestFraction(G4double newFraction)
955 { 1040 {
956 if( (newFraction > 1.e-16) && (newFraction < 1041 if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
957 { 1042 {
958 fSmallestFraction= newFraction; 1043 fSmallestFraction= newFraction;
959 } 1044 }
960 else 1045 else
961 { 1046 {
962 std::ostringstream message; << 1047 G4cerr << "Warning: SmallestFraction not changed. " << G4endl
963 message << "Smallest Fraction not changed. << 1048 << " Proposed value was " << newFraction << G4endl
964 << " Proposed value was " << newF << 1049 << " Value must be between 1.e-8 and 1.e-16" << G4endl;
965 << " Value must be between 1.e-8 <<
966 G4Exception("G4MagInt_Driver::SetSmallestF <<
967 "GeomField1001", JustWarning, <<
968 } 1050 }
969 } <<
970 <<
971 void G4MagInt_Driver:: <<
972 GetDerivatives(const G4FieldTrack& y_curr, G4d <<
973 { <<
974 G4double ytemp[G4FieldTrack::ncompSVEC]; <<
975 y_curr.DumpToArray(ytemp); <<
976 pIntStepper->RightHandSide(ytemp, dydx); <<
977 // Avoid virtual call for GetStepper <<
978 // Was: GetStepper()->ComputeRightHandSi <<
979 } <<
980 <<
981 void G4MagInt_Driver::GetDerivatives(const G4F <<
982 G4double <<
983 G4double <<
984 { <<
985 G4double ytemp[G4FieldTrack::ncompSVEC]; <<
986 track.DumpToArray(ytemp); <<
987 pIntStepper->RightHandSide(ytemp, dydx, fi <<
988 } <<
989 <<
990 G4EquationOfMotion* G4MagInt_Driver::GetEquati <<
991 { <<
992 return pIntStepper->GetEquationOfMotion(); <<
993 } <<
994 <<
995 void G4MagInt_Driver::SetEquationOfMotion(G4Eq <<
996 { <<
997 pIntStepper->SetEquationOfMotion(equation) <<
998 } <<
999 <<
1000 const G4MagIntegratorStepper* G4MagInt_Driver <<
1001 { <<
1002 return pIntStepper; <<
1003 } <<
1004 <<
1005 G4MagIntegratorStepper* G4MagInt_Driver::GetS <<
1006 { <<
1007 return pIntStepper; <<
1008 } <<
1009 <<
1010 void G4MagInt_Driver:: <<
1011 RenewStepperAndAdjust(G4MagIntegratorStepper* <<
1012 { <<
1013 pIntStepper = pItsStepper; <<
1014 ReSetParameters(); <<
1015 } <<
1016 <<
1017 void G4MagInt_Driver::StreamInfo( std::ostrea <<
1018 { <<
1019 os << "State of G4MagInt_Driver: " << std <<
1020 os << " Max number of Steps = " << fMaxN <<
1021 << " (base # = " << fMaxStepBase << <<
1022 os << " Safety factor = " << safet <<
1023 os << " Power - shrink = " << pshrn <<
1024 os << " Power - grow = " << pgrow <<
1025 os << " threshold (errcon) = " << errco <<
1026 <<
1027 os << " fMinimumStep = " << fMini <<
1028 os << " Smallest Fraction = " << fSmal <<
1029 <<
1030 os << " No Integrat Vars = " << fNoIn <<
1031 os << " Min No Vars = " << fMinN <<
1032 os << " Num-Vars = " << fNoVa <<
1033 <<
1034 os << " verbose level = " << fVerb <<
1035 os << " Reintegrates = " << DoesR <<
1036 } <<
1037 <<
1038 void PrintInfo( const G4MagInt_Driver & magDr <<
1039 { <<
1040 os << "State of G4MagInt_Driver: " << std <<
1041 os << " Max number of Steps = " << magDr <<
1042 // << " (base # = " << magDrv.fMaxSt <<
1043 os << " Safety factor = " << magDr <<
1044 os << " Power - shrink = " << magDr <<
1045 os << " Power - grow = " << magDr <<
1046 os << " threshold (errcon) = " << magDr <<
1047 <<
1048 os << " fMinimumStep = " << magDr <<
1049 os << " Smallest Fraction = " << magDr <<
1050 <<
1051 /***** <<
1052 os << " No Integrat Vars = " << magDr <<
1053 os << " Min No Vars = " << magDr <<
1054 os << " Num-Vars = " << magDr <<
1055 *****/ <<
1056 os << " verbose level = " << magDr <<
1057 os << " Reintegrates = " << magDr <<
1058 } 1051 }
1059 1052