LArG4_module.cc

Go to the documentation of this file.

#include "nug4/G4Base/G4Helper.h"

// C++ Includes
#include <cassert>
#include <map>
#include <set>
#include <sstream>
#include <sys/stat.h>

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Persistency/Common/PtrMaker.h"
#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Utilities/Exception.h"
#include "cetlib/search_path.h"
#include "cetlib_except/exception.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// art extensions
#include "nurandom/RandomUtils/NuRandomService.h"

// LArSoft Includes
#include "larcore/Geometry/AuxDetGeometry.h"
#include "larcore/Geometry/Geometry.h"
#include "larcore/Geometry/WireReadout.h"
#include "larcorealg/CoreUtils/ParticleFilters.h" // util::PositionInVolumeFilter
#include "lardata/DetectorInfoServices/DetectorClocksService.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardataalg/MCDumpers/MCDumpers.h" // sim::dump namespace
#include "lardataobj/MCBase/MCParticleLite.h"
#include "lardataobj/Simulation/AuxDetSimChannel.h"
#include "lardataobj/Simulation/GeneratedParticleInfo.h"
#include "lardataobj/Simulation/OpDetBacktrackerRecord.h"
#include "lardataobj/Simulation/SimChannel.h"
#include "lardataobj/Simulation/SimPhotons.h"
#include "larsim/LegacyLArG4/AllPhysicsLists.h"
#include "larsim/LegacyLArG4/AuxDetReadout.h"
#include "larsim/LegacyLArG4/AuxDetReadoutGeometry.h"
#include "larsim/LegacyLArG4/IonizationAndScintillation.h"
#include "larsim/LegacyLArG4/LArStackingAction.h"
#include "larsim/LegacyLArG4/LArVoxelReadout.h"
#include "larsim/LegacyLArG4/LArVoxelReadoutGeometry.h"
#include "larsim/LegacyLArG4/MaterialPropertyLoader.h"
#include "larsim/LegacyLArG4/OpDetPhotonTable.h"
#include "larsim/LegacyLArG4/OpDetReadoutGeometry.h"
#include "larsim/LegacyLArG4/OpDetSensitiveDetector.h"
#include "larsim/LegacyLArG4/ParticleListAction.h"
#include "larsim/PhotonPropagation/PhotonVisibilityService.h"
#include "larsim/Simulation/LArG4Parameters.h"
#include "nug4/G4Base/UserActionManager.h"
#include "nug4/ParticleNavigation/ParticleList.h"
#include "nusimdata/SimulationBase/MCTruth.h"

// G4 Includes
#include "Geant4/G4LogicalVolumeStore.hh"
#include "Geant4/G4RunManager.hh"
#include "Geant4/G4SDManager.hh"
#include "Geant4/G4VSensitiveDetector.hh"
#include "Geant4/G4VUserDetectorConstruction.hh"

//For energy depositions
#include "lardataobj/Simulation/SimEnergyDeposit.h"

// Boost includes
#include "boost/algorithm/string.hpp"

namespace larg4 {

  // Forward declarations within namespace.
  class LArVoxelListAction;

  class LArG4 : public art::EDProducer {
  public:
    explicit LArG4(fhicl::ParameterSet const& pset);

  private:
    void produce(art::Event& evt) override;
    void beginJob() override;
    void beginRun(art::Run& run) override;

    std::unique_ptr<g4b::G4Helper> fG4Help{nullptr}; 
    larg4::ParticleListAction* fparticleListAction{
      nullptr}; 

    std::string fG4PhysListName; 
    std::string fG4MacroPath;    
    bool fCheckOverlaps;         
    bool fMakeMCParticles;       
    bool
      fStoreDroppedMCParticles; 
    bool fdumpParticleList;     
    bool fdumpSimChannels;      
    bool fUseLitePhotons;
    bool fStoreReflected{false};
    int fSmartStacking;          
    double fOffPlaneMargin = 0.; 
    std::vector<std::string> fInputLabels;
    std::vector<std::string>
      fKeepParticlesInVolumes; 

    bool fSparsifyTrajectories; 

    CLHEP::HepRandomEngine& fEngine; 

    detinfo::DetectorPropertiesData fDetProp; 
    AllPhysicsLists fAllPhysicsLists;
    LArVoxelReadoutGeometry* fVoxelReadoutGeometry{
      nullptr}; 

    std::unique_ptr<util::PositionInVolumeFilter> CreateParticleVolumeFilter(
      std::set<std::string> const& vol_names) const;
  };

} // namespace LArG4

namespace {

  // ---------------------------------------------------------------------------
  template <typename T>
  std::vector<T>& append(std::vector<T>& dest, std::vector<T>&& source)
  {
    if (empty(dest))
      dest = std::move(source);
    else {
      dest.insert(dest.end(), std::move_iterator{begin(source)}, std::move_iterator{end(source)});
      source = std::vector<T>{}; // ensure the old memory is released
    }
    return dest;
  }
  // ---------------------------------------------------------------------------

} // local namespace

namespace larg4 {

  //----------------------------------------------------------------------
  // Constructor
  LArG4::LArG4(fhicl::ParameterSet const& pset)
    : art::EDProducer{pset}
    , fG4PhysListName(pset.get<std::string>("G4PhysListName", "larg4::PhysicsList"))
    , fCheckOverlaps(pset.get<bool>("CheckOverlaps", false))
    , fMakeMCParticles(pset.get<bool>("MakeMCParticles", true))
    , fStoreDroppedMCParticles(pset.get<bool>("StoreDroppedMCParticles", false))
    , 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>()->registerAndSeedEngine(
        createEngine(0, "HepJamesRandom", "propagation"),
        "HepJamesRandom",
        "propagation",
        pset,
        "PropagationSeed"))
    , fDetProp{art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob()}
    , fAllPhysicsLists{fDetProp}
  {
    MF_LOG_DEBUG("LArG4") << "Debug: LArG4()";

    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>()->registerAndSeedEngine(
      createEngine(0, "G4Engine", "GEANT"), "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 (fStoreDroppedMCParticles) { produces<std::vector<sim::MCParticleLite>>(); }
    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";
  }

  //----------------------------------------------------------------------
  void LArG4::beginJob()
  {
    fG4Help = std::make_unique<g4b::G4Helper>(fG4MacroPath, fG4PhysListName);

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

    art::ServiceHandle<geo::Geometry const> geom;
    fG4Help->ConstructDetector(geom->GDMLFile());
    auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();

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

    LArVoxelReadoutGeometry::Setup_t readoutGeomSetupData{
      geom.get(), &wireReadoutGeom, lgp.get(), {&fEngine, fOffPlaneMargin}};

    fVoxelReadoutGeometry =
      new LArVoxelReadoutGeometry("LArVoxelReadoutGeometry", readoutGeomSetupData);
    pworlds.push_back(fVoxelReadoutGeometry);
    pworlds.push_back(new OpDetReadoutGeometry(
      geom->Cryostat().OpDetGeoName(), "OpDetReadoutGeometry", fUseLitePhotons));
    pworlds.push_back(new AuxDetReadoutGeometry(
      art::ServiceHandle<geo::AuxDetGeometry const>()->GetProviderPtr(), "AuxDetReadoutGeometry"));

    fG4Help->SetParallelWorlds(pworlds);

    // 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 particles and trajectories produced in the
    // detector.
    fparticleListAction = new larg4::ParticleListAction(lgp->ParticleKineticEnergyCut(),
                                                        lgp->StoreTrajectories(),
                                                        lgp->KeepEMShowerDaughters(),
                                                        fMakeMCParticles,
                                                        fStoreDroppedMCParticles);
    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::beginRun(art::Run& run)
  {
    // 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::CreateParticleVolumeFilter(
    std::set<std::string> const& vol_names) const
  {
    // 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::produce(art::Event& evt)
  {
    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 droppedPartCol =
      fStoreDroppedMCParticles ? std::make_unique<std::vector<sim::MCParticleLite>>() : 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;
    auto const* auxDetGeom = art::ServiceHandle<geo::AuxDetGeometry const>()->GetProviderPtr();

    // 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))
      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)

        if (fStoreDroppedMCParticles && droppedPartCol) {
          // Request a list of dropped particles. Store them in MCParticleLite format
          sim::ParticleList droppedParticleList = fparticleListAction->YieldDroppedList();
          droppedPartCol->reserve(droppedParticleList.size());

          for (auto const& partPair : droppedParticleList) {
            simb::MCParticle& p = *(partPair.second);
            if (ParticleListAction::isDropped(&p)) continue;
            if (p.StatusCode() != 1) continue;

            sim::MCParticleLite mini_mcp(p);
            mini_mcp.Origin(mct->Origin());

            droppedPartCol->push_back(std::move(mini_mcp));
          } // for(droppedParticleList)
        }

        // 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(geo::CryostatID(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,
                                                             ide.origTrackID);
                  } // 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(auxDetGeom->NAuxDets());
    for (unsigned int a = 0; a < auxDetGeom->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
      auto const& ad = auxDetGeom->AuxDet(a);
      for (size_t sv = 0; sv < ad.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 (droppedPartCol) {
      std::cout << "LArG4 dropped particles length = " << droppedPartCol->size() << std::endl;
      evt.put(std::move(droppedPartCol));
    }
    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()

} // namespace LArG4

DEFINE_ART_MODULE(larg4::LArG4)