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