Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // class G4HelixMixedStepper 27 // 28 // Class description: 29 // 30 // G4HelixMixedStepper split the Method used f 31 // 32 // If Stepping Angle ( h / R_curve) < pi/3 33 // use Stepper for small step(Classical 34 // Else use HelixExplicitEuler Stepper 35 // 36 // Created: T.Nikitina, CERN - 18.05.2007, der 37 // ------------------------------------------- 38 39 #include "G4HelixMixedStepper.hh" 40 #include "G4PhysicalConstants.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" 61 #include "G4LineSection.hh" 62 63 // ------------------------------------------- 64 G4HelixMixedStepper:: 65 G4HelixMixedStepper(G4Mag_EqRhs* EqRhs, 66 G4int stepperNumber 67 G4double angleThreshol 68 : G4MagHelicalStepper(EqRhs) 69 { 70 if( angleThreshold < 0.0 ) 71 { 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 } 89 90 // ------------------------------------------- 91 G4HelixMixedStepper::~G4HelixMixedStepper() 92 { 93 delete fRK4Stepper; 94 if (fVerbose>0) { PrintCalls(); } 95 } 96 97 // ------------------------------------------- 98 void G4HelixMixedStepper::Stepper( const G4do 99 const G4do 100 G4do 101 G4do 102 G4do 103 { 104 // Estimation of the Stepping Angle 105 // 106 G4ThreeVector Bfld; 107 MagFieldEvaluate(yInput, Bfld); 108 109 G4double Bmag = Bfld.mag(); 110 const G4double* pIn = yInput+3; 111 G4ThreeVector initVelocity = G4ThreeVector( 112 G4double velocityVal = initVelocity.mag(); 113 114 const G4double R_1 = std::abs(GetInverseCurv 115 // curv = inverse Radius 116 G4double Ang_curve = R_1 * Step; 117 // SetAngCurve(Ang_curve); 118 // SetCurve(std::abs(1/R_1)); 119 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 145 // 1. Do first half step and full step 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 158 // 3. Estimate the integration error 159 // should be (nearly) zero if Bfield= c 160 // 161 for(G4int i=0; i<nvar; ++i) 162 { 163 yErr[i] = yOut[i] - yTemp2[i]; 164 } 165 } 166 } 167 168 // ------------------------------------------- 169 void G4HelixMixedStepper::DumbStepper( const G 170 G 171 G 172 G 173 { 174 AdvanceHelix(yIn, Bfld, h, yOut); 175 } 176 177 // ------------------------------------------- 178 G4double G4HelixMixedStepper::DistChord() cons 179 { 180 // Implementation : must check whether h/R > 181 // If( h/R < pi) use G4LineSection::DistL 182 // Else DistChord=R_helix 183 // 184 G4double distChord; 185 G4double Ang_curve=GetAngCurve(); 186 187 if(Ang_curve<=pi) 188 { 189 distChord=GetRadHelix()*(1-std::cos(0.5*An 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 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 } 317 318 if(fVerbose>0) 319 { 320 G4cout << " chosen as stepper for small st 321 << G4endl; 322 } 323 324 return pStepper; 325 } 326