MagneticField/src/CMSFieldManager.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "SimG4Core/MagneticField/interface/CMSFieldManager.h"
0003 #include "SimG4Core/MagneticField/interface/Field.h"
0004
0005 #include <CLHEP/Units/SystemOfUnits.h>
0006 #include "G4ChordFinder.hh"
0007 #include "G4MagIntegratorStepper.hh"
0008 #include "G4PropagatorInField.hh"
0009 #include "G4Region.hh"
0010 #include "G4RegionStore.hh"
0011 #include "G4Track.hh"
0012
0013 CMSFieldManager::CMSFieldManager()
0014     : G4FieldManager(),
0015       m_currChordFinder(nullptr),
0016       m_chordFinder(nullptr),
0017       m_chordFinderMonopole(nullptr),
0018       m_propagator(nullptr),
0019       m_dChord(0.001),
0020       m_dChordTracker(0.001),
0021       m_dOneStep(0.001),
0022       m_dOneStepTracker(0.0001),
0023       m_dIntersection(0.0001),
0024       m_dInterTracker(1e-6),
0025       m_Rmax2(1.e+6),
0026       m_Zmax(3.e+3),
0027       m_stepMax(1000000.),
0028       m_energyThTracker(1.e+7),
0029       m_energyThreshold(0.0),
0030       m_dChordSimple(0.1),
0031       m_dOneStepSimple(0.1),
0032       m_dIntersectionSimple(0.01),
0033       m_stepMaxSimple(1000.),
0034       m_cfTracker(false),
0035       m_cfVacuum(false) {}
0036
0037 CMSFieldManager::~CMSFieldManager() {
0038   if (m_chordFinder != m_currChordFinder) {
0039     delete m_chordFinder;
0040   }
0041   if (m_chordFinderMonopole != m_currChordFinder) {
0042     delete m_chordFinderMonopole;
0043   }
0044 }
0045
0046 void CMSFieldManager::InitialiseForVolume(const edm::ParameterSet &p,
0047                                           sim::Field *field,
0048                                           G4ChordFinder *cf,
0049                                           G4ChordFinder *cfmon,
0050                                           const std::string &vol,
0051                                           const std::string &type,
0052                                           const std::string &stepper,
0053                                           double delta,
0054                                           G4PropagatorInField *pf) {
0055   double minstep = p.getParameter<double>("MinStep") * CLHEP::mm;
0056   double minEpsStep = p.getUntrackedParameter<double>("MinimumEpsilonStep", 0.00001) * CLHEP::mm;
0057   double maxEpsStep = p.getUntrackedParameter<double>("MaximumEpsilonStep", 0.01) * CLHEP::mm;
0058   int maxLC = (int)p.getUntrackedParameter<double>("MaximumLoopCounts", 1000.);
0059
0060   // double
0061   m_dChord = p.getParameter<double>("DeltaChord") * CLHEP::mm;
0062   m_dChordTracker = p.getParameter<double>("DeltaChord") * CLHEP::mm;
0063   m_dOneStep = p.getParameter<double>("DeltaOneStep") * CLHEP::mm;
0064   m_dOneStepTracker = p.getParameter<double>("DeltaOneStepTracker") * CLHEP::mm;
0065   m_dIntersection = p.getParameter<double>("DeltaIntersection") * CLHEP::mm;
0066   m_dInterTracker = p.getParameter<double>("DeltaIntersectionTracker") * CLHEP::mm;
0067   m_stepMax = p.getParameter<double>("MaxStep") * CLHEP::cm;
0068
0069   m_energyThreshold = p.getParameter<double>("EnergyThSimple") * CLHEP::GeV;
0070   m_energyThTracker = p.getParameter<double>("EnergyThTracker") * CLHEP::GeV;
0071
0072   double rmax = p.getParameter<double>("RmaxTracker") * CLHEP::mm;
0073   m_Rmax2 = rmax * rmax;
0074   m_Zmax = p.getParameter<double>("ZmaxTracker") * CLHEP::mm;
0075
0076   m_dChordSimple = p.getParameter<double>("DeltaChordSimple") * CLHEP::mm;
0077   m_dOneStepSimple = p.getParameter<double>("DeltaOneStepSimple") * CLHEP::mm;
0078   m_dIntersectionSimple = p.getParameter<double>("DeltaIntersectionSimple") * CLHEP::mm;
0079   m_stepMaxSimple = p.getParameter<double>("MaxStepSimple") * CLHEP::cm;
0080
0081   edm::LogVerbatim("SimG4CoreApplication")
0082       << " New CMSFieldManager: LogicalVolume:      <" << vol << ">\n"
0083       << "               Stepper:                   <" << stepper << ">\n"
0084       << "               Field type                 <" << type << ">\n"
0085       << "               Field const delta           " << delta << " mm\n"
0086       << "               MaximumLoopCounts           " << maxLC << "\n"
0087       << "               MinimumEpsilonStep          " << minEpsStep << "\n"
0088       << "               MaximumEpsilonStep          " << maxEpsStep << "\n"
0089       << "               MinStep                     " << minstep << " mm\n"
0090       << "               MaxStep                     " << m_stepMax / CLHEP::cm << " cm\n"
0091       << "               DeltaChord                  " << m_dChord << " mm\n"
0092       << "               DeltaOneStep                " << m_dOneStep << " mm\n"
0093       << "               DeltaIntersection           " << m_dIntersection << " mm\n"
0094       << "               DeltaInterTracker           " << m_dInterTracker << " mm\n"
0095       << "               EnergyThresholdSimple       " << m_energyThreshold / CLHEP::MeV << " MeV\n"
0096       << "               EnergyThresholdTracker      " << m_energyThTracker / CLHEP::MeV << " MeV\n"
0097       << "               DeltaChordSimple            " << m_dChordSimple << " mm\n"
0098       << "               DeltaOneStepSimple          " << m_dOneStepSimple << " mm\n"
0099       << "               DeltaIntersectionSimple     " << m_dIntersectionSimple << " mm\n"
0100       << "               MaxStepInVacuum             " << m_stepMaxSimple / CLHEP::cm << " cm";
0101
0102   // initialisation of chord finders
0103   m_chordFinder = cf;
0104   m_chordFinderMonopole = cfmon;
0105
0106   m_chordFinderMonopole->SetDeltaChord(m_dChord);
0107
0108   // initialisation of field manager
0109   theField.reset(field);
0110   SetDetectorField(field);
0111   SetMinimumEpsilonStep(minEpsStep);
0112   SetMaximumEpsilonStep(maxEpsStep);
0113
0114   // propagater in field
0115   m_propagator = pf;
0116   pf->SetMaxLoopCount(maxLC);
0117   pf->SetMinimumEpsilonStep(minEpsStep);
0118   pf->SetMaximumEpsilonStep(maxEpsStep);
0119
0120   // initial initialisation the default chord finder
0121   setMonopoleTracking(false);
0122
0123   // define regions
0124   std::vector<std::string> rnames = p.getParameter<std::vector<std::string>>("VacRegions");
0125   if (!rnames.empty()) {
0126     std::stringstream ss;
0127     std::vector<G4Region *> *rs = G4RegionStore::GetInstance();
0128     for (auto &regnam : rnames) {
0129       for (auto &reg : *rs) {
0130         if (regnam == reg->GetName()) {
0131           m_regions.push_back(reg);
0132           ss << "  " << regnam;
0133         }
0134       }
0135     }
0136     edm::LogVerbatim("SimG4CoreApplication") << "Simple field integration in G4Regions:\n" << ss.str() << "\n";
0137   }
0138 }
0139
0140 void CMSFieldManager::ConfigureForTrack(const G4Track *track) {
0141   // run time parameters per track
0142   if (track->GetKineticEnergy() > m_energyThTracker && isInsideTracker(track)) {
0143     if (!m_cfTracker) {
0144       setChordFinderForTracker();
0145     }
0146
0147   } else if (track->GetKineticEnergy() <= m_energyThreshold || isInsideVacuum(track)) {
0148     if (!m_cfVacuum) {
0149       setChordFinderForVacuum();
0150     }
0151
0152   } else if (m_cfTracker || m_cfVacuum) {
0153     // restore defaults
0154     setDefaultChordFinder();
0155   }
0156 }
0157
0158 void CMSFieldManager::setMonopoleTracking(G4bool flag) {
0159   if (flag) {
0160     if (m_currChordFinder != m_chordFinderMonopole) {
0161       if (m_cfTracker || m_cfVacuum) {
0162         setDefaultChordFinder();
0163       }
0164       m_currChordFinder = m_chordFinderMonopole;
0165       SetChordFinder(m_currChordFinder);
0166     }
0167   } else {
0168     setDefaultChordFinder();
0169   }
0170   SetFieldChangesEnergy(flag);
0171   m_currChordFinder->ResetStepEstimate();
0172 }
0173
0174 bool CMSFieldManager::isInsideVacuum(const G4Track *track) {
0175   if (!m_regions.empty()) {
0176     const G4Region *reg = track->GetVolume()->GetLogicalVolume()->GetRegion();
0177     for (auto &areg : m_regions) {
0178       if (reg == areg) {
0179         return true;
0180       }
0181     }
0182   }
0183   return false;
0184 }
0185
0186 bool CMSFieldManager::isInsideTracker(const G4Track *track) {
0187   const G4ThreeVector &pos = track->GetPosition();
0188   const double x = pos.x();
0189   const double y = pos.y();
0190   return (x * x + y * y < m_Rmax2 && std::abs(pos.z()) < m_Zmax);
0191 }
0192
0193 void CMSFieldManager::setDefaultChordFinder() {
0194   if (m_currChordFinder != m_chordFinder) {
0195     m_currChordFinder = m_chordFinder;
0196     SetChordFinder(m_currChordFinder);
0197   }
0198   m_currChordFinder->SetDeltaChord(m_dChord);
0199   SetDeltaOneStep(m_dOneStep);
0200   SetDeltaIntersection(m_dIntersection);
0201   m_propagator->SetLargestAcceptableStep(m_stepMax);
0202   m_cfVacuum = m_cfTracker = false;
0203 }
0204
0205 void CMSFieldManager::setChordFinderForTracker() {
0206   if (m_currChordFinder != m_chordFinder) {
0207     m_currChordFinder = m_chordFinder;
0208     SetChordFinder(m_currChordFinder);
0209   }
0210   m_currChordFinder->SetDeltaChord(m_dChordTracker);
0211   SetDeltaOneStep(m_dOneStepTracker);
0212   SetDeltaIntersection(m_dInterTracker);
0213   m_propagator->SetLargestAcceptableStep(m_stepMax);
0214   m_cfVacuum = false;
0215   m_cfTracker = true;
0216 }
0217
0218 void CMSFieldManager::setChordFinderForVacuum() {
0219   if (m_currChordFinder != m_chordFinder) {
0220     m_currChordFinder = m_chordFinder;
0221     SetChordFinder(m_currChordFinder);
0222   }
0223   m_currChordFinder->SetDeltaChord(m_dChordSimple);
0224   SetDeltaOneStep(m_dOneStepSimple);
0225   SetDeltaIntersection(m_dIntersectionSimple);
0226   m_propagator->SetLargestAcceptableStep(m_stepMaxSimple);
0227   m_cfVacuum = true;
0228   m_cfTracker = false;
0229 }