FieldBuilder.cc

Go to the documentation of this file.

#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "SimG4Core/MagneticField/interface/CMSFieldManager.h"
#include "SimG4Core/MagneticField/interface/Field.h"
#include "SimG4Core/MagneticField/interface/FieldBuilder.h"
#include "SimG4Core/MagneticField/interface/FieldStepper.h"
#include "SimG4Core/MagneticField/interface/MonopoleEquation.h"

#include <CLHEP/Units/SystemOfUnits.h>
#include "G4ChordFinder.hh"
#include "G4ClassicalRK4.hh"
#include "G4FieldManager.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4TMagFieldEquation.hh"
#include "CMSTMagFieldEquation.h"
#include "G4PropagatorInField.hh"

using namespace sim;

FieldBuilder::FieldBuilder(const MagneticField *f, const edm::ParameterSet &p) : theTopVolume(nullptr), thePSet(p) {
  theDelta = p.getParameter<double>("delta") * CLHEP::mm;
  theField = new Field(f, theDelta);
  theFieldEquation = nullptr;
}

FieldBuilder::~FieldBuilder() {}

void FieldBuilder::build(CMSFieldManager *fM, G4PropagatorInField *fP) {
  edm::ParameterSet thePSetForGMFM = thePSet.getParameter<edm::ParameterSet>("ConfGlobalMFM");
  std::string volName = thePSetForGMFM.getParameter<std::string>("Volume");
  edm::ParameterSet volPSet = thePSetForGMFM.getParameter<edm::ParameterSet>(volName);

  configureForVolume(volName, volPSet, fM, fP);

  edm::LogVerbatim("SimG4CoreMagneticField") << " FieldBuilder::build: Global magnetic field is used";
}

void FieldBuilder::configureForVolume(const std::string &volName,
                                      edm::ParameterSet &volPSet,
                                      CMSFieldManager *fM,
                                      G4PropagatorInField *fP) {
  G4LogicalVolumeStore *theStore = G4LogicalVolumeStore::GetInstance();
  for (auto vol : *theStore) {
    if ((std::string)vol->GetName() == volName) {
      theTopVolume = vol;
      break;
    }
  }

  std::string fieldType = volPSet.getParameter<std::string>("Type");
  std::string stepper = volPSet.getParameter<std::string>("Stepper");

  edm::ParameterSet stpPSet = volPSet.getParameter<edm::ParameterSet>("StepperParam");
  double minStep = stpPSet.getParameter<double>("MinStep") * CLHEP::mm;

  if (stepper == "CMSTDormandPrince45") {
    theFieldEquation = new CMSTMagFieldEquation<sim::Field>(theField);
  } else if (stepper == "G4TDormandPrince45") {
    theFieldEquation = new G4TMagFieldEquation<sim::Field>(theField);
  } else {
    theFieldEquation = new G4Mag_UsualEqRhs(theField);
  }

  FieldStepper *dStepper = new FieldStepper(theFieldEquation, theDelta, stepper);
  G4ChordFinder *cf = new G4ChordFinder(theField, minStep, dStepper);

  MonopoleEquation *monopoleEquation = new MonopoleEquation(theField);
  G4MagIntegratorStepper *mStepper = new G4ClassicalRK4(monopoleEquation, 8);
  G4ChordFinder *cfmon = new G4ChordFinder(theField, minStep, mStepper);

  fM->InitialiseForVolume(stpPSet, theField, cf, cfmon, volName, fieldType, stepper, theDelta, fP);
}