SdBaselineSimulatorOG/SdBaselineSimulator.cc

Go to the documentation of this file.

#include "SdBaselineSimulator.h"

#include <evt/Event.h>

#include <det/Detector.h>

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

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>
#include <sevt/PMT.h>
#include <sevt/PMTSimData.h>

#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>

#include <utl/RandomEngine.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/Accumulator.h>

#include <CLHEP/Random/Randomize.h>

#include <sstream>

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


namespace SdBaselineSimulatorOG {

  VModule::ResultFlag
  SdBaselineSimulator::Init()
  {
    fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);
    return eSuccess;
  }


  VModule::ResultFlag
  SdBaselineSimulator::Run(evt::Event& event)
  {
    const auto rnd = &fRandomEngine->GetEngine();

    auto& sEvent = event.GetSEvent();

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

      // skip station if no particles
      if (!station.HasSimData())
        continue;

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

      // make sure all traces in a station start and end exactly at the same times
      Accumulator::Safe<Accumulator::Min<int>> minTime;
      Accumulator::Safe<Accumulator::Max<int>> maxTime;
      for (auto& pmt : station.PMTsRange(sdet::PMTConstants::eAnyType)) {
        if (!pmt.HasSimData())
          continue;
        const auto& pmtSim = pmt.GetSimData();
        if (!pmtSim.HasFADCTrace())
          continue;
        const auto& trace = pmtSim.GetFADCTraceD(sdet::PMTConstants::eHighGain, StationConstants::eTotal);
        minTime(trace.GetStart());
        maxTime(trace.GetStop());
      }

      const auto saturationValue = dStation.GetSaturationValue();

      // add padding before and after
      const int startTime = minTime.GetMin() - traceLength/3;  // this should be ok also for UUB
      const int stopTime = maxTime.GetMax() + traceLength;

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

        if (!pmt.HasSimData())
          continue;
        auto& pmtSim = pmt.GetSimData();

        if (!pmtSim.HasFADCTrace())
          continue;

        if (!pmtSim.HasFADCTrace(StationConstants::eTotalNoSaturation))
          pmtSim.MakeFADCTrace(StationConstants::eTotalNoSaturation);

        if (!pmt.HasCalibData())
          pmt.MakeCalibData();
        auto& pmtCalib = pmt.GetCalibData();

        const auto& dPMT = dStation.GetPMT(pmt);

        for (auto gain : { sdet::PMTConstants::eHighGain, sdet::PMTConstants::eLowGain }) {
          const double nominalBaseline = dPMT.GetBaseline(gain);
          // choose uniform random phase between the adc counts
          const double baseline = CLHEP::RandFlat::shoot(rnd, nominalBaseline-0.5, nominalBaseline+0.5);
          const double noise = dPMT.GetBaselineRMS(gain);
          pmtCalib.SetBaseline(baseline, noise, gain);

          auto& traceD = pmtSim.GetFADCTraceD(gain, StationConstants::eTotal);
          auto& noSat = pmtSim.GetFADCTrace(gain, StationConstants::eTotalNoSaturation);
          auto& trace = pmtSim.GetFADCTrace(gain, StationConstants::eTotal);
          for (int i = startTime; i <= stopTime; ++i) {
            const double noisyBaseline = CLHEP::RandGauss::shoot(rnd, baseline, noise);
            const double signal = traceD.At(i) + noisyBaseline;
            const double saturatedSignal = min(max(signal, 0.), double(saturationValue));
            traceD.SetTime(i, saturatedSignal);
            const int signalI = round(signal);
            const int positiveSignalI = max(signalI, 0);
            noSat.SetTime(i, positiveSignalI);
            const int saturatedSignalI = min(positiveSignalI, int(saturationValue));
            trace.SetTime(i, saturatedSignalI);
          }
        }

      }

    }

    return eSuccess;
  }

}