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