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