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 // History: >> 37 // Derived from ExactHelicalStepper 18/05/07 >> 38 // 37 // ------------------------------------------- 39 // ------------------------------------------------------------------------- 38 40 39 #include "G4HelixMixedStepper.hh" 41 #include "G4HelixMixedStepper.hh" 40 #include "G4PhysicalConstants.hh" 42 #include "G4PhysicalConstants.hh" 41 #include "G4ClassicalRK4.hh" 43 #include "G4ClassicalRK4.hh" 42 #include "G4CashKarpRKF45.hh" 44 #include "G4CashKarpRKF45.hh" 43 #include "G4SimpleRunge.hh" 45 #include "G4SimpleRunge.hh" 44 #include "G4HelixImplicitEuler.hh" 46 #include "G4HelixImplicitEuler.hh" 45 #include "G4HelixExplicitEuler.hh" 47 #include "G4HelixExplicitEuler.hh" 46 #include "G4HelixSimpleRunge.hh" 48 #include "G4HelixSimpleRunge.hh" 47 #include "G4ExactHelixStepper.hh" 49 #include "G4ExactHelixStepper.hh" 48 #include "G4ExplicitEuler.hh" 50 #include "G4ExplicitEuler.hh" 49 #include "G4ImplicitEuler.hh" 51 #include "G4ImplicitEuler.hh" 50 #include "G4SimpleHeum.hh" 52 #include "G4SimpleHeum.hh" 51 #include "G4RKG3_Stepper.hh" 53 #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 54 60 #include "G4ThreeVector.hh" 55 #include "G4ThreeVector.hh" 61 #include "G4LineSection.hh" 56 #include "G4LineSection.hh" 62 << 57 G4HelixMixedStepper::G4HelixMixedStepper(G4Mag_EqRhs *EqRhs,G4int fStepperNumber) 63 // ------------------------------------------- << 64 G4HelixMixedStepper:: << 65 G4HelixMixedStepper(G4Mag_EqRhs* EqRhs, << 66 G4int stepperNumber << 67 G4double angleThreshol << 68 : G4MagHelicalStepper(EqRhs) 58 : G4MagHelicalStepper(EqRhs) >> 59 69 { 60 { 70 if( angleThreshold < 0.0 ) << 61 SetVerbose(1); fNumCallsRK4=0; fNumCallsHelix=0; 71 { << 62 if(!fStepperNumber) fStepperNumber=4; 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 63 fRK4Stepper = SetupStepper(EqRhs, fStepperNumber); 88 } 64 } 89 65 90 // ------------------------------------------- << 91 G4HelixMixedStepper::~G4HelixMixedStepper() << 92 { << 93 delete fRK4Stepper; << 94 if (fVerbose>0) { PrintCalls(); } << 95 } << 96 66 97 // ------------------------------------------- << 67 G4HelixMixedStepper::~G4HelixMixedStepper() { 98 void G4HelixMixedStepper::Stepper( const G4do << 68 99 const G4do << 69 delete(fRK4Stepper); 100 G4do << 70 if (fVerbose>0){ PrintCalls();}; 101 G4do << 71 } 102 G4do << 72 void G4HelixMixedStepper::Stepper( const G4double yInput[7], >> 73 const G4double dydx[7], >> 74 G4double Step, >> 75 G4double yOut[7], >> 76 G4double yErr[]) >> 77 103 { 78 { 104 // Estimation of the Stepping Angle << 79 105 // << 80 //Estimation of the Stepping Angle >> 81 106 G4ThreeVector Bfld; 82 G4ThreeVector Bfld; 107 MagFieldEvaluate(yInput, Bfld); << 83 MagFieldEvaluate(yInput, Bfld); 108 << 84 109 G4double Bmag = Bfld.mag(); 85 G4double Bmag = Bfld.mag(); 110 const G4double* pIn = yInput+3; << 86 const G4double *pIn = yInput+3; 111 G4ThreeVector initVelocity = G4ThreeVector( << 87 G4ThreeVector initVelocity= G4ThreeVector( pIn[0], pIn[1], pIn[2]); 112 G4double velocityVal = initVelocity.mag(); << 88 G4double velocityVal = initVelocity.mag(); 113 << 89 G4double R_1; 114 const G4double R_1 = std::abs(GetInverseCurv << 90 G4double Ang_curve; 115 // curv = inverse Radius << 91 116 G4double Ang_curve = R_1 * Step; << 92 R_1=std::abs(GetInverseCurve(velocityVal,Bmag)); 117 // SetAngCurve(Ang_curve); << 93 Ang_curve=R_1*Step; 118 // SetCurve(std::abs(1/R_1)); << 94 SetAngCurve(Ang_curve); 119 << 95 SetCurve(std::abs(1/R_1)); 120 if(Ang_curve < fAngle_threshold) << 96 121 { << 97 122 ++fNumCallsRK4; << 98 if(Ang_curve<0.33*pi){ 123 fRK4Stepper->Stepper(yInput,dydx,Step,yOut << 99 fNumCallsRK4++; 124 } << 100 fRK4Stepper->Stepper(yInput,dydx,Step,yOut,yErr); 125 else << 126 { << 127 constexpr G4int nvar = 6 ; << 128 constexpr G4int nvarMax = 8 ; << 129 G4double yTemp[nvarMax], yIn[nvarMax] << 130 G4ThreeVector Bfld_midpoint; << 131 101 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 102 145 // 1. Do first half step and full step << 103 } 146 // << 104 else{ 147 AdvanceHelix(yIn, Bfld, halfS, yTemp, yTem << 105 fNumCallsHelix++; 148 << 106 const G4int nvar = 6 ; 149 MagFieldEvaluate(yTemp, Bfld_midpoint) ; << 107 G4int i; 150 << 108 G4double yTemp[7], yIn[7] ; 151 // 2. Do second half step - with revised f << 109 G4double yTemp2[7]; 152 // NOTE: Could avoid this call if 'Bfld_m << 110 G4ThreeVector Bfld_midpoint; 153 // or diff 'almost' zero << 111 // Saving yInput because yInput and yOut can be aliases for same array 154 // << 112 for(i=0;i<nvar;i++) yIn[i]=yInput[i]; 155 AdvanceHelix(yTemp, Bfld_midpoint, halfS, << 113 156 // Not requesting y at s=2*h (halfS) << 114 G4double h = Step * 0.5; 157 << 115 // Do two half steps and full step 158 // 3. Estimate the integration error << 116 AdvanceHelix(yIn, Bfld, h, yTemp,yTemp2); 159 // should be (nearly) zero if Bfield= c << 117 MagFieldEvaluate(yTemp, Bfld_midpoint) ; 160 // << 118 AdvanceHelix(yTemp, Bfld_midpoint, h, yOut); 161 for(G4int i=0; i<nvar; ++i) << 119 // Error estimation 162 { << 120 for(i=0;i<nvar;i++) { 163 yErr[i] = yOut[i] - yTemp2[i]; << 121 yErr[i] = yOut[i] - yTemp2[i] ; >> 122 >> 123 } 164 } 124 } 165 } << 125 >> 126 >> 127 >> 128 166 } 129 } 167 130 168 // ------------------------------------------- << 131 void 169 void G4HelixMixedStepper::DumbStepper( const G << 132 G4HelixMixedStepper::DumbStepper( const G4double yIn[], 170 G << 133 G4ThreeVector Bfld, 171 G << 134 G4double h, 172 G << 135 G4double yOut[]) 173 { 136 { 174 AdvanceHelix(yIn, Bfld, h, yOut); << 137 175 } << 138 >> 139 AdvanceHelix(yIn, Bfld, h, yOut); 176 140 177 // ------------------------------------------- << 141 178 G4double G4HelixMixedStepper::DistChord() cons << 142 >> 143 } >> 144 >> 145 G4double G4HelixMixedStepper::DistChord() const 179 { 146 { 180 // Implementation : must check whether h/R > 147 // Implementation : must check whether h/R > 2 pi !! 181 // If( h/R < pi) use G4LineSection::DistL 148 // If( h/R < pi) use G4LineSection::DistLine 182 // Else DistChord=R_helix 149 // Else DistChord=R_helix 183 // 150 // 184 G4double distChord; 151 G4double distChord; 185 G4double Ang_curve=GetAngCurve(); 152 G4double Ang_curve=GetAngCurve(); 186 << 153 187 if(Ang_curve<=pi) << 154 188 { << 155 if(Ang_curve<=pi){ 189 distChord=GetRadHelix()*(1-std::cos(0.5*An << 156 distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve)); 190 } << 157 } 191 else << 158 else 192 { << 159 if(Ang_curve<twopi){ 193 if(Ang_curve<twopi) << 160 distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve))); 194 { << 161 } 195 distChord=GetRadHelix()*(1+std::cos(0.5* << 162 else{ 196 } << 163 distChord=2.*GetRadHelix(); 197 else << 164 } 198 { << 165 199 distChord=2.*GetRadHelix(); << 166 200 } << 167 201 } << 202 << 203 return distChord; 168 return distChord; >> 169 204 } 170 } 205 << 206 // ------------------------------------------- 171 // --------------------------------------------------------------------------- 207 void G4HelixMixedStepper::PrintCalls() 172 void G4HelixMixedStepper::PrintCalls() 208 { 173 { 209 G4cout << "In HelixMixedStepper::Number of c << 174 G4cout<<"In HelixMixedStepper::Number of calls to smallStepStepper = "<<fNumCallsRK4 210 << fNumCallsRK4 << 175 <<" and Number of calls to Helix = "<<fNumCallsHelix<<G4endl; 211 << " and Number of calls to Helix = << 212 } 176 } 213 177 214 // ------------------------------------------- << 178 215 G4MagIntegratorStepper* << 179 216 G4HelixMixedStepper::SetupStepper(G4Mag_EqRhs* << 180 G4MagIntegratorStepper* G4HelixMixedStepper:: SetupStepper(G4Mag_EqRhs* pE, G4int StepperNumber) 217 { 181 { 218 G4MagIntegratorStepper* pStepper; 182 G4MagIntegratorStepper* pStepper; 219 if (fVerbose>0) { G4cout << " G4HelixMixedSt << 183 if (fVerbose>0)G4cout<<"In G4HelixMixedStepper Stepper for small steps is "; 220 } << 221 switch ( StepperNumber ) 184 switch ( StepperNumber ) 222 { << 185 { 223 // Robust, classic method << 186 case 0: pStepper = new G4ExplicitEuler( pE ); if (fVerbose>0)G4cout<<"G4ExplicitEuler"<<G4endl; break; 224 case 4: << 187 case 1: pStepper = new G4ImplicitEuler( pE ); if (fVerbose>0)G4cout<<"G4ImplicitEuler"<<G4endl; break; 225 pStepper = new G4ClassicalRK4( pE ); << 188 case 2: pStepper = new G4SimpleRunge( pE ); if (fVerbose>0)G4cout<<"G4SimpleRunge"<<G4endl; break; 226 if (fVerbose>0) { G4cout << "G4Classic << 189 case 3: pStepper = new G4SimpleHeum( pE ); if (fVerbose>0)G4cout<<"G4SimpleHeum"<<G4endl;break; 227 break; << 190 case 4: pStepper = new G4ClassicalRK4( pE ); if (fVerbose>0)G4cout<<"G4ClassicalRK4"<<G4endl; break; 228 << 191 case 5: pStepper = new G4HelixExplicitEuler( pE ); if (fVerbose>0)G4cout<<"G4HelixExplicitEuler"<<G4endl; break; 229 // Steppers with embedded estimation of << 192 case 6: pStepper = new G4HelixImplicitEuler( pE ); if (fVerbose>0)G4cout<<"G4HelixImplicitEuler"<<G4endl; break; 230 case 8: << 193 case 7: pStepper = new G4HelixSimpleRunge( pE ); if (fVerbose>0)G4cout<<"G4HelixSimpleRunge"<<G4endl; break; 231 pStepper = new G4CashKarpRKF45( pE ); << 194 case 8: pStepper = new G4CashKarpRKF45( pE ); if (fVerbose>0)G4cout<<"G4CashKarpRKF45"<<G4endl; break; 232 if (fVerbose>0) { G4cout << "G4CashKar << 195 case 9: pStepper = new G4ExactHelixStepper( pE ); if (fVerbose>0)G4cout<<"G4ExactHelixStepper"<<G4endl; break; 233 break; << 196 case 10: pStepper = new G4RKG3_Stepper( pE ); if (fVerbose>0)G4cout<<"G4RKG3_Stepper"<<G4endl; break; 234 case 13: << 197 235 pStepper = new G4NystromRK4( pE ); << 198 default: pStepper = new G4ClassicalRK4( pE );G4cout<<"Default G4ClassicalRK4"<<G4endl; break; 236 if (fVerbose>0) { G4cout << "G4Nystrom << 199 237 break; << 200 } 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 } << 317 << 318 if(fVerbose>0) << 319 { << 320 G4cout << " chosen as stepper for small st << 321 << G4endl; << 322 } << 323 << 324 return pStepper; 201 return pStepper; 325 } 202 } 326 203