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 // class G4HelixMixedStepper 26 // class G4HelixMixedStepper 27 // 27 // 28 // Class description: 28 // Class description: 29 // 29 // 30 // G4HelixMixedStepper split the Method used f 30 // G4HelixMixedStepper split the Method used for Integration in two: 31 // 31 // 32 // If Stepping Angle ( h / R_curve) < pi/3 32 // If Stepping Angle ( h / R_curve) < pi/3 33 // use Stepper for small step(Classical 33 // use Stepper for small step(ClassicalRK4 by default) 34 // Else use HelixExplicitEuler Stepper 34 // Else use HelixExplicitEuler Stepper 35 // 35 // 36 // Created: T.Nikitina, CERN - 18.05.2007, der 36 // Created: T.Nikitina, CERN - 18.05.2007, derived from G4ExactHelicalStepper 37 // ------------------------------------------- 37 // ------------------------------------------------------------------------- 38 38 39 #include "G4HelixMixedStepper.hh" 39 #include "G4HelixMixedStepper.hh" 40 #include "G4PhysicalConstants.hh" 40 #include "G4PhysicalConstants.hh" 41 #include "G4ClassicalRK4.hh" 41 #include "G4ClassicalRK4.hh" 42 #include "G4CashKarpRKF45.hh" 42 #include "G4CashKarpRKF45.hh" 43 #include "G4SimpleRunge.hh" 43 #include "G4SimpleRunge.hh" 44 #include "G4HelixImplicitEuler.hh" 44 #include "G4HelixImplicitEuler.hh" 45 #include "G4HelixExplicitEuler.hh" 45 #include "G4HelixExplicitEuler.hh" 46 #include "G4HelixSimpleRunge.hh" 46 #include "G4HelixSimpleRunge.hh" 47 #include "G4ExactHelixStepper.hh" 47 #include "G4ExactHelixStepper.hh" 48 #include "G4ExplicitEuler.hh" 48 #include "G4ExplicitEuler.hh" 49 #include "G4ImplicitEuler.hh" 49 #include "G4ImplicitEuler.hh" 50 #include "G4SimpleHeum.hh" 50 #include "G4SimpleHeum.hh" 51 #include "G4RKG3_Stepper.hh" 51 #include "G4RKG3_Stepper.hh" 52 #include "G4NystromRK4.hh" 52 #include "G4NystromRK4.hh" 53 53 54 // Additional potential steppers 54 // Additional potential steppers 55 #include "G4DormandPrince745.hh" 55 #include "G4DormandPrince745.hh" 56 #include "G4BogackiShampine23.hh" 56 #include "G4BogackiShampine23.hh" 57 #include "G4BogackiShampine45.hh" 57 #include "G4BogackiShampine45.hh" 58 #include "G4TsitourasRK45.hh" 58 #include "G4TsitourasRK45.hh" 59 59 60 #include "G4ThreeVector.hh" 60 #include "G4ThreeVector.hh" 61 #include "G4LineSection.hh" 61 #include "G4LineSection.hh" 62 62 63 // ------------------------------------------- 63 // --------------------------------------------------------------------------- 64 G4HelixMixedStepper:: 64 G4HelixMixedStepper:: 65 G4HelixMixedStepper(G4Mag_EqRhs* EqRhs, 65 G4HelixMixedStepper(G4Mag_EqRhs* EqRhs, 66 G4int stepperNumber 66 G4int stepperNumber, 67 G4double angleThreshol 67 G4double angleThreshold) 68 : G4MagHelicalStepper(EqRhs) 68 : G4MagHelicalStepper(EqRhs) 69 { 69 { 70 if( angleThreshold < 0.0 ) 70 if( angleThreshold < 0.0 ) 71 { 71 { 72 fAngle_threshold = (1.0/3.0)*pi; 72 fAngle_threshold = (1.0/3.0)*pi; 73 } 73 } 74 else 74 else 75 { 75 { 76 fAngle_threshold = angleThreshold; 76 fAngle_threshold = angleThreshold; 77 } 77 } 78 78 79 if(stepperNumber<0) 79 if(stepperNumber<0) 80 { 80 { 81 // stepperNumber = 4; // Default is RK4 81 // stepperNumber = 4; // Default is RK4 (original) 82 stepperNumber = 745; // Default is Dorm 82 stepperNumber = 745; // Default is DormandPrince745 (ie DoPri5) 83 // stepperNumber = 8; // Default is Cash 83 // stepperNumber = 8; // Default is CashKarp 84 } 84 } 85 85 86 fStepperNumber = stepperNumber; // Store th 86 fStepperNumber = stepperNumber; // Store the choice 87 fRK4Stepper = SetupStepper(EqRhs, fStepper 87 fRK4Stepper = SetupStepper(EqRhs, fStepperNumber); 88 } 88 } 89 89 90 // ------------------------------------------- 90 // --------------------------------------------------------------------------- 91 G4HelixMixedStepper::~G4HelixMixedStepper() 91 G4HelixMixedStepper::~G4HelixMixedStepper() 92 { 92 { 93 delete fRK4Stepper; 93 delete fRK4Stepper; 94 if (fVerbose>0) { PrintCalls(); } 94 if (fVerbose>0) { PrintCalls(); } 95 } 95 } 96 96 97 // ------------------------------------------- 97 // --------------------------------------------------------------------------- 98 void G4HelixMixedStepper::Stepper( const G4do 98 void G4HelixMixedStepper::Stepper( const G4double yInput[], // [7] 99 const G4do 99 const G4double dydx[], // [7] 100 G4do 100 G4double Step, 101 G4do 101 G4double yOut[], // [7] 102 G4do 102 G4double yErr[]) 103 { 103 { 104 // Estimation of the Stepping Angle 104 // Estimation of the Stepping Angle 105 // 105 // 106 G4ThreeVector Bfld; 106 G4ThreeVector Bfld; 107 MagFieldEvaluate(yInput, Bfld); 107 MagFieldEvaluate(yInput, Bfld); 108 108 109 G4double Bmag = Bfld.mag(); 109 G4double Bmag = Bfld.mag(); 110 const G4double* pIn = yInput+3; 110 const G4double* pIn = yInput+3; 111 G4ThreeVector initVelocity = G4ThreeVector( 111 G4ThreeVector initVelocity = G4ThreeVector( pIn[0], pIn[1], pIn[2] ); 112 G4double velocityVal = initVelocity.mag(); 112 G4double velocityVal = initVelocity.mag(); 113 113 114 const G4double R_1 = std::abs(GetInverseCurv 114 const G4double R_1 = std::abs(GetInverseCurve(velocityVal,Bmag)); 115 // curv = inverse Radius 115 // curv = inverse Radius 116 G4double Ang_curve = R_1 * Step; 116 G4double Ang_curve = R_1 * Step; 117 // SetAngCurve(Ang_curve); 117 // SetAngCurve(Ang_curve); 118 // SetCurve(std::abs(1/R_1)); 118 // SetCurve(std::abs(1/R_1)); 119 119 120 if(Ang_curve < fAngle_threshold) 120 if(Ang_curve < fAngle_threshold) 121 { 121 { 122 ++fNumCallsRK4; 122 ++fNumCallsRK4; 123 fRK4Stepper->Stepper(yInput,dydx,Step,yOut 123 fRK4Stepper->Stepper(yInput,dydx,Step,yOut,yErr); 124 } 124 } 125 else 125 else 126 { 126 { 127 constexpr G4int nvar = 6 ; 127 constexpr G4int nvar = 6 ; 128 constexpr G4int nvarMax = 8 ; 128 constexpr G4int nvarMax = 8 ; 129 G4double yTemp[nvarMax], yIn[nvarMax] 129 G4double yTemp[nvarMax], yIn[nvarMax], yTemp2[nvarMax]; 130 G4ThreeVector Bfld_midpoint; 130 G4ThreeVector Bfld_midpoint; 131 131 132 SetAngCurve(Ang_curve); 132 SetAngCurve(Ang_curve); 133 SetCurve(std::abs(1.0/R_1)); 133 SetCurve(std::abs(1.0/R_1)); 134 ++fNumCallsHelix; 134 ++fNumCallsHelix; 135 135 136 // Saving yInput because yInput and yOut c 136 // Saving yInput because yInput and yOut can be aliases for same array 137 // 137 // 138 for(G4int i=0; i<nvar; ++i) 138 for(G4int i=0; i<nvar; ++i) 139 { 139 { 140 yIn[i]=yInput[i]; 140 yIn[i]=yInput[i]; 141 } 141 } 142 142 143 G4double halfS = Step * 0.5; 143 G4double halfS = Step * 0.5; 144 144 145 // 1. Do first half step and full step 145 // 1. Do first half step and full step 146 // 146 // 147 AdvanceHelix(yIn, Bfld, halfS, yTemp, yTem 147 AdvanceHelix(yIn, Bfld, halfS, yTemp, yTemp2); // yTemp2 for s=2*h (halfS) 148 148 149 MagFieldEvaluate(yTemp, Bfld_midpoint) ; 149 MagFieldEvaluate(yTemp, Bfld_midpoint) ; 150 150 151 // 2. Do second half step - with revised f 151 // 2. Do second half step - with revised field 152 // NOTE: Could avoid this call if 'Bfld_m 152 // NOTE: Could avoid this call if 'Bfld_midpoint == Bfld' 153 // or diff 'almost' zero 153 // or diff 'almost' zero 154 // 154 // 155 AdvanceHelix(yTemp, Bfld_midpoint, halfS, 155 AdvanceHelix(yTemp, Bfld_midpoint, halfS, yOut); 156 // Not requesting y at s=2*h (halfS) 156 // Not requesting y at s=2*h (halfS) 157 157 158 // 3. Estimate the integration error 158 // 3. Estimate the integration error 159 // should be (nearly) zero if Bfield= c 159 // should be (nearly) zero if Bfield= constant 160 // 160 // 161 for(G4int i=0; i<nvar; ++i) 161 for(G4int i=0; i<nvar; ++i) 162 { 162 { 163 yErr[i] = yOut[i] - yTemp2[i]; 163 yErr[i] = yOut[i] - yTemp2[i]; 164 } 164 } 165 } 165 } 166 } 166 } 167 167 168 // ------------------------------------------- 168 // --------------------------------------------------------------------------- 169 void G4HelixMixedStepper::DumbStepper( const G 169 void G4HelixMixedStepper::DumbStepper( const G4double yIn[], 170 G 170 G4ThreeVector Bfld, 171 G 171 G4double h, 172 G 172 G4double yOut[] ) 173 { 173 { 174 AdvanceHelix(yIn, Bfld, h, yOut); 174 AdvanceHelix(yIn, Bfld, h, yOut); 175 } 175 } 176 176 177 // ------------------------------------------- 177 // --------------------------------------------------------------------------- 178 G4double G4HelixMixedStepper::DistChord() cons 178 G4double G4HelixMixedStepper::DistChord() const 179 { 179 { 180 // Implementation : must check whether h/R > 180 // Implementation : must check whether h/R > 2 pi !! 181 // If( h/R < pi) use G4LineSection::DistL 181 // If( h/R < pi) use G4LineSection::DistLine 182 // Else DistChord=R_helix 182 // Else DistChord=R_helix 183 // 183 // 184 G4double distChord; 184 G4double distChord; 185 G4double Ang_curve=GetAngCurve(); 185 G4double Ang_curve=GetAngCurve(); 186 186 187 if(Ang_curve<=pi) 187 if(Ang_curve<=pi) 188 { 188 { 189 distChord=GetRadHelix()*(1-std::cos(0.5*An 189 distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve)); 190 } 190 } 191 else 191 else 192 { 192 { 193 if(Ang_curve<twopi) 193 if(Ang_curve<twopi) 194 { 194 { 195 distChord=GetRadHelix()*(1+std::cos(0.5* 195 distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve))); 196 } 196 } 197 else 197 else 198 { 198 { 199 distChord=2.*GetRadHelix(); 199 distChord=2.*GetRadHelix(); 200 } 200 } 201 } 201 } 202 202 203 return distChord; 203 return distChord; 204 } 204 } 205 205 206 // ------------------------------------------- 206 // --------------------------------------------------------------------------- 207 void G4HelixMixedStepper::PrintCalls() 207 void G4HelixMixedStepper::PrintCalls() 208 { 208 { 209 G4cout << "In HelixMixedStepper::Number of c 209 G4cout << "In HelixMixedStepper::Number of calls to smallStepStepper = " 210 << fNumCallsRK4 210 << fNumCallsRK4 211 << " and Number of calls to Helix = 211 << " and Number of calls to Helix = " << fNumCallsHelix << G4endl; 212 } 212 } 213 213 214 // ------------------------------------------- 214 // --------------------------------------------------------------------------- 215 G4MagIntegratorStepper* 215 G4MagIntegratorStepper* 216 G4HelixMixedStepper::SetupStepper(G4Mag_EqRhs* 216 G4HelixMixedStepper::SetupStepper(G4Mag_EqRhs* pE, G4int StepperNumber) 217 { 217 { 218 G4MagIntegratorStepper* pStepper; 218 G4MagIntegratorStepper* pStepper; 219 if (fVerbose>0) { G4cout << " G4HelixMixedSt 219 if (fVerbose>0) { G4cout << " G4HelixMixedStepper: "; 220 } 220 } 221 switch ( StepperNumber ) 221 switch ( StepperNumber ) 222 { 222 { 223 // Robust, classic method 223 // Robust, classic method 224 case 4: 224 case 4: 225 pStepper = new G4ClassicalRK4( pE ); 225 pStepper = new G4ClassicalRK4( pE ); 226 if (fVerbose>0) { G4cout << "G4Classic 226 if (fVerbose>0) { G4cout << "G4ClassicalRK4"; } 227 break; 227 break; 228 228 229 // Steppers with embedded estimation of 229 // Steppers with embedded estimation of error 230 case 8: 230 case 8: 231 pStepper = new G4CashKarpRKF45( pE ); 231 pStepper = new G4CashKarpRKF45( pE ); 232 if (fVerbose>0) { G4cout << "G4CashKar 232 if (fVerbose>0) { G4cout << "G4CashKarpRKF45"; } 233 break; 233 break; 234 case 13: 234 case 13: 235 pStepper = new G4NystromRK4( pE ); 235 pStepper = new G4NystromRK4( pE ); 236 if (fVerbose>0) { G4cout << "G4Nystrom 236 if (fVerbose>0) { G4cout << "G4NystromRK4"; } 237 break; 237 break; 238 238 239 // Lowest order RK Stepper - experimenta 239 // Lowest order RK Stepper - experimental 240 case 1: 240 case 1: 241 pStepper = new G4ImplicitEuler( pE ); 241 pStepper = new G4ImplicitEuler( pE ); 242 if (fVerbose>0) { G4cout << "G4Implici 242 if (fVerbose>0) { G4cout << "G4ImplicitEuler"; } 243 break; 243 break; 244 244 245 // Lower order RK Steppers - ok overall, 245 // Lower order RK Steppers - ok overall, good for uneven fields 246 case 2: 246 case 2: 247 pStepper = new G4SimpleRunge( pE ); 247 pStepper = new G4SimpleRunge( pE ); 248 if (fVerbose>0) { G4cout << "G4SimpleR 248 if (fVerbose>0) { G4cout << "G4SimpleRunge"; } 249 break; 249 break; 250 case 3: 250 case 3: 251 pStepper = new G4SimpleHeum( pE ); 251 pStepper = new G4SimpleHeum( pE ); 252 if (fVerbose>0) { G4cout << "G4SimpleH 252 if (fVerbose>0) { G4cout << "G4SimpleHeum"; } 253 break; 253 break; 254 case 23: 254 case 23: 255 pStepper = new G4BogackiShampine23( pE 255 pStepper = new G4BogackiShampine23( pE ); 256 if (fVerbose>0) { G4cout << "G4Bogacki 256 if (fVerbose>0) { G4cout << "G4BogackiShampine23"; } 257 break; 257 break; 258 258 259 // Higher order RK Steppers 259 // Higher order RK Steppers 260 // for smoother fields and high accuracy 260 // for smoother fields and high accuracy requirements 261 case 45: 261 case 45: 262 pStepper = new G4BogackiShampine45( pE 262 pStepper = new G4BogackiShampine45( pE ); 263 if (fVerbose>0) { G4cout << "G4Bogacki 263 if (fVerbose>0) { G4cout << "G4BogackiShampine45"; } 264 break; 264 break; 265 case 145: 265 case 145: 266 pStepper = new G4TsitourasRK45( pE ); 266 pStepper = new G4TsitourasRK45( pE ); 267 if (fVerbose>0) { G4cout << "G4Tsitour 267 if (fVerbose>0) { G4cout << "G4TsitourasRK45"; } 268 break; 268 break; 269 case 745: 269 case 745: 270 pStepper = new G4DormandPrince745( pE 270 pStepper = new G4DormandPrince745( pE ); 271 if (fVerbose>0) { G4cout << "G4Dormand 271 if (fVerbose>0) { G4cout << "G4DormandPrince745"; } 272 break; 272 break; 273 273 274 // Helical Steppers 274 // Helical Steppers 275 case 6: 275 case 6: 276 pStepper = new G4HelixImplicitEuler( p 276 pStepper = new G4HelixImplicitEuler( pE ); 277 if (fVerbose>0) { G4cout << "G4HelixIm 277 if (fVerbose>0) { G4cout << "G4HelixImplicitEuler"; } 278 break; 278 break; 279 case 7: 279 case 7: 280 pStepper = new G4HelixSimpleRunge( pE 280 pStepper = new G4HelixSimpleRunge( pE ); 281 if (fVerbose>0) { G4cout << "G4HelixSi 281 if (fVerbose>0) { G4cout << "G4HelixSimpleRunge"; } 282 break; 282 break; 283 case 5: 283 case 5: 284 pStepper = new G4HelixExplicitEuler( p 284 pStepper = new G4HelixExplicitEuler( pE ); 285 if (fVerbose>0) { G4cout << "G4HelixEx 285 if (fVerbose>0) { G4cout << "G4HelixExplicitEuler"; } 286 break; // Since Helix Explicit is use 286 break; // Since Helix Explicit is used for long steps, 287 // this is useful only to measu 287 // this is useful only to measure overhead. 288 // Exact Helix - likely good only for ca 288 // Exact Helix - likely good only for cases of 289 // i) uniform field (potentia 289 // i) uniform field (potentially over small distances) 290 // ii) segmented uniform field 290 // ii) segmented uniform field (maybe) 291 case 9: 291 case 9: 292 pStepper = new G4ExactHelixStepper( pE 292 pStepper = new G4ExactHelixStepper( pE ); 293 if (fVerbose>0) { G4cout << "G4ExactHe 293 if (fVerbose>0) { G4cout << "G4ExactHelixStepper"; } 294 break; 294 break; 295 case 10: 295 case 10: 296 pStepper = new G4RKG3_Stepper( pE ); 296 pStepper = new G4RKG3_Stepper( pE ); 297 if (fVerbose>0) { G4cout << "G4RKG3_St 297 if (fVerbose>0) { G4cout << "G4RKG3_Stepper"; } 298 break; 298 break; 299 299 300 // Low Order Steppers - not good except 300 // Low Order Steppers - not good except for very weak fields 301 case 11: 301 case 11: 302 pStepper = new G4ExplicitEuler( pE ); 302 pStepper = new G4ExplicitEuler( pE ); 303 if (fVerbose>0) { G4cout << "G4Explici 303 if (fVerbose>0) { G4cout << "G4ExplicitEuler"; } 304 break; 304 break; 305 case 12: 305 case 12: 306 pStepper = new G4ImplicitEuler( pE ); 306 pStepper = new G4ImplicitEuler( pE ); 307 if (fVerbose>0) { G4cout << "G4Implici 307 if (fVerbose>0) { G4cout << "G4ImplicitEuler"; } 308 break; 308 break; 309 309 310 case 0: 310 case 0: 311 case -1: 311 case -1: 312 default: 312 default: 313 pStepper = new G4DormandPrince745( pE 313 pStepper = new G4DormandPrince745( pE ); // Was G4ClassicalRK4( pE ); 314 if (fVerbose>0) { G4cout << "G4Dormand 314 if (fVerbose>0) { G4cout << "G4DormandPrince745 (Default)"; } 315 break; 315 break; 316 } 316 } 317 317 318 if(fVerbose>0) 318 if(fVerbose>0) 319 { 319 { 320 G4cout << " chosen as stepper for small st 320 G4cout << " chosen as stepper for small steps in G4HelixMixedStepper." 321 << G4endl; 321 << G4endl; 322 } 322 } 323 323 324 return pStepper; 324 return pStepper; 325 } 325 } 326 326