MagneticField/src/FieldStepper.cc

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 }