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