larg4::LArG4 Class Reference

Runs Geant4 simulation and propagation of electrons and photons to readout. More...

Inheritance diagram for larg4::LArG4:

Public Member Functions
	LArG4 (fhicl::ParameterSet const &pset)
Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDProducer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept
	Producer (fhicl::ParameterSet const &)
	Producer (Producer const &)=delete
	Producer (Producer &&)=delete
Producer &	operator= (Producer const &)=delete
Producer &	operator= (Producer &&)=delete
void	doBeginJob (SharedResources const &resources)
void	doEndJob ()
void	doRespondToOpenInputFile (FileBlock const &fb)
void	doRespondToCloseInputFile (FileBlock const &fb)
void	doRespondToOpenOutputFiles (FileBlock const &fb)
void	doRespondToCloseOutputFiles (FileBlock const &fb)
bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
Public Member Functions inherited from art::Modifier
	~Modifier () noexcept
	Modifier ()
	Modifier (Modifier const &)=delete
	Modifier (Modifier &&)=delete
Modifier &	operator= (Modifier const &)=delete
Modifier &	operator= (Modifier &&)=delete
Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept
	ModuleBase ()
ModuleDescription const &	moduleDescription () const
void	setModuleDescription (ModuleDescription const &)
std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const
void	sortConsumables (std::string const &current_process_name)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Member Functions
void	produce (art::Event &evt) override
void	beginJob () override
void	beginRun (art::Run &run) override
std::unique_ptr< util::PositionInVolumeFilter >	CreateParticleVolumeFilter (std::set< std::string > const &vol_names) const
	Pointer used for correctly updating the clock data state. More...

Private Attributes
std::unique_ptr< g4b::G4Helper >	fG4Help {nullptr}
	G4 interface object. More...
larg4::ParticleListAction *	fparticleListAction
	Geant4 user action to particle information. More...
std::string	fG4PhysListName
	predefined physics list to use if not making a custom one More...
std::string	fG4MacroPath
bool	fCheckOverlaps
	Whether to use the G4 overlap checker. More...
bool	fMakeMCParticles
	Whether to keep a sim::MCParticle list. More...
bool	fdumpParticleList
	Whether each event's sim::ParticleList will be displayed. More...
bool	fdumpSimChannels
	Whether each event's sim::Channel will be displayed. More...
bool	fUseLitePhotons
bool	fStoreReflected {false}
int	fSmartStacking
	Whether to instantiate and use class to. More...
double	fOffPlaneMargin = 0.
std::vector< std::string >	fInputLabels
std::vector< std::string >	fKeepParticlesInVolumes
	Only write particles that have trajectories through these volumes. More...
bool	fSparsifyTrajectories
	Sparsify MCParticle Trajectories. More...
CLHEP::HepRandomEngine &	fEngine
detinfo::DetectorPropertiesData	fDetProp
	Must outlive fAllPhysicsLists! More...
AllPhysicsLists	fAllPhysicsLists
LArVoxelReadoutGeometry *	fVoxelReadoutGeometry

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer
using	WorkerType = WorkerT< EDProducer >
Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >
Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)
Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Runs Geant4 simulation and propagation of electrons and photons to readout.

This module collects generated particles from one or more generators and processes them through Geant4.

Input

The module reads the particles to process from simb::MCTruth records. Each particle generator is required to produce a vector of such records: std::vector<simb::MCTruth>.

The module allows two operation modes:

1. process specific generators: the label of the generator modules to be processed is specified explicitly in LArG4 configuration
2. process all truth information generated so far: no generator is specified in the LArG4 module configuration, and the module will process all data products of type std::vector<simb::MCTruth>, in a non-specified order

For each simb::MCTruth, a Geant4 run is started. The interface with Geant4 is via a helper class provided by nug4. Only the particles in the truth record which have status code (simb::MCParticle::StatusCode()) equal to 1 are processed. These particles are called, in LArG4 jargon, primaries.

Output

The LArG4 module produces:

• a collection of sim::SimChannel: each sim::SimChannel represents the set of energy depositions in liquid argon which drifted and were observed on a certain channel; it includes physics effects like attenuation, diffusion, electric field distortion, etc. Information of the generating Geant4 "track" is retained;
• a collection of sim::SimPhotons or sim::SimPhotonsLite: each sim::SimPhotons represents the set of individual photons reaching a channel of the optical detector; it includes physics effects as well as quantum efficiency of the detector (to reduce data size early in the process); sim::SimPhotonsLite drops the information of the single photons and stores only collective information (e.g. their number).
• a collection of sim::OpDetBacktrackerRecord (to be documented)
• a collection of sim::AuxDetSimChannel (to be documented)
• a collection of simb::MCParticle: the particles generated in the interaction of the primary particles with the material in the world are stored, but minor filtering by geometry and by physics is possible. An association of them with the originating simb::MCTruth object is also produced.

Notes on the conventions

• all and the particles in the truth record (simb::MCTruth) which have status code (simb::MCParticle::StatusCode()) equal to 1 are passed to Geant4. These particles are called, in LArG4 jargon, primaries. The interface with Geant4 is via a helper class provided by nug4.
• normally, information about each particle that Geant4 propagates (which Geant4 calls tracks), primary or not, is saved as an individual simb::MCParticle object into the output particle list. Each simb::MCParticle includes a Geant4-like track ID which is also recorded into each sim::IDE deposited by that particle. This information can be used to track all the deposition from a particle, or to backtrack the particle responsible of a deposition (but see below...). Note that the stored track ID may be different than the one Geant4 used (and, in particular, it's guaranteed to be unique within a sim::LArG4 instance output).
• there are options (some set in sim::LArG4Parameters service) which allow for Geant4 tracks not to be saved as simb::MCParticle (e.g. ParticleKineticEnergyCut, KeepEMShowerDaughters). When these particles have deposited energy, their sim::IDE will report the ID of the first parent Geant4 track which is saved in the simb::MCParticle list, but with its sign flipped. Therefore, when tracking or backtracking (see above), comparisons should be performed using the absolute value of the sim::IDE (e.g. std::abs(ide.trackID)).

Timing

The LArG4 module produces sim::SimChannel objects from generated simb::MCParticle. Each particle ("primary") is assigned the time taken from its vertex (a 4-vector), which is expected to be represented in nanoseconds. The sim::SimChannel object is a collection of sim::IDE in time. The position in the sim::IDE is the location where some ionization occurred. The time associated to a sim::IDE is stored in tick units. The time it represents is the time when the ionization happened, which is the time of the primary particle plus the propagation time to the ionization location, plus the drift time, which the ionized electrons take to reach the anode wire. This time is then shifted to the frame of the electronics time via detinfo::DetectorClocks::G4ToElecTime(), which adds a configurable time offset. The time is converted into ticks via detinfo::DetectorClocks::TPCClock(), and this is the final value associated to the sim::IDE. For a more complete overview, see https://cdcvs.fnal.gov/redmine/projects/larsoft/wiki/Simulation#Simulation-Timing

Randomness

The random number generators used by this process are:

• 'GEANT' instance: used by Geant4
• 'propagation' instance: used in electron propagation

Configuration parameters

• G4PhysListName (string, default: "larg4::PhysicsList"): whether to use the G4 overlap checker, which catches different issues than ROOT
• CheckOverlaps (bool, default: false): whether to use the G4 overlap checker
• DumpParticleList (bool, default: false): whether to print all MCParticles tracked; requires MakeMCParticles being true
• DumpSimChannels (bool, default: false): whether to print all depositions on each SimChannel
• SmartStacking (int, default: 0): whether to use class to dictate how tracks are put on stack (nonzero is on)
• MakeMCParticles (flag, default: true): keep a list of the particles seen in the detector, and eventually save it; you almost always want this on
• KeepParticlesInVolumes (list of strings, default: empty): list of volumes in which to keep simb::MCParticle objects (empty keeps all); requires MakeMCParticles being true
• GeantCommandFile (string, required): G4 macro file to pass to G4Helper for setting G4 command
• Seed (integer, not defined by default): if defined, override the seed for random number generator used in Geant4 simulation (which is obtained from NuRandomService by default)
• PropagationSeed (integer, not defined by default): if defined, override the seed for the random generator used for electrons propagation to the wire planes (obtained from the NuRandomService by default)
• InputLabels (list of strings, default: process all truth): optional list of generator labels whose produced simb::MCTruth will be simulated; if not specified, all simb::MCTruth vector data products are simulated
• ChargeRecoveryMargin (double, default: 0): sets the maximum distance from a plane for the wire charge recovery to occur, in centimeters; for details on how it works, see larg4::LArVoxelReadout::SetOffPlaneChargeRecoveryMargin(). A value of 0 effectively disables this feature. All TPCs will have the same margin applied.

Simulation details

Source of the operational parameters

Some of the physical properties have their values set in FHiCL configuration (e.g. detinfo::LArParameters). Then, GEANT4 is informed of them via larg4::MaterialPropertyLoader. The material property table in GEANT4 is then used by other LArSoft components to discover the parameter values.

Among the parameters registered to GEANT4, the scintillation yields, i.e. how many scintillation photons are produced on average by 1 MeV of deposited energy, are also stored by type of ioniziong particle. These scintillation yields do include a prescale factor (that may include, for example, the photomultiplier quantum efficiency), from the ScintPreScale parameter of detinfo::LArPropertiesStandard or equivalent.

Reflectivity to optical photons

Two models are supported for the simulation of (scintillation) light crossing detector surfaces:

1. the standard one from GEANT4, implemented in G4OpBoundaryProcess
2. a simplified one, implemented in larg4::OpBoundaryProcessSimple

The model is chosen according to the value of detinfo::DetectorProperties::SimpleBoundary(), and the choice is currently exerted by larg4::OpticalPhysics.

The simplified model is faster and simpler: it only deals with absorption and reflection (both specular and diffues). This is the "default" model used in most contexts.

GEANT4 model is more complete and slower. It may take some art to fully configure all the properties of the materials at the sides of the surfaces. The price is a detailed simulation that includes among others refraction and wavelength shifting.

Scintillation

When using the fast optical simulation, which is the "standard" running mode, energy depositions from GEANT4 are "converted" into a number of scintillation photons by the global larg4::IonizationAndScintillation object instance, which internally utilizes the algorithm set up via configuration parameter IonAndScintCalculator in LArG4Parameters service (at the time of writing, "Separate" is supported and "NEST" is accepted too). The number of scintillation photons per energy unit is read from GEANT4 material properties table. It includes already quantum efficiency ("prescale") and it may depend on the type of ionizing particle, depending on the configuration (LArPropertiesStandard parameter ScintByParticleType). This value ("yield") is used as the average of a Poisson distribution from which the actual number of scintillation photons is extracted case by case. The implementation larg4::ISCalculationSeparate may also include medium saturation effects as well, if configured, but only if the scintillation yield is set not to depend on the type of ionizing particle. The number of scintillation photons is then distributed between the fast and slow component by a yield ratio also set in the material parameters, and the single photons are distributed in time accordingly to their component.

Definition at line 314 of file LArG4_module.cc.

Constructor & Destructor Documentation

larg4::LArG4::LArG4 ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 404 of file LArG4_module.cc.

    : art::EDProducer{pset}
    , fG4PhysListName(pset.get<std::string>("G4PhysListName", "larg4::PhysicsList"))
    , fCheckOverlaps(pset.get<bool>("CheckOverlaps", false))
    , fMakeMCParticles(pset.get<bool>("MakeMCParticles", true))
    , fdumpParticleList(pset.get<bool>("DumpParticleList", false))
    , fdumpSimChannels(pset.get<bool>("DumpSimChannels", false))
    , fSmartStacking(pset.get<int>("SmartStacking", 0))
    , fOffPlaneMargin(pset.get<double>("ChargeRecoveryMargin", 0.0))
    , fKeepParticlesInVolumes(pset.get<std::vector<std::string>>("KeepParticlesInVolumes", {}))
    , fSparsifyTrajectories(pset.get<bool>("SparsifyTrajectories", false))
    , fEngine(art::ServiceHandle<rndm::NuRandomService> {}
                ->createEngine(*this, "HepJamesRandom", "propagation", pset, "PropagationSeed"))
    , fDetProp{art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob()}
    , fAllPhysicsLists{fDetProp}
  {
    MF_LOG_DEBUG("LArG4") << "Debug: LArG4()";
    art::ServiceHandle<art::RandomNumberGenerator const> rng;

    if (!fMakeMCParticles) { // configuration option consistency
      if (fdumpParticleList) {
        throw art::Exception(art::errors::Configuration)
          << "Option `DumpParticleList` can't be set if `MakeMCParticles` is unset.\n";
      }
      if (!fKeepParticlesInVolumes.empty()) {
        throw art::Exception(art::errors::Configuration)
          << "Option `KeepParticlesInVolumes` can't be set if `MakeMCParticles` is unset.\n";
      }
    } // if

    if (pset.has_key("Seed")) {
      throw art::Exception(art::errors::Configuration)
        << "The configuration of LArG4 module has the discontinued 'Seed' parameter.\n"
           "Seeds are now controlled by two parameters: 'GEANTSeed' and 'PropagationSeed'.";
    }
    // setup the random number service for Geant4, the "G4Engine" label is a
    // special tag setting up a global engine for use by Geant4/CLHEP;
    // obtain the random seed from NuRandomService,
    // unless overridden in configuration with key "Seed" or "GEANTSeed"
    // FIXME: THIS APPEARS TO BE A NO-OP; IS IT NEEDED?
    (void)art::ServiceHandle<rndm::NuRandomService>()->createEngine(
      *this, "G4Engine", "GEANT", pset, "GEANTSeed");

    //get a list of generators to use, otherwise, we'll end up looking for anything that's
    //made an MCTruth object
    bool useInputLabels =
      pset.get_if_present<std::vector<std::string>>("InputLabels", fInputLabels);
    if (!useInputLabels) fInputLabels.resize(0);

    art::ServiceHandle<sim::LArG4Parameters const> lgp;
    fUseLitePhotons = lgp->UseLitePhotons();

    if (!lgp->NoPhotonPropagation()) {
      try {
        art::ServiceHandle<phot::PhotonVisibilityService const> pvs;
        fStoreReflected = pvs->StoreReflected();
      }
      catch (art::Exception const& e) {
        // If the service is not configured, then just keep the default
        // false for reflected light. If reflected photons are simulated
        // without PVS they will show up in the regular SimPhotons collection
        if (e.categoryCode() != art::errors::ServiceNotFound) throw;
      }

      if (!fUseLitePhotons) {
        produces<std::vector<sim::SimPhotons>>();
        if (fStoreReflected) { produces<std::vector<sim::SimPhotons>>("Reflected"); }
      }
      else {
        produces<std::vector<sim::SimPhotonsLite>>();
        produces<std::vector<sim::OpDetBacktrackerRecord>>();
        if (fStoreReflected) {
          produces<std::vector<sim::SimPhotonsLite>>("Reflected");
          produces<std::vector<sim::OpDetBacktrackerRecord>>("Reflected");
        }
      }
    }

    if (lgp->FillSimEnergyDeposits()) {
      produces<std::vector<sim::SimEnergyDeposit>>("TPCActive");
      produces<std::vector<sim::SimEnergyDeposit>>("Other");
    }

    if (fMakeMCParticles) {
      produces<std::vector<simb::MCParticle>>();
      produces<art::Assns<simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo>>();
    }
    if (!lgp->NoElectronPropagation()) produces<std::vector<sim::SimChannel>>();
    produces<std::vector<sim::AuxDetSimChannel>>();

    // constructor decides if initialized value is a path or an environment variable
    cet::search_path sp("FW_SEARCH_PATH");

    sp.find_file(pset.get<std::string>("GeantCommandFile"), fG4MacroPath);
    struct stat sb;
    if (fG4MacroPath.empty() || stat(fG4MacroPath.c_str(), &sb) != 0)
      // failed to resolve the file name
      throw cet::exception("NoG4Macro") << "G4 macro file " << fG4MacroPath << " not found!\n";
  }

Member Function Documentation

void larg4::LArG4::beginJob ( )

overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 506 of file LArG4_module.cc.

  {
    fG4Help = std::make_unique<g4b::G4Helper>(fG4MacroPath, fG4PhysListName);

    if (fCheckOverlaps) fG4Help->SetOverlapCheck(true);

    art::ServiceHandle<geo::Geometry const> geom;
    fG4Help->ConstructDetector(geom->GDMLFile());

    // Get the logical volume store and assign material properties
    larg4::MaterialPropertyLoader mpl;
    auto const detProp =
      art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
    mpl.GetPropertiesFromServices(detProp);
    mpl.UpdateGeometry(G4LogicalVolumeStore::GetInstance());

    // Tell the detector about the parallel LAr voxel geometry.
    std::vector<G4VUserParallelWorld*> pworlds;

    // Intialize G4 physics and primary generator action
    fG4Help->InitPhysics();

    // create the ionization and scintillation calculator;
    // this is a singleton (!) so it does not make sense
    // to create it in LArVoxelReadoutGeometry
    IonizationAndScintillation::CreateInstance(detProp, fEngine);

    // make a parallel world for each TPC in the detector
    LArVoxelReadoutGeometry::Setup_t readoutGeomSetupData;
    readoutGeomSetupData.readoutSetup.offPlaneMargin = fOffPlaneMargin;
    readoutGeomSetupData.readoutSetup.propGen = &fEngine;

    fVoxelReadoutGeometry =
      new LArVoxelReadoutGeometry("LArVoxelReadoutGeometry", readoutGeomSetupData);
    pworlds.push_back(fVoxelReadoutGeometry);
    pworlds.push_back(
      new OpDetReadoutGeometry(geom->OpDetGeoName(), "OpDetReadoutGeometry", fUseLitePhotons));
    pworlds.push_back(new AuxDetReadoutGeometry("AuxDetReadoutGeometry"));

    fG4Help->SetParallelWorlds(pworlds);

    // moved up
    // Intialize G4 physics and primary generator action
    fG4Help->InitPhysics();

    // Use the UserActionManager to handle all the Geant4 user hooks.
    g4b::UserActionManager* uaManager = g4b::UserActionManager::Instance();

    // User-action class for accumulating LAr voxels.
    art::ServiceHandle<sim::LArG4Parameters const> lgp;

    // User-action class for accumulating particles and trajectories
    // produced in the detector.
    fparticleListAction = new larg4::ParticleListAction(lgp->ParticleKineticEnergyCut(),
                                                        lgp->StoreTrajectories(),
                                                        lgp->KeepEMShowerDaughters(),
                                                        fMakeMCParticles);
    uaManager->AddAndAdoptAction(fparticleListAction);

    // UserActionManager is now configured so continue G4 initialization
    fG4Help->SetUserAction();

    // With an enormous detector with lots of rock ala LAr34 (nee LAr20)
    // we need to be smarter about stacking.
    if (fSmartStacking > 0) {
      G4UserStackingAction* stacking_action = new LArStackingAction(fSmartStacking);
      fG4Help->GetRunManager()->SetUserAction(stacking_action);
    }
  }

void larg4::LArG4::beginRun ( art::Run &  run )

overrideprivatevirtual

Reimplemented from art::EDProducer.

Definition at line 577 of file LArG4_module.cc.

  {
    // prepare the filter object (null if no filtering)
    std::set<std::string> volnameset(fKeepParticlesInVolumes.begin(),
                                     fKeepParticlesInVolumes.end());
    fparticleListAction->ParticleFilter(CreateParticleVolumeFilter(volnameset));
  }

std::unique_ptr< util::PositionInVolumeFilter > larg4::LArG4::CreateParticleVolumeFilter ( std::set< std::string > const &  vol_names ) const

private

Pointer used for correctly updating the clock data state.

Configures and returns a particle filter

Definition at line 586 of file LArG4_module.cc.

  {
    // if we don't have favourite volumes, don't even bother creating a filter
    if (empty(vol_names)) return {};

    auto const& geom = *art::ServiceHandle<geo::Geometry const>();

    std::vector<std::vector<TGeoNode const*>> node_paths = geom.FindAllVolumePaths(vol_names);

    // collection of interesting volumes
    util::PositionInVolumeFilter::AllVolumeInfo_t GeoVolumePairs;
    GeoVolumePairs.reserve(node_paths.size()); // because we are obsessed

    //for each interesting volume, follow the node path and collect
    //total rotations and translations
    for (size_t iVolume = 0; iVolume < node_paths.size(); ++iVolume) {
      std::vector<TGeoNode const*> path = node_paths[iVolume];

      auto pTransl = new TGeoTranslation(0., 0., 0.);
      auto pRot = new TGeoRotation();
      for (TGeoNode const* node : path) {
        TGeoTranslation thistranslate(*node->GetMatrix());
        TGeoRotation thisrotate(*node->GetMatrix());
        pTransl->Add(&thistranslate);
        *pRot = *pRot * thisrotate;
      }

      // for some reason, pRot and pTransl don't have tr and rot bits set
      // correctly make new translations and rotations so bits are set correctly
      auto pTransl2 = new TGeoTranslation(
        pTransl->GetTranslation()[0], pTransl->GetTranslation()[1], pTransl->GetTranslation()[2]);
      double phi = 0., theta = 0., psi = 0.;
      pRot->GetAngles(phi, theta, psi);
      auto pRot2 = new TGeoRotation();
      pRot2->SetAngles(phi, theta, psi);

      auto pTransf = new TGeoCombiTrans(*pTransl2, *pRot2);
      GeoVolumePairs.emplace_back(node_paths[iVolume].back()->GetVolume(), pTransf);
    }

    return std::make_unique<util::PositionInVolumeFilter>(std::move(GeoVolumePairs));
  } // CreateParticleVolumeFilter()

void larg4::LArG4::produce ( art::Event &  evt )

overrideprivatevirtual

The main routine of this module: Fetch the primary particles from the event, simulate their evolution in the detector, and produce the detector response.

Implements art::EDProducer.

Definition at line 630 of file LArG4_module.cc.

  {
    MF_LOG_DEBUG("LArG4") << "produce()";
    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
    auto const detProp =
      art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
    LArVoxelReadoutGeometry::Sentry const set_for_event{fVoxelReadoutGeometry, clockData, detProp};

    // loop over the lists and put the particles and voxels into the event as
    // collections
    auto scCol = std::make_unique<std::vector<sim::SimChannel>>();
    auto adCol = std::make_unique<std::vector<sim::AuxDetSimChannel>>();
    auto tpassn =
      fMakeMCParticles ?
      std::make_unique<art::Assns<simb::MCTruth, simb::MCParticle, sim::GeneratedParticleInfo>>() :
      nullptr;
    auto partCol =
      fMakeMCParticles ?
      std::make_unique<std::vector<simb::MCParticle>>() :
      nullptr;
    auto PhotonCol = std::make_unique<std::vector<sim::SimPhotons>>();
    auto PhotonColRefl = std::make_unique<std::vector<sim::SimPhotons>>();
    auto LitePhotonCol = std::make_unique<std::vector<sim::SimPhotonsLite>>();
    auto LitePhotonColRefl = std::make_unique<std::vector<sim::SimPhotonsLite>>();
    auto cOpDetBacktrackerRecordCol = std::make_unique<std::vector<sim::OpDetBacktrackerRecord>>();
    auto cOpDetBacktrackerRecordColRefl =
      std::make_unique<std::vector<sim::OpDetBacktrackerRecord>>();

    std::optional<art::PtrMaker<simb::MCParticle>> makeMCPartPtr;
    if (fMakeMCParticles) makeMCPartPtr.emplace(evt);

    // for energy deposits
    auto edepCol_TPCActive = std::make_unique<std::vector<sim::SimEnergyDeposit>>();
    auto edepCol_Other = std::make_unique<std::vector<sim::SimEnergyDeposit>>();

    // Fetch the lists of LAr voxels and particles.
    art::ServiceHandle<sim::LArG4Parameters const> lgp;
    art::ServiceHandle<geo::Geometry const> geom;

    // Clear the detected photon table
    OpDetPhotonTable::Instance()->ClearTable(geom->NOpDets());
    if (lgp->FillSimEnergyDeposits()) OpDetPhotonTable::Instance()->ClearEnergyDeposits();

    // reset the track ID offset as we have a new collection of interactions
    fparticleListAction->ResetTrackIDOffset();

    //look to see if there is any MCTruth information for this
    //event
    std::vector<art::Handle<std::vector<simb::MCTruth>>> mclists;
    if (empty(fInputLabels))
      //evt.getManyByType(mclists);
      mclists = evt.getMany<std::vector<simb::MCTruth>>();
    else {
      mclists.resize(fInputLabels.size());
      for (size_t i = 0; i < fInputLabels.size(); i++)
        evt.getByLabel(fInputLabels[i], mclists[i]);
    }

    unsigned int nGeneratedParticles = 0;

    // Need to process Geant4 simulation for each interaction separately.
    for (size_t mcl = 0; mcl < mclists.size(); ++mcl) {

      art::Handle<std::vector<simb::MCTruth>> mclistHandle = mclists[mcl];

      for (size_t m = 0; m < mclistHandle->size(); ++m) {
        art::Ptr<simb::MCTruth> mct(mclistHandle, m);

        MF_LOG_DEBUG("LArG4") << *(mct.get());

        // The following tells Geant4 to track the particles in this interaction.
        fG4Help->G4Run(mct);

        if (!partCol) continue;
        assert(tpassn);

        // receive the particle list
        sim::ParticleList particleList = fparticleListAction->YieldList();

        for (auto const& partPair : particleList) {
          simb::MCParticle& p = *(partPair.second);
          ++nGeneratedParticles;

          // if the particle has been marked as dropped, we don't save it
          // (as of LArSoft ~v5.6 this does not ever happen because
          // ParticleListAction has already taken care of deleting them)
          if (ParticleListAction::isDropped(&p)) continue;

          sim::GeneratedParticleInfo const truthInfo{
            fparticleListAction->GetPrimaryTruthIndex(p.TrackId())};
          if (!truthInfo.hasGeneratedParticleIndex() && (p.Mother() == 0)) {
            // this means it's primary but with no information; logic error!!
            art::Exception error(art::errors::LogicError);
            error << "Failed to match primary particle:\n";
            sim::dump::DumpMCParticle(error, p, "  ");
            error << "\nwith particles from the truth record '"
                  << mclistHandle.provenance()->inputTag() << "':\n";
            sim::dump::DumpMCTruth(error, *mct, 2U, "  "); // 2 points per line
            error << "\n";
            throw error;
          }

          if (fSparsifyTrajectories) p.SparsifyTrajectory();

          partCol->push_back(std::move(p));

          tpassn->addSingle(mct, (*makeMCPartPtr)(partCol->size() - 1), truthInfo);

        } // for(particleList)

        // Has the user request a detailed dump of the output objects?
        if (fdumpParticleList) {
          mf::LogInfo("LArG4") << "Dump sim::ParticleList; size()=" << particleList.size() << "\n"
                               << particleList;
        }
      }

    } // end loop over interactions

    // get the electrons from the LArVoxelReadout sensitive detector
    // Get the sensitive-detector manager.
    G4SDManager* sdManager = G4SDManager::GetSDMpointer();

    // Find the sensitive detector with the name "LArVoxelSD".
    auto theOpDetDet = dynamic_cast<OpDetSensitiveDetector*>(
      sdManager->FindSensitiveDetector("OpDetSensitiveDetector"));

    // Store the contents of the detected photon table
    //
    if (theOpDetDet) {

      if (!lgp->NoPhotonPropagation()) {

        for (int Reflected = 0; Reflected <= 1; Reflected++) {
          if (Reflected && !fStoreReflected) continue;

          if (!fUseLitePhotons) {
            MF_LOG_DEBUG("Optical") << "Storing OpDet Hit Collection in Event";
            std::vector<sim::SimPhotons>& ThePhotons =
              OpDetPhotonTable::Instance()->GetPhotons(Reflected);
            if (Reflected)
              PhotonColRefl->reserve(ThePhotons.size());
            else
              PhotonCol->reserve(ThePhotons.size());
            for (auto& it : ThePhotons) {
              if (Reflected)
                PhotonColRefl->push_back(std::move(it));
              else
                PhotonCol->push_back(std::move(it));
            }
          }
          else {
            MF_LOG_DEBUG("Optical") << "Storing OpDet Hit Collection in Event";

            std::map<int, std::map<int, int>> ThePhotons =
              OpDetPhotonTable::Instance()->GetLitePhotons(Reflected);

            if (size(ThePhotons) > 0) {
              LitePhotonCol->reserve(ThePhotons.size());
              for (auto const& [opChannel, detectedPhotons] : ThePhotons) {
                sim::SimPhotonsLite ph;
                ph.OpChannel = opChannel;
                ph.DetectedPhotons = detectedPhotons;
                if (Reflected)
                  LitePhotonColRefl->push_back(std::move(ph));
                else
                  LitePhotonCol->push_back(std::move(ph));
              }
            }
          }
          if (Reflected)
            *cOpDetBacktrackerRecordColRefl =
              OpDetPhotonTable::Instance()->YieldReflectedOpDetBacktrackerRecords();
          else
            *cOpDetBacktrackerRecordCol =
              OpDetPhotonTable::Instance()->YieldOpDetBacktrackerRecords();
        }
      } //end if no photon propagation

      if (lgp->FillSimEnergyDeposits()) {
        // we steal the only existing copy of the energy deposit map. Oink!
        auto edepMap = OpDetPhotonTable::Instance()->YieldSimEnergyDeposits();
        for (auto& [volumeName, edepCol] : edepMap) {
          // note: constant reference to a (smart) pointer to non-const data
          auto const& destColl =
            boost::contains(volumeName, "TPCActive") ? edepCol_TPCActive : edepCol_Other;
          append(*destColl, std::move(edepCol));
        } // for
      }
    } //end if theOpDetDet

    if (!lgp->NoElectronPropagation()) {

      // only put the sim::SimChannels into the event once, not once for every
      // MCTruth in the event

      std::set<LArVoxelReadout*> ReadoutList; // to be cleared later on

      for (unsigned int c = 0; c < geom->Ncryostats(); ++c) {

        // map to keep track of which channels we already have SimChannels for in scCol
        // remake this map on each cryostat as channels ought not to be shared between
        // cryostats, just between TPC's

        std::map<unsigned int, unsigned int> channelToscCol;

        unsigned int ntpcs = geom->Cryostat(c).NTPC();
        for (unsigned int t = 0; t < ntpcs; ++t) {
          std::string name("LArVoxelSD");
          std::ostringstream sstr;
          sstr << name << "_Cryostat" << c << "_TPC" << t;

          // try first to find the sensitive detector specific for this TPC;
          // do not bother writing on screen if there is none (yet)
          G4VSensitiveDetector* sd = sdManager->FindSensitiveDetector(sstr.str(), false);
          // if there is none, catch the general one (called just "LArVoxelSD")
          if (!sd) sd = sdManager->FindSensitiveDetector(name, false);
          // If this didn't work, then a sensitive detector with
          // the name "LArVoxelSD" does not exist.
          if (!sd) {
            throw cet::exception("LArG4")
              << "Sensitive detector for cryostat " << c << " TPC " << t << " not found (neither '"
              << sstr.str() << "' nor '" << name << "' exist)\n";
          }

          // Convert the G4VSensitiveDetector* to a LArVoxelReadout*.
          auto larVoxelReadout = dynamic_cast<LArVoxelReadout*>(sd);

          // If this didn't work, there is a "LArVoxelSD" detector, but
          // it's not a LArVoxelReadout object.
          if (!larVoxelReadout) {
            throw cet::exception("LArG4")
              << "Sensitive detector '" << sd->GetName() << "' is not a LArVoxelReadout object\n";
          }

          LArVoxelReadout::ChannelMap_t& channels = larVoxelReadout->GetSimChannelMap(c, t);
          if (!empty(channels)) {
            MF_LOG_DEBUG("LArG4") << "now put " << channels.size() << " SimChannels from C=" << c
                                  << " T=" << t << " into the event";
          }

          for (auto ch_pair : channels) {
            sim::SimChannel& sc = ch_pair.second;

            // push sc onto scCol but only if we haven't already put something in scCol for this channel.
            // if we have, then merge the ionization deposits.  Skip the check if we only have one TPC

            if (ntpcs > 1) {
              unsigned int ichan = sc.Channel();
              auto itertest = channelToscCol.find(ichan);
              if (itertest == channelToscCol.end()) {
                channelToscCol[ichan] = scCol->size();
                scCol->emplace_back(std::move(sc));
              }
              else {
                unsigned int idtest = itertest->second;
                auto const& tdcideMap = sc.TDCIDEMap();
                for (auto const& tdcide : tdcideMap) {
                  for (auto const& ide : tdcide.second) {
                    double xyz[3] = {ide.x, ide.y, ide.z};
                    scCol->at(idtest).AddIonizationElectrons(
                      ide.trackID, tdcide.first, ide.numElectrons, xyz, ide.energy);
                  } // end loop to add ionization electrons to  scCol->at(idtest)
                }   // end loop over tdc to vector<sim::IDE> map
              }     // end if check to see if we've put SimChannels in for ichan yet or not
            }
            else {
              scCol->emplace_back(std::move(sc));
            } // end of check if we only have one TPC (skips check for multiple simchannels if we have just one TPC)
          }   // end loop over simchannels for this TPC

          // mark it for clearing
          ReadoutList.insert(const_cast<LArVoxelReadout*>(larVoxelReadout));

        } // end loop over tpcs
      }   // end loop over cryostats

      for (LArVoxelReadout* larVoxelReadout : ReadoutList) {
        larVoxelReadout->ClearSimChannels();
      }
    } //endif electron prop

    // only put the sim::AuxDetSimChannels into the event once, not once for every
    // MCTruth in the event

    adCol->reserve(geom->NAuxDets());
    for (unsigned int a = 0; a < geom->NAuxDets(); ++a) {

      // there should always be at least one senstive volume because
      // we make one for the full aux det if none are specified in the
      // gdml file - see AuxDetGeo.cxx
      for (size_t sv = 0; sv < geom->AuxDet(a).NSensitiveVolume(); ++sv) {

        // N.B. this name convention is used when creating the
        //      AuxDetReadout SD in AuxDetReadoutGeometry
        std::stringstream name;
        name << "AuxDetSD_AuxDet" << a << "_" << sv;
        G4VSensitiveDetector* sd = sdManager->FindSensitiveDetector(name.str().c_str());
        if (!sd) {
          throw cet::exception("LArG4")
            << "Sensitive detector '" << name.str() << "' does not exist\n";
        }

        // Convert the G4VSensitiveDetector* to a AuxDetReadout*.
        larg4::AuxDetReadout* auxDetReadout = dynamic_cast<larg4::AuxDetReadout*>(sd);

        MF_LOG_DEBUG("LArG4") << "now put the AuxDetSimTracks in the event";

        const sim::AuxDetSimChannel adsc = auxDetReadout->GetAuxDetSimChannel();
        adCol->push_back(adsc);
        auxDetReadout->clear();
      }
    } // Loop over AuxDets

    if (partCol) {
      mf::LogInfo("LArG4") << "Geant4 simulated " << nGeneratedParticles << " MC particles, we keep "
                           << partCol->size() << " .";
    }

    if (fdumpSimChannels) {
      mf::LogVerbatim("DumpSimChannels")
        << "Event " << evt.id() << ": " << scCol->size() << " channels with signal";
      unsigned int nChannels = 0;
      for (const sim::SimChannel& sc : *scCol) {
        mf::LogVerbatim out("DumpSimChannels");
        out << " #" << nChannels << ": ";
        // dump indenting with "    ", but not on the first line
        sc.Dump(out, "  ");
        ++nChannels;
      } // for
    }   // if dump SimChannels

    if (!lgp->NoElectronPropagation()) evt.put(std::move(scCol));

    evt.put(std::move(adCol));
    if (partCol) evt.put(std::move(partCol));
    if (tpassn) evt.put(std::move(tpassn));
    if (!lgp->NoPhotonPropagation()) {
      if (!fUseLitePhotons) {
        evt.put(std::move(PhotonCol));
        if (fStoreReflected) evt.put(std::move(PhotonColRefl), "Reflected");
      }
      else {
        evt.put(std::move(LitePhotonCol));
        evt.put(std::move(cOpDetBacktrackerRecordCol));
        if (fStoreReflected) {
          evt.put(std::move(LitePhotonColRefl), "Reflected");
          evt.put(std::move(cOpDetBacktrackerRecordColRefl), "Reflected");
        }
      }
    }

    if (lgp->FillSimEnergyDeposits()) {
      evt.put(std::move(edepCol_TPCActive), "TPCActive");
      evt.put(std::move(edepCol_Other), "Other");
    }
    return;
  } // LArG4::produce()

Member Data Documentation

AllPhysicsLists larg4::LArG4::fAllPhysicsLists

private

Definition at line 352 of file LArG4_module.cc.

bool larg4::LArG4::fCheckOverlaps

private

Whether to use the G4 overlap checker.

Definition at line 333 of file LArG4_module.cc.

detinfo::DetectorPropertiesData larg4::LArG4::fDetProp

private

Must outlive fAllPhysicsLists!

Definition at line 351 of file LArG4_module.cc.

bool larg4::LArG4::fdumpParticleList

private

Whether each event's sim::ParticleList will be displayed.

Definition at line 335 of file LArG4_module.cc.

bool larg4::LArG4::fdumpSimChannels

private

Whether each event's sim::Channel will be displayed.

Definition at line 336 of file LArG4_module.cc.

CLHEP::HepRandomEngine& larg4::LArG4::fEngine

private

Random-number engine for IonizationAndScintillation initialization

Definition at line 348 of file LArG4_module.cc.

std::unique_ptr<g4b::G4Helper> larg4::LArG4::fG4Help {nullptr}

private

G4 interface object.

Definition at line 326 of file LArG4_module.cc.

std::string larg4::LArG4::fG4MacroPath

private

directory path for Geant4 macro file to be executed before main MC processing.

Definition at line 331 of file LArG4_module.cc.

std::string larg4::LArG4::fG4PhysListName

private

predefined physics list to use if not making a custom one

Definition at line 330 of file LArG4_module.cc.

std::vector<std::string> larg4::LArG4::fInputLabels

private

Definition at line 342 of file LArG4_module.cc.

std::vector<std::string> larg4::LArG4::fKeepParticlesInVolumes

private

Only write particles that have trajectories through these volumes.

Definition at line 344 of file LArG4_module.cc.

bool larg4::LArG4::fMakeMCParticles

private

Whether to keep a sim::MCParticle list.

Definition at line 334 of file LArG4_module.cc.

double larg4::LArG4::fOffPlaneMargin = 0.

private

Off-plane charge recovery margin dictate how tracks are put on stack.

Definition at line 340 of file LArG4_module.cc.

larg4::ParticleListAction* larg4::LArG4::fparticleListAction

private

Initial value:

{
      nullptr}

Geant4 user action to particle information.

Definition at line 327 of file LArG4_module.cc.

int larg4::LArG4::fSmartStacking

private

Whether to instantiate and use class to.

Definition at line 339 of file LArG4_module.cc.

bool larg4::LArG4::fSparsifyTrajectories

private

Sparsify MCParticle Trajectories.

Definition at line 346 of file LArG4_module.cc.

bool larg4::LArG4::fStoreReflected {false}

private

Definition at line 338 of file LArG4_module.cc.

bool larg4::LArG4::fUseLitePhotons

private

Definition at line 337 of file LArG4_module.cc.

LArVoxelReadoutGeometry* larg4::LArG4::fVoxelReadoutGeometry

private

Initial value:

{
      nullptr}

Definition at line 353 of file LArG4_module.cc.

---

The documentation for this class was generated from the following file:

• larsim/larsim/LegacyLArG4/LArG4_module.cc