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