File indexing completed on 2024-04-06 12:30:26
0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "SimG4Core/MagneticField/interface/FieldStepper.h"
0003 #include "SimG4Core/MagneticField/interface/Field.h"
0004
0005 #include "G4BogackiShampine45.hh"
0006 #include "G4CashKarpRKF45.hh"
0007 #include "G4TCashKarpRKF45.hh"
0008 #include "G4ClassicalRK4.hh"
0009 #include "G4TClassicalRK4.hh"
0010 #include "G4DormandPrince745.hh"
0011 #include "G4TDormandPrince45.hh"
0012 #include "CMSTDormandPrince45.h"
0013 #include "G4HelixExplicitEuler.hh"
0014 #include "G4HelixHeum.hh"
0015 #include "G4HelixImplicitEuler.hh"
0016 #include "G4HelixSimpleRunge.hh"
0017 #include "G4ImplicitEuler.hh"
0018 #include "G4Mag_UsualEqRhs.hh"
0019 #include "G4TMagFieldEquation.hh"
0020 #include "CMSTMagFieldEquation.h"
0021 #include "G4NystromRK4.hh"
0022 #include "G4SimpleHeum.hh"
0023 #include "G4SimpleRunge.hh"
0024 #include "G4TsitourasRK45.hh"
0025
0026 FieldStepper::FieldStepper(G4Mag_UsualEqRhs *eq, double del, const std::string &nam)
0027 : G4MagIntegratorStepper(eq, 6), theEquation(eq), theDelta(del) {
0028 selectStepper(nam);
0029 }
0030
0031 FieldStepper::~FieldStepper() {}
0032
0033 void FieldStepper::Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) {
0034 theStepper->Stepper(y, dydx, h, yout, yerr);
0035 }
0036
0037 G4double FieldStepper::DistChord() const { return theStepper->DistChord(); }
0038
0039 G4int FieldStepper::IntegratorOrder() const { return theStepper->IntegratorOrder(); }
0040
0041 void FieldStepper::selectStepper(const std::string &ss) {
0042 if (ss == "G4ClassicalRK4")
0043 theStepper = new G4ClassicalRK4(theEquation);
0044 else if (ss == "G4TClassicalRK4")
0045 theStepper = new G4TClassicalRK4<G4Mag_UsualEqRhs, 8>(theEquation);
0046 else if (ss == "G4NystromRK4")
0047 theStepper = new G4NystromRK4(theEquation, theDelta);
0048 else if (ss == "G4SimpleRunge")
0049 theStepper = new G4SimpleRunge(theEquation);
0050 else if (ss == "G4SimpleHeum")
0051 theStepper = new G4SimpleHeum(theEquation);
0052 else if (ss == "G4CashKarpRKF45")
0053 theStepper = new G4CashKarpRKF45(theEquation);
0054 else if (ss == "G4TCashKarpRKF45")
0055 theStepper = new G4TCashKarpRKF45<G4Mag_UsualEqRhs>(theEquation);
0056 else if (ss == "G4DormandPrince745")
0057 theStepper = new G4DormandPrince745(theEquation);
0058 else if (ss == "G4TDormandPrince45")
0059 theStepper = new G4TDormandPrince45<G4TMagFieldEquation<sim::Field>>(
0060 dynamic_cast<G4TMagFieldEquation<sim::Field> *>(theEquation));
0061 else if (ss == "CMSTDormandPrince45")
0062 theStepper = new CMSTDormandPrince45<CMSTMagFieldEquation<sim::Field>>(
0063 dynamic_cast<CMSTMagFieldEquation<sim::Field> *>(theEquation));
0064 else if (ss == "G4BogackiShampine45")
0065 theStepper = new G4BogackiShampine45(theEquation);
0066 else if (ss == "G4TsitourasRK45")
0067 theStepper = new G4TsitourasRK45(theEquation);
0068 else if (ss == "G4ImplicitEuler")
0069 theStepper = new G4ImplicitEuler(theEquation);
0070 else if (ss == "G4HelixExplicitEuler")
0071 theStepper = new G4HelixExplicitEuler(theEquation);
0072 else if (ss == "G4HelixImplicitEuler")
0073 theStepper = new G4HelixImplicitEuler(theEquation);
0074 else if (ss == "G4HelixSimpleRunge")
0075 theStepper = new G4HelixSimpleRunge(theEquation);
0076 else if (ss == "G4HelixHeum")
0077 theStepper = new G4HelixHeum(theEquation);
0078 else {
0079 edm::LogWarning("SimG4CoreMagneticField")
0080 << " FieldStepper <" << ss << "> is not known, defaulting to G4ClassicalRK4 ";
0081 theStepper = new G4ClassicalRK4(theEquation);
0082 }
0083 edm::LogVerbatim("SimG4CoreMagneticField") << "### FieldStepper: <" << ss << ">";
0084 }
