PDFastSimPVS_module.cc

Go to the documentation of this file.

// Class:       PDFastSimPVS
// Plugin Type: producer
// File:        PDFastSimPVS_module.cc
// Description:
// - acts on sim::SimEnergyDeposit from LArG4Main,
// - simulate (fast, photon visibility service) the OpDet response to optical photons
// Input: 'sim::SimEnergyDeposit'
// Output: 'sim::OpDetBacktrackerRecord'
//Fast simulation of propagating the photons created from SimEnergyDeposits.

//This module does a fast simulation of propagating the photons created from SimEnergyDeposits,
//This simulation is done using the PhotonLibrary, which stores the visibilities of each optical channel
//with respect to each optical voxel in the TPC volume, to avoid propagating single photons using Geant4.
//At the end of this module a collection of the propagated photons either as
//'sim::OpDetBacktrackerRecord' are placed into the art event.

//The steps this module takes are:
//  - to take number of photon and the vertex information from 'sim::SimEnergyDeposits',
//  - use the PhotonLibrary (visibilities) to determine the amount of visible photons at each optical channel,
//  - visible photons: the number of photons times the visibility at the middle of the Geant4 step for a given optical channel.
//  - other photon information is got from 'sim::SimEnergyDeposits'
//  - add 'sim::OpDetBacktrackerRecord' to event
// Aug. 19 by Mu Wei
// Restructured Nov. 21 by P. Green

// LArSoft libraries
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "lardataobj/Simulation/OpDetBacktrackerRecord.h"
#include "lardataobj/Simulation/SimEnergyDeposit.h"
#include "lardataobj/Simulation/SimPhotons.h"
#include "larsim/PhotonPropagation/PhotonVisibilityService.h"
#include "larsim/PhotonPropagation/PropagationTimeModel.h"
#include "larsim/PhotonPropagation/ScintTimeTools/ScintTime.h"
#include "larsim/Simulation/LArG4Parameters.h"

#include "nurandom/RandomUtils/NuRandomService.h"

// Art libraries
#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/Utilities/make_tool.h"
#include "canvas/Utilities/Exception.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "CLHEP/Random/RandPoissonQ.h"
#include "CLHEP/Units/SystemOfUnits.h"

#include <iostream>
#include <memory>
#include <utility>
#include <vector>

namespace phot {
  class PDFastSimPVS : public art::EDProducer {
  public:
    explicit PDFastSimPVS(fhicl::ParameterSet const&);
    void produce(art::Event&) override;
    void AddOpDetBTR(std::vector<sim::OpDetBacktrackerRecord>& opbtr,
                     std::vector<int>& ChannelMap,
                     const sim::OpDetBacktrackerRecord& btr) const;

  private:
    art::ServiceHandle<PhotonVisibilityService const> fPVS;
    const bool fDoFastComponent;
    const bool fDoSlowComponent;
    const bool fIncludePropTime;
    const bool fUseLitePhotons;
    const bool fStoreReflected;
    const size_t fNOpChannels;
    const art::InputTag fSimTag;
    CLHEP::HepRandomEngine& fPhotonEngine;
    std::unique_ptr<CLHEP::RandPoissonQ> fRandPoissPhot;
    CLHEP::HepRandomEngine& fScintTimeEngine;
    std::unique_ptr<ScintTime> fScintTime; // Tool to retrive timing of scintillation
    std::unique_ptr<PropagationTimeModel> fPropTimeModel;
  };

  //......................................................................
  PDFastSimPVS::PDFastSimPVS(fhicl::ParameterSet const& pset)
    : art::EDProducer{pset}
    , fPVS(art::ServiceHandle<PhotonVisibilityService const>())
    , fDoFastComponent(pset.get<bool>("DoFastComponent", true))
    , fDoSlowComponent(pset.get<bool>("DoSlowComponent", true))
    , fIncludePropTime(pset.get<bool>("IncludePropTime", false))
    , fUseLitePhotons(art::ServiceHandle<sim::LArG4Parameters const>()->UseLitePhotons())
    , fStoreReflected(fPVS->StoreReflected())
    , fNOpChannels(fPVS->NOpChannels())
    , fSimTag{pset.get<art::InputTag>("SimulationLabel")}
    , fPhotonEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(
        createEngine(0, "HepJamesRandom", "photon"),
        "HepJamesRandom",
        "photon",
        pset,
        "SeedPhoton"))
    , fRandPoissPhot(std::make_unique<CLHEP::RandPoissonQ>(fPhotonEngine))
    , fScintTimeEngine(art::ServiceHandle<rndm::NuRandomService>()->registerAndSeedEngine(
        createEngine(0, "HepJamesRandom", "scinttime"),
        "HepJamesRandom",
        "scinttime",
        pset,
        "SeedScintTime"))
    , fScintTime{art::make_tool<phot::ScintTime>(pset.get<fhicl::ParameterSet>("ScintTimeTool"))}
  {
    mf::LogInfo("PDFastSimPVS") << "Initializing PDFastSimPVS." << std::endl;
    fhicl::ParameterSet VUVTimingParams;
    fhicl::ParameterSet VISTimingParams;
    // Validate configuration options
    if (fIncludePropTime &&
        !pset.get_if_present<fhicl::ParameterSet>("VUVTiming", VUVTimingParams)) {
      throw art::Exception(art::errors::Configuration)
        << "Propagation time simulation requested, but VUVTiming not specified."
        << "\n";
    }
    if (fIncludePropTime && fStoreReflected &&
        !pset.get_if_present<fhicl::ParameterSet>("VISTiming", VISTimingParams)) {
      throw art::Exception(art::errors::Configuration)
        << "Reflected light propagation time simulation requested, but VISTiming not specified."
        << "\n";
    }
    // Initialise the Scintillation Time
    fScintTime->initRand(fScintTimeEngine);

    // propagation time model
    if (fIncludePropTime)
      fPropTimeModel = std::make_unique<PropagationTimeModel>(
        VUVTimingParams, VISTimingParams, fScintTimeEngine, fStoreReflected);

    if (fUseLitePhotons) {
      mf::LogInfo("PDFastSimPVS") << "Use Lite Photon." << std::endl;
      produces<std::vector<sim::SimPhotonsLite>>();
      produces<std::vector<sim::OpDetBacktrackerRecord>>();

      if (fStoreReflected) {
        mf::LogInfo("PDFastSimPVS") << "Store Reflected Photons";
        produces<std::vector<sim::SimPhotonsLite>>("Reflected");
        produces<std::vector<sim::OpDetBacktrackerRecord>>("Reflected");
      }
    }
    else {
      mf::LogInfo("PDFastSimPVS") << "Use Sim Photon.";
      produces<std::vector<sim::SimPhotons>>();
      if (fStoreReflected) {
        mf::LogInfo("PDFastSimPVS") << "Store Reflected Photons";
        produces<std::vector<sim::SimPhotons>>("Reflected");
      }
    }
    mf::LogInfo("PDFastSimPVS") << "PDFastSimPVS Initialization finish" << std::endl;
  }

  //......................................................................
  void PDFastSimPVS::produce(art::Event& event)
  {
    std::vector<int> PDChannelToSOCMapDirect(fNOpChannels, -1);  // Where each OpChan is.
    std::vector<int> PDChannelToSOCMapReflect(fNOpChannels, -1); // Where each OpChan is.

    // SimPhotonsLite
    auto phlit = std::make_unique<std::vector<sim::SimPhotonsLite>>();
    auto opbtr = std::make_unique<std::vector<sim::OpDetBacktrackerRecord>>();
    auto phlit_ref = std::make_unique<std::vector<sim::SimPhotonsLite>>();
    auto opbtr_ref = std::make_unique<std::vector<sim::OpDetBacktrackerRecord>>();
    auto& dir_phlitcol(*phlit);
    auto& ref_phlitcol(*phlit_ref);
    // SimPhotons
    auto phot = std::make_unique<std::vector<sim::SimPhotons>>();
    auto phot_ref = std::make_unique<std::vector<sim::SimPhotons>>();
    auto& dir_photcol(*phot);
    auto& ref_photcol(*phot_ref);
    if (fUseLitePhotons) {
      dir_phlitcol.resize(fNOpChannels);
      ref_phlitcol.resize(fNOpChannels);
      for (unsigned int i = 0; i < fNOpChannels; i++) {
        dir_phlitcol[i].OpChannel = i;
        ref_phlitcol[i].OpChannel = i;
      }
    }
    else { // SimPhotons
      dir_photcol.resize(fNOpChannels);
      ref_photcol.resize(fNOpChannels);
      for (unsigned int i = 0; i < fNOpChannels; i++) {
        dir_photcol[i].fOpChannel = i;
        ref_photcol[i].fOpChannel = i;
      }
    }

    art::Handle<std::vector<sim::SimEnergyDeposit>> edepHandle;
    if (!event.getByLabel(fSimTag, edepHandle)) {
      mf::LogWarning("PDFastSimPVS") << "PDFastSimPVS Module Cannot getByLabel: " << fSimTag;
      return;
    }

    int num_points = 0;
    auto const& edeps = edepHandle;
    for (auto const& edepi : *edeps) {
      num_points++;

      int nphot_fast = edepi.NumFPhotons();
      int nphot_slow = edepi.NumSPhotons();
      if (!((nphot_fast > 0 && fDoFastComponent) || (nphot_slow > 0 && fDoSlowComponent))) continue;

      auto const& prt = edepi.MidPoint();

      MappedCounts_t Visibilities = fPVS->GetAllVisibilities(prt);
      MappedCounts_t Visibilities_Ref;
      if (fStoreReflected) {
        Visibilities_Ref = fPVS->GetAllVisibilities(prt, true);
        if (!Visibilities_Ref)
          mf::LogWarning("PDFastSimPVS") << "Fail to get visibilities for reflected photons.";
      }

      if (!Visibilities) {
        //throw cet::exception("PDFastSimPVS")
        mf::LogWarning("PDFastSimPVS")
          << "There is no entry in the PhotonLibrary for this position in space. Position: "
          << edepi.MidPoint() << "\n Move to next point";
        continue;
      }

      int trackID = edepi.TrackID();
      int nphot = edepi.NumPhotons();
      double edeposit = edepi.Energy() / nphot;
      double pos[3] = {edepi.MidPointX(), edepi.MidPointY(), edepi.MidPointZ()};

      geo::Point_t const ScintPoint = {pos[0], pos[1], pos[2]};

      // loop through direct photons then reflected photons cases
      size_t DoReflected = (fStoreReflected) ? 1 : 0;
      for (size_t Reflected = 0; Reflected <= DoReflected; ++Reflected) {
        for (unsigned int channel = 0; channel < fNOpChannels; ++channel) {

          double visibleFraction = (Reflected) ? Visibilities_Ref[channel] : Visibilities[channel];
          if (visibleFraction < 1e-9) continue; // voxel is not visible at this optical channel

          // number of detected photons
          int ndetected_fast =
            (nphot_fast > 0) ? fRandPoissPhot->fire(nphot_fast * visibleFraction) : 0;
          int ndetected_slow =
            (nphot_slow > 0) ? fRandPoissPhot->fire(nphot_slow * visibleFraction) : 0;
          if (!((ndetected_fast > 0 && fDoFastComponent) ||
                (ndetected_slow > 0 && fDoSlowComponent)))
            continue;

          // calculate propagation times if included, does not matter whether fast or slow photon
          std::vector<double> transport_time;
          if (fIncludePropTime) {
            transport_time.resize(ndetected_fast + ndetected_slow);
            fPropTimeModel->propagationTime(transport_time, ScintPoint, channel, Reflected);
          }

          // SimPhotonsLite case
          if (fUseLitePhotons) {
            sim::OpDetBacktrackerRecord tmpbtr(channel);
            if (ndetected_fast > 0 && fDoFastComponent) {
              for (int i = 0; i < ndetected_fast; ++i) {
                // calculate the time at which each photon is seen
                double dtime = edepi.StartT() + fScintTime->fastScintTime();
                if (fIncludePropTime) dtime += transport_time[i];
                int time(dtime > static_cast<double>(std::numeric_limits<int>::max()) ?
                           std::numeric_limits<int>::max() :
                           static_cast<int>(std::round(dtime)));
                if (Reflected)
                  ++ref_phlitcol[channel].DetectedPhotons[time];
                else
                  ++dir_phlitcol[channel].DetectedPhotons[time];
                tmpbtr.AddScintillationPhotons(trackID, time, 1, pos, edeposit);
              }
            }
            if ((ndetected_slow > 0) && fDoSlowComponent) {
              for (int i = 0; i < ndetected_slow; ++i) {
                // calculate the time at which each photon is seen
                double dtime = edepi.StartT() + fScintTime->slowScintTime();
                if (fIncludePropTime) dtime += transport_time[ndetected_fast + i];
                int time(dtime > static_cast<double>(std::numeric_limits<int>::max()) ?
                           std::numeric_limits<int>::max() :
                           static_cast<int>(std::round(dtime)));
                if (Reflected)
                  ++ref_phlitcol[channel].DetectedPhotons[time];
                else
                  ++dir_phlitcol[channel].DetectedPhotons[time];
                tmpbtr.AddScintillationPhotons(trackID, time, 1, pos, edeposit);
              }
            }
            if (Reflected)
              AddOpDetBTR(*opbtr_ref, PDChannelToSOCMapReflect, tmpbtr);
            else
              AddOpDetBTR(*opbtr, PDChannelToSOCMapDirect, tmpbtr);
          }

          // SimPhotons case
          else {
            sim::OnePhoton photon;
            photon.SetInSD = false;
            photon.InitialPosition = edepi.End();
            photon.MotherTrackID = edepi.TrackID();
            if (Reflected)
              photon.Energy = 2.9 * CLHEP::eV; // 430 nm
            else
              photon.Energy = 9.7 * CLHEP::eV; // 128 nm
            // TODO: un-hardcode and add another energy for Xe scintillation
            if (ndetected_fast > 0 && fDoFastComponent) {
              for (int i = 0; i < ndetected_fast; ++i) {
                // calculates the time at which the photon was produced
                double dtime = edepi.StartT() + fScintTime->fastScintTime();
                if (fIncludePropTime) dtime += transport_time[i];
                int time = static_cast<int>(std::round(dtime));
                photon.Time = time;
                if (Reflected)
                  ref_photcol[channel].insert(ref_photcol[channel].end(), 1, photon);
                else
                  dir_photcol[channel].insert(dir_photcol[channel].end(), 1, photon);
              }
            }
            if (ndetected_slow > 0 && fDoSlowComponent) {
              for (int i = 0; i < ndetected_slow; ++i) {
                double dtime = edepi.StartT() + fScintTime->slowScintTime();
                if (fIncludePropTime) dtime += transport_time[ndetected_fast + i];
                int time = static_cast<int>(std::round(dtime));
                photon.Time = time;
                if (Reflected)
                  ref_photcol[channel].insert(ref_photcol[channel].end(), 1, photon);
                else
                  dir_photcol[channel].insert(dir_photcol[channel].end(), 1, photon);
              }
            }
          }
        }
      }
    }

    if (fUseLitePhotons) {
      event.put(move(phlit));
      event.put(move(opbtr));
      if (fStoreReflected) {
        event.put(move(phlit_ref), "Reflected");
        event.put(move(opbtr_ref), "Reflected");
      }
    }
    else {
      event.put(move(phot));
      if (fStoreReflected) { event.put(move(phot_ref), "Reflected"); }
    }

    return;
  }

  //......................................................................
  void PDFastSimPVS::AddOpDetBTR(std::vector<sim::OpDetBacktrackerRecord>& opbtr,
                                 std::vector<int>& ChannelMap,
                                 const sim::OpDetBacktrackerRecord& btr) const
  {
    int iChan = btr.OpDetNum();
    if (ChannelMap[iChan] < 0) {
      ChannelMap[iChan] = opbtr.size();
      opbtr.emplace_back(std::move(btr));
    }
    else {
      size_t idtest = ChannelMap[iChan];
      auto const& timePDclockSDPsMap = btr.timePDclockSDPsMap();
      for (auto const& timePDclockSDP : timePDclockSDPsMap) {
        for (auto const& sdp : timePDclockSDP.second) {
          double xyz[3] = {sdp.x, sdp.y, sdp.z};
          opbtr.at(idtest).AddScintillationPhotons(
            sdp.trackID, timePDclockSDP.first, sdp.numPhotons, xyz, sdp.energy);
        }
      }
    }
  }

} // namespace

DEFINE_ART_MODULE(phot::PDFastSimPVS)