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