SdPMTSimulatorOG/SdPMTSimulator.cc

Go to the documentation of this file.

#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <sevt/PMT.h>
#include <sevt/PMTSimData.h>
#include <utl/ErrorLogger.h>
#include <utl/RandomEngine.h>
#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>
#include <utl/Reader.h>
#include <utl/Accumulator.h>
#include <utl/SaveCurrentTDirectory.h>
#include <utl/config.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <sdet/PMT.h>
#include <sdet/PMTConstants.h>

#include "SdPMTSimulator.h"

using namespace SdPMTSimulatorOG;
using namespace std;
using namespace utl;
using namespace fwk;
using namespace sevt;

using CLHEP::RandGeneral;


VModule::ResultFlag
SdPMTSimulator::Init()
{
  auto topB = CentralConfig::GetInstance()->GetTopBranch("SdPMTSimulator");
  vector<int> components;
  if (topB.GetChild("fundamentalComponents"))
    topB.GetChild("fundamentalComponents").GetData(components);
  for (const auto c : components) {
    if (StationConstants::eTotal <= c && c <= StationConstants::eLastSource)
      fFundamentalComponents.insert(StationConstants::SignalComponent(c));
    else {
      ostringstream err;
      err << "Invalid signal component: " << c;
      ERROR(err);
      throw InvalidConfigurationException(err.str());
    }
  }
  if (fFundamentalComponents.empty())
    INFO("Will not use any fundamental components.");
  else {
    ostringstream info;
    info << "Will use " << fFundamentalComponents.size() << " fundamental components:";
    for (const auto c : fFundamentalComponents)
      info << ' ' << c << "='" << GetSignalComponentName(c) << '\'';
    INFO(info);
  }

  auto elecB = topB.GetChild("electronicsParameters");

  elecB.GetChild("kFactor").GetData(fKFactor);

  fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

  // Read the voltage values for PMT saturation
  const AttributeMap useAtt = { { "use", "yes" } };

  auto satB = elecB.GetChild("saturationCurve", useAtt);
  if (satB) {
    auto vIn = satB.GetChild("inputVoltage").Get<vector<double>>();
    auto vOut = satB.GetChild("outputVoltage").Get<vector<double>>();
    satB.GetChild("Current2VoltageMultiplier").GetData(fCurrent2VoltageMultiplier);
    if (vIn.size() != vOut.size()) {
      ERROR("Lists of input/output voltages for saturation curve must have same size");
      return eFailure;
    }
    ostringstream info;
    info << "saturation curve [V]: ";
    for (unsigned int i = 0, n = vIn.size(); i < n; ++i)
      info << vIn[i] / volt << " -> " << vOut[i] / volt << "   ";
    INFO(info);
    // convert voltages to "current" (i.e. smeared&amplified pe per ns bin):
    const double v2c = 1 / fCurrent2VoltageMultiplier;
    for (unsigned int i = 0, n = vIn.size(); i < n; ++i) {
      vIn[i] *= v2c;
      vOut[i] *= v2c;
    }
    // input/output of this function is "current" (i.e. pe per ns bin)
    fSaturationFunction = new TabulatedFunction(vIn, vOut);
  } else
    INFO("No PMT saturation simulation");

  return eSuccess;
}


void
SdPMTSimulator::FillPulseShape(const utl::TabulatedFunction& rawPulse)
{
  const unsigned int n = rawPulse.GetNPoints();
  if (n < 2)
    return;

  /* The raw, measured pulse (accessed through SDetector PMT) is
  normalized and time=0 is set to its start. This process is
  performed by iterating over the raw pulse. No assumptions are
  made as to whether the measurement was performed in equally
  spaced time intervals or as to the polarity (negative or
  positive pulse).
  During the iteration, we find the smallest interval of
  time between measurements. We will use one order of magnitude
  lower than this for binning during integration. */

  const double start = rawPulse.XFront();
  const double stop = rawPulse.XBack();

  utl::Accumulator::Min<double> minTimeDifference(stop - start);
  for (unsigned int i = 1; i < n; ++i)
    minTimeDifference(rawPulse[i].X() - rawPulse[i-1].X());

  const double timeBinning = minTimeDifference.GetMin() / 10;

  double sum = 0;
  for (double t = start; t < stop; t += timeBinning)
    sum += rawPulse.InterpolateY(t, 1);

  const double norm = sum * timeBinning;

  // During the third iteration, we perform the normalization and
  // appropriately adjust the times as we fill fMeanPulseShape.
  fMeanPulseShape.Clear();
  for (const auto& xy : rawPulse)
    fMeanPulseShape.PushBack(xy.X() - start, xy.Y() / norm);

  // We iterate over the newly filled tabulated function once to
  // save the point at which the pulse is at a maximumum.
  double currentMax = 0;
  for (const auto& th : fMeanPulseShape) {
    const double h = th.Y();
    if (h > currentMax) {
      currentMax = h;
      fPulseMaxTime = th.X();
    }
  }
}



void
SdPMTSimulator::ConvertPEToBaseSignal(const double scaling,
                                      const TimeDistributionI& pe,
                                      TimeDistributionD& base)
{
  if (!pe.GetBinning())
    return;

  const double dt = pe.GetBinning();
  const double pulseStopTime = fMeanPulseShape.XBack();
  vector<double> meanPulseShape;
  meanPulseShape.reserve(int(pulseStopTime / dt) + 1);
  double meanPulseShapeNorm = 0;

  for (double time = 0; time < pulseStopTime; time += dt) {
    const double pulse = fMeanPulseShape.Y(time);
    meanPulseShapeNorm += pulse;
    meanPulseShape.push_back(pulse);
  }

  const double normGain = fPMTGain / meanPulseShapeNorm;
  const double scaleFactorMin = normGain * fMinCharge;
  // DV: instead of changing mu = scaling * tn.get<1>() in the for-loop below
  // to a Poi(mu) we opt for scale * charge (effectively via scaleFactorMax) so
  // that no additional fluctuations are introduced to already increased
  // fluctuations due to the station-level thinning
  const double scaleFactorMax = normGain * fMaxCharge * scaling;

  for (const auto& tn : pe.SparseRange()) {

    double charge = 0;
    for (int i = 0, n = tn.get<1>(); i < n; ++i)
      charge += fChargeDistribution->fire();

    const double realCharge = scaleFactorMin + scaleFactorMax * charge;

    // note: get<0>() is a slot index while offset has to be in ns
    const double offset = tn.get<0>() * dt - fPulseMaxTime;

    // place scaled pulse shape centered with max at current time bin
    int i = -1;
    for (double time = 0; time <= pulseStopTime; time += dt) {
      // center pulse max on pe release time.
      // (may not be the most accurate thing to do ...)
      base.AddTime(time + offset, realCharge * meanPulseShape[++i]);
    }

  }
}


void
SdPMTSimulator::SimulateSaturation(TimeDistributionD& baseSignal,
                                   const TimeDistributionD& totalBaseSignal)
{
  for (const auto& ts : baseSignal.SparseRange()) {
    const auto t = ts.get<0>();
    auto& sig = baseSignal[t]; // op[] is the only way to get a reference
    const double tot = totalBaseSignal.At(t);
    const double maxDist = max(sig, tot);
    const double scaling = (maxDist > 0) ? fSaturationFunction->Y(maxDist) / maxDist : 1;
    sig *= scaling;
  }
}


VModule::ResultFlag
SdPMTSimulator::Run(evt::Event& event)
{
  if (!event.HasSEvent()) {
    ERROR("Event does not have an SEvent");
    return eFailure;
  }

  const bool haveFund = !fFundamentalComponents.empty();

  auto& sEvent = event.GetSEvent();

  for (auto& station : sEvent.StationsRange()) {

    if (!station.HasSimData())
      continue;

    auto& sSim = station.GetSimData();

    const auto& simulatorSignature = sSim.GetSimulatorSignature();

    // fKFactor is the relative ratio of (QE*CE) for small PMT and standard PMT
    // it must be applied to Small PMT signal in case of standard killing of photons
    // (G4StationSimulatorFullOG)
    // no kFactor needed for G4StationSimulatorFullTrackOG

    const double kFactor =
      (simulatorSignature == "G4StationSimulatorFullTrackOG") ? 1 : fKFactor;

    const auto& dStation = det::Detector::GetInstance().GetSDetector().GetStation(station);

    const bool isUUB = dStation.IsUUB();

    // scale PE up if the particles underwent Station-based thinning
    const double scaling =
      (sSim.GetThinning() == 1) ?
      1 : double(sSim.GetTotalSimCandidateParticleCount()) / sSim.GetTotalSimParticleCount();

    for (auto& pmt : station.PMTsRange(sdet::PMTConstants::eAnyType)) {

      if (!pmt.HasSimData())
        continue;

      auto& pSim = pmt.GetSimData();
      if (!pSim.HasPETimeDistribution() || pSim.HasBaseSignal())
        continue;

      // Getting efficiencies gain pulse shape and charge measurements from SDetector PMT
      const auto& dPMT = dStation.GetPMT(pmt.GetId());

      fPMTGain = dPMT.GetWorkingGain();

      if (isUUB && dPMT.GetType() == sdet::PMTConstants::eWaterCherenkovSmall)
        fPMTGain *= kFactor;

      const auto& rawPulse = dPMT.GetPulseShape();
      const auto& rawChargeProbDist = dPMT.GetChargeProbabilityDistribution();
      fMinCharge = dPMT.GetMinCharge();
      fMaxCharge = dPMT.GetMaxCharge();

      // Filling pulse shape and charge probability distribution
      FillPulseShape(rawPulse);
      delete fChargeDistribution;
      fChargeDistribution =
        new RandGeneral(fRandomEngine->GetEngine(), &rawChargeProbDist.front(), rawChargeProbDist.size());

      pSim.MakeBaseSignal();
      auto& totalSignal = pSim.GetBaseSignal();

      for (const auto peComp : pSim.PETimeDistributionsRange()) {

        const auto comp = StationConstants::SignalComponent(peComp.GetLabel());

        if (comp == StationConstants::eTotalNoSaturation || // later
            (haveFund && comp == StationConstants::eTotal))
          continue;

        if (!pSim.HasBaseSignal(comp))
          pSim.MakeBaseSignal(comp);
        auto& base = pSim.GetBaseSignal(comp);

        ConvertPEToBaseSignal(scaling, peComp.GetTimeDistribution(), base);

        if (haveFund && fFundamentalComponents.count(comp))
          totalSignal += base;

      }

      // eTotalNoSaturation made out of eTotal before saturation
      if (!pSim.HasBaseSignal(StationConstants::eTotalNoSaturation))
        pSim.MakeBaseSignal(StationConstants::eTotalNoSaturation);
      auto& totalNoSaturationSignal = pSim.GetBaseSignal(StationConstants::eTotalNoSaturation);
      totalNoSaturationSignal = totalSignal;

      if (fSaturationFunction)
        SimulateSaturation(totalSignal, totalNoSaturationSignal);

    }

    sSim.SetUsedWeight(scaling);

  }

  return eSuccess;
}


VModule::ResultFlag
SdPMTSimulator::Finish()
{
  delete fChargeDistribution;
  fChargeDistribution = nullptr;
  delete fSaturationFunction;
  fSaturationFunction = nullptr;

  return eSuccess;
}