OpMCDigi_module.cc

Go to the documentation of this file.

// \file OpMCDigi.h
// \author Ben Jones and Christie Chiu, MIT, Sept 2012
//   bjpjones@mit.edu, cschiu@mit.edu
//
// This module starts from MC truth sim::OnePhoton objects
// and produces a digitized waveform.
//
// It is assumed that the electronics response is linear,
// and the 1PE waveform can be described by a discreate
// response shape.  The many PE response is then the linear
// superposition of the relevant function at the appropriate
// arrival times.
//

// 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 "fhiclcpp/ParameterSet.h"

// LArSoft includes
#include "larana/OpticalDetector/OpDigiProperties.h"
#include "larana/OpticalDetector/OpDetResponseInterface.h"
#include "larsim/Simulation/SimListUtils.h"
#include "larsim/Simulation/LArG4Parameters.h"
#include "lardataobj/Simulation/SimPhotons.h"
#include "lardataobj/RawData/OpDetPulse.h"

// CLHEP includes
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandPoisson.h"

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

// C++ language includes
#include <cstring>

namespace opdet {

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

  private:
    void produce(art::Event&) override;

    // The parameters we'll read from the .fcl file.
    std::string fInputModule;              // Input tag for OpDet collection
    float fSampleFreq;                     // in MHz
    float fTimeBegin;                      // in us
    float fTimeEnd;                        // in us
    //float fQE;                             // quantum efficiency of opdet
    float fSaturationScale;                // adc count w/ saturation occurs

    float fDarkRate;                      // Noise rate in Hz

    std::vector<double> fSinglePEWaveform;

    CLHEP::HepRandomEngine& fEngine;
    CLHEP::RandFlat    fFlatRandom;
    CLHEP::RandPoisson fPoissonRandom;

    void AddTimedWaveform (int time, std::vector<double>& OldPulse, std::vector<double>& NewPulse);
  };
}

// Debug flag; only used during code development.
// const bool debug = true;

namespace opdet {


  OpMCDigi::OpMCDigi(fhicl::ParameterSet const& pset)
    : EDProducer{pset}
    , fInputModule{pset.get<std::string>("InputModule")}
      //, fQE{pset.get<double>("QE")}
    , fSaturationScale{pset.get<float>("SaturationScale")}
    , fDarkRate{pset.get<float>("DarkRate")}
      // create a default random engine; obtain the random seed from NuRandomService,
      // unless overridden in configuration with key "Seed"
    , fEngine(art::ServiceHandle<rndm::NuRandomService>{}->createEngine(*this, pset, "Seed"))
    , fFlatRandom{fEngine}
    , fPoissonRandom{fEngine}
  {
    produces<std::vector< raw::OpDetPulse> >();

    art::ServiceHandle<OpDigiProperties> odp;
    fSampleFreq = odp->SampleFreq();
    fTimeBegin  = odp->TimeBegin();
    fTimeEnd    = odp->TimeEnd();
    fSinglePEWaveform = odp->SinglePEWaveform();
  }

  //-------------------------------------------------


  void OpMCDigi::AddTimedWaveform (int binTime, std::vector<double>& OldPulse, std::vector<double>& NewPulse)
  {

    if( (binTime + NewPulse.size() ) > OldPulse.size()) {
      OldPulse.resize(binTime + NewPulse.size());
    }

    // Add shifted NewWaveform to Waveform at pointer
    for(size_t i = 0; i!=NewPulse.size(); ++i) {
      OldPulse.at(binTime+i) += NewPulse.at(i);
    }
  }


  //-------------------------------------------------

  void OpMCDigi::produce(art::Event& evt)
  {
    auto StoragePtr = std::make_unique<std::vector<raw::OpDetPulse>>();

    bool const fUseLitePhotons = art::ServiceHandle<sim::LArG4Parameters const>{}->UseLitePhotons();

    // Service for determining opdet responses
    art::ServiceHandle<opdet::OpDetResponseInterface const> odresponse;

    double const TimeBegin_ns  = fTimeBegin  *  1000;
    double const TimeEnd_ns    = fTimeEnd    *  1000;
    double const SampleFreq_ns = fSampleFreq /  1000;

    int const nSamples = ( TimeEnd_ns-TimeBegin_ns)*SampleFreq_ns;
    int const NOpChannels = odresponse->NOpChannels();


    // This vector will store all the waveforms we will make
    std::vector<std::vector<double> > PulsesFromDetPhotons(NOpChannels,std::vector<double>(nSamples,0.0));

    if(!fUseLitePhotons) {
      // Read in the Sim Photons
      sim::SimPhotonsCollection ThePhotCollection = sim::SimListUtils::GetSimPhotonsCollection(evt,fInputModule);
      // For every OpDet:
      for(auto const& pr : ThePhotCollection) {
        const sim::SimPhotons& ThePhot=pr.second;

        int const Ch = ThePhot.OpChannel();
        int readoutCh;

        // For every photon in the hit:
        for(const sim::OnePhoton& Phot: ThePhot) {
          // Sample a random subset according to QE
          if(!odresponse->detected(Ch, Phot, readoutCh)) {
            continue;
          }

          // Convert photon arrival time to the appropriate bin,
          // dictated by fSampleFreq. Photon arrival time is in ns,
          // beginning time in us, and sample frequency in MHz. Notice
          // that we have to accommodate for the beginning time
          if((Phot.Time > TimeBegin_ns) && (Phot.Time < TimeEnd_ns)) {
            auto const binTime = static_cast<int>((Phot.Time - TimeBegin_ns) * SampleFreq_ns);
            AddTimedWaveform( binTime, PulsesFromDetPhotons[readoutCh], fSinglePEWaveform );
          }
        } // for each Photon in SimPhotons
      }
    }
    else {
      auto const photons = *evt.getValidHandle<std::vector<sim::SimPhotonsLite>>("largeant");
      // For every OpDet:
      for (auto const& photon : photons) {
        int const Ch=photon.OpChannel;
        int readoutCh;

        std::map<int, int> PhotonsMap = photon.DetectedPhotons;

        // For every photon in the hit:
        for(auto const& pr : photon.DetectedPhotons) {
          for(int i = 0; i < pr.second; i++) {
            // Sample a random subset according to QE
            if(odresponse->detectedLite(Ch, readoutCh)) {
              // Convert photon arrival time to the appropriate bin, dictated by fSampleFreq.
              // Photon arrival time is in ns, beginning time in us, and sample frequency in MHz.
              // Notice that we have to accommodate for the beginning time
              if((pr.first > TimeBegin_ns) && (pr.first < TimeEnd_ns)) {
                auto const binTime = static_cast<int>((pr.first - TimeBegin_ns) * SampleFreq_ns);
                AddTimedWaveform( binTime, PulsesFromDetPhotons[readoutCh], fSinglePEWaveform );
              }
            } // random QE cut
          }
        } // for each Photon in SimPhotons
      }
    }

    // Create vector of output objects, add dark noise and apply
    //  saturation

    std::vector<raw::OpDetPulse*> ThePulses(NOpChannels);
    for(int iCh=0; iCh!=NOpChannels; ++iCh) {
      PulsesFromDetPhotons[iCh].resize((TimeEnd_ns - TimeBegin_ns) * SampleFreq_ns);

      // Add dark noise
      double const MeanDarkPulses = fDarkRate * (fTimeEnd-fTimeBegin) / 1000000;
      unsigned const int NumberOfPulses = fPoissonRandom.fire(MeanDarkPulses);

      for(size_t i=0; i!=NumberOfPulses; ++i) {
        double const PulseTime = (fTimeEnd-fTimeBegin)*fFlatRandom.fire(1.0);
        int const binTime = static_cast<int>(PulseTime * fSampleFreq);

        AddTimedWaveform(binTime, PulsesFromDetPhotons[iCh], fSinglePEWaveform);
      }

      // Apply saturation for large signals
      for(size_t i=0; i!=PulsesFromDetPhotons[iCh].size(); ++i) {
        if(PulsesFromDetPhotons[iCh].at(i)>fSaturationScale) PulsesFromDetPhotons[iCh].at(i) = fSaturationScale;
      }

      // Produce ADC pulse of integers rather than doubles

      std::vector<short> shortvec;

      for(size_t i=0; i!=PulsesFromDetPhotons[iCh].size(); ++i) {
        // Throw randoms to fairly sample +ve and -ve side of doubles
        int ThisSample = PulsesFromDetPhotons[iCh].at(i);
        if(ThisSample>0) {
          if(fFlatRandom.fire(1.0) > (ThisSample - int(ThisSample)))
            shortvec.push_back(int(ThisSample));
          else
            shortvec.push_back(int(ThisSample)+1);
        }
        else {
          if(fFlatRandom.fire(1.0) >  (int(ThisSample)-ThisSample))
            shortvec.push_back(int(ThisSample));
          else
            shortvec.push_back(int(ThisSample)-1);
        }
      }

      StoragePtr->emplace_back(iCh, shortvec ,0, fTimeBegin);

    } // for each OpDet in SimPhotonsCollection

    evt.put(std::move(StoragePtr));
  }
}

DEFINE_ART_MODULE(opdet::OpMCDigi)