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