Deprecated/UpgradeASCIITests/SdPMTSimulatorASCII/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/config.h>

// Root stuff for diagnostic histograms
#include <TH1.h>
#include <TFile.h>

#include "SdPMTSimulator.h"

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

using CLHEP::RandGeneral;


SdPMTSimulator::SdPMTSimulator() :
  fChargeDist(0),
  fLimitStationsPerCycle(false),
  fNStationsPerCycle(numeric_limits<unsigned int>::max()),
  fAveragePulseShape(0),
  fWantRootHistos(false),
  fChargeHisto(0),
  fAnodePulseHisto(0),
  fSaturationFunction(0),
  fRandomEngine(0)
{
}

SdPMTSimulator::~SdPMTSimulator()
{
  delete fChargeDist;
  delete fAveragePulseShape;
  delete fChargeHisto;
  delete fAnodePulseHisto;
  delete fSaturationFunction;
}


VModule::ResultFlag
SdPMTSimulator::Init()
{
  Branch topB = CentralConfig::GetInstance()->GetTopBranch("SdPMTSimulatorASCII");
  vector<int> components;

  // WCD 
  if (topB.GetChild("fundamentalComponents"))
    topB.GetChild("fundamentalComponents").GetData(components);
  for (unsigned int i = 0; i != components.size(); ++i)
    fFundamentalComponents.push_back(StationConstants::SignalComponent(components[i]));
  
  // ASCII
  components.clear();
  if (topB.GetChild("fundamentalComponents_ASCII"))
    topB.GetChild("fundamentalComponents_ASCII").GetData(components);
  for (unsigned int i = 0; i != components.size(); ++i)
    fFundamentalComponents_ASCII.push_back(StationConstants::SignalComponent(components[i]));

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

  // Set up the average anode pulse shape from data in file
  Branch anodePulseBranch = elecB.GetChild("anodePulseShape");

//#warning XML for this two quantities should be replaced by the tabulated function
  vector<double> anodeTimeBins;
  anodePulseBranch.GetChild("timeBins").GetData(anodeTimeBins);
  vector<double> anodePulseHeight;
  anodePulseBranch.GetChild("pulseSize").GetData(anodePulseHeight);

  if (anodeTimeBins.empty() || anodeTimeBins.size() != anodePulseHeight.size()) {
    ERROR("Malformed <anodePulseShape> configuration.");
    return eFailure;
  }

  {
    const unsigned int naph = anodePulseHeight.size();

    // invert the raw pulse shape from scope trace and ignore below baseline
    for (unsigned int i = 0; i < naph; ++i)
      anodePulseHeight[i] = (anodePulseHeight[i] >= 0) ? 0 : -anodePulseHeight[i];

    double anodeIntegral = 0;
    for (unsigned int i1 = 0, i2 = 1; i2 < naph; ++i1, ++i2) {
      const double dt = anodeTimeBins[i2] - anodeTimeBins[i1];
      anodeIntegral += 0.5*(anodePulseHeight[i1] + anodePulseHeight[i2])*dt;
    }

    // normalize the pulse shape
    for (unsigned int i1 = 0, i2 = 1; i2 < naph; ++i1, ++i2) {
      const double dt = 1/*anodeTimeBins[i2] - anodeTimeBins[i1]*/;
      anodePulseHeight[i1] *= dt / anodeIntegral;
    }

    // shift the times so that they start at zero
    const double offset = anodeTimeBins[0];
    for (unsigned int i = 0; i < naph; ++i)
      anodeTimeBins[i] -= offset;
  }

  fAveragePulseShape = new TabulatedFunction(anodeTimeBins, anodePulseHeight);

  // Find the maximum pulse height and corresponding bin. This is used later
  // to center the pulse maximum on the PE release time bin
  double currentMax = 0;
  for (TabulatedFunction::Iterator fIt = fAveragePulseShape->Begin();
       fIt != fAveragePulseShape->End(); ++fIt)
    if (fIt->Y() >= currentMax) {
      currentMax = fIt->Y();
      fMaxPulseHeight = fIt;
    }

  // Set up the SPE charge distribution
  Branch chargeB = elecB.GetChild("chargeDistributionSPE");

  chargeB.GetChild("minCharge").GetData(fMinCharge);
  chargeB.GetChild("maxCharge").GetData(fMaxCharge);

  vector<double> chargeProbability;
  chargeB.GetChild("relativeProbability").GetData(chargeProbability);

  const unsigned int chargeProbabilityN = chargeProbability.size();

  fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);
  fChargeDist = new RandGeneral(fRandomEngine->GetEngine(), &chargeProbability[0], chargeProbabilityN);

  elecB.GetChild("PMTGainWCD").GetData(fPMTGain_WCD);
  elecB.GetChild("PMTGainASCII").GetData(fPMTGain_ASCII);
    
  topB.GetChild("RootHistogramControl").GetData(fWantRootHistos);

  // Fill diagnostic histograms, if requested
  if (fWantRootHistos) {
    topB.GetChild("RootHistogramFilename").GetData(fRootHistoFilename);
    INFO("Will write diagnostic histograms for this module");

    // check SPE distribution
    fChargeHisto = new TH1D("chargeDist", "charge distribution", 20,0,4);
    for (int itest = 0; itest < 100000; ++itest) {
      const double x = fChargeDist->fire();
      const double val = fMinCharge + (fMaxCharge - fMinCharge)*x;
      fChargeHisto->Fill(val);
    }

    // check average pulse shape
    const int nPoints = fAveragePulseShape->GetNPoints();
    fAnodePulseHisto = new TH1D("anodePulseHisto", "anode pulse histo",
                                nPoints,
                                (*fAveragePulseShape)[0].X()/nanosecond,
                                (*fAveragePulseShape)[nPoints-1].X()/nanosecond);

    for (TabulatedFunction::Iterator it = fAveragePulseShape->Begin();
         it != fAveragePulseShape->End(); ++it)
      fAnodePulseHisto->Fill(it->X()/nanosecond, it->Y());
  }

  // Read the voltage values for PMT saturation
  AttributeMap useAtt;
  useAtt["use"] = string("yes");

  Branch satB = elecB.GetChild("saturationCurve", useAtt);
  if (satB) {
    vector<double> vIn;
    satB.GetChild("inputVoltage").GetData(vIn);
    vector<double> vOut;
    satB.GetChild("outputVoltage").GetData(vOut);
    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; i < vIn.size(); ++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; i < vIn.size(); ++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");

  Branch limitStations = topB.GetChild("LimitStationsPerCycle", useAtt);
  if (limitStations) {
    limitStations.GetData(fNStationsPerCycle);
    fLimitStationsPerCycle = true;
  }

  fFirstCycle = true;
  return eSuccess;
}


void
SdPMTSimulator::ConvertPEToBaseSignal(const TimeDistributionI& pe, TimeDistributionD& base, double PMTGain)
{
  const double dt = pe.GetBinning();
  const double pulseStopTime = (*fAveragePulseShape)[fAveragePulseShape->GetNPoints()-1].X();
  double avgPulseShapeNorm = 0;
  vector<double> avgPulseShape;
  avgPulseShape.reserve(int(pulseStopTime / dt) + 1);
  for (double time = 0; time < pulseStopTime; time += dt) {
    const double pulse = fAveragePulseShape->Y(time);
    avgPulseShapeNorm += pulse;
    avgPulseShape.push_back(pulse);
  }

  const double pulseMaxTime = fMaxPulseHeight->X();
  const double scaleFactorMin = PMTGain * fMinCharge / avgPulseShapeNorm;
  const double scaleFactorMax = PMTGain * fMaxCharge / avgPulseShapeNorm;

  for (TimeDistributionI::SparseIterator peIt = pe.SparseBegin(), end = pe.SparseEnd();
       peIt != end; ++peIt) {

//#warning implement central-limit normal distribution for large npe
    double charge = 0;
    for (int i = 0, n = peIt->get<1>(); i < n; ++i)
      charge += fChargeDist->fire();

    const double realCharge = scaleFactorMin + scaleFactorMax*charge;

    const double offset = peIt->get<0>() - pulseMaxTime;

    // 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 * avgPulseShape[++i]);

  }
}


void
SdPMTSimulator::SimulateSaturation(TimeDistributionD& baseSignal,
                                   const TimeDistributionD& totalBaseSignal)
{
  for (int i = baseSignal.GetStart(), n = baseSignal.GetStop(); i <= n; ++i) {
    double& sig = baseSignal[i];
    const double tot = totalBaseSignal[i];
    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;
  }

  SEvent& sEvent = event.GetSEvent();

  if (fStationIterator == sEvent.StationsEnd()) {
    ostringstream info;
    info << "all stations processed.";
    INFO(info);
    fFirstCycle = true;
    fStationIterator = sEvent.StationsBegin();
    // this check is to be backward compatible with
    // module sequences that don't put SdPMTSimulator in a <loop>
    if (fLimitStationsPerCycle) {
      return eBreakLoop;
    }
  }

  if (fFirstCycle) {
    fStationIterator = sEvent.StationsBegin();
  }

  unsigned int fStationsProcessed = 0;

  for (; fStationIterator != sEvent.StationsEnd(); ++fStationIterator) {

    if (fStationsProcessed == fNStationsPerCycle) {
      fFirstCycle = false;
      return eSuccess;
    }

    for (Station::PMTIterator pIt = fStationIterator->PMTsBegin(sdet::PMTConstants::eAnyType);
         pIt != fStationIterator->PMTsEnd(sdet::PMTConstants::eAnyType); ++pIt) {

      if (!pIt->HasSimData())
        continue;

          const StationConstants::SignalComponent totalComp = StationConstants::eTotal;
          const unsigned int nFundComponents=((pIt->GetType() == sdet::PMTConstants::eWaterCherenkovLarge) ? fFundamentalComponents.size() : fFundamentalComponents_ASCII.size() );
          const double PMTGain=((pIt->GetType() == sdet::PMTConstants::eWaterCherenkovLarge) ? fPMTGain_WCD : fPMTGain_ASCII);

          PMTSimData& pSim = pIt->GetSimData();
          if (!pSim.HasPETimeDistribution(totalComp))
            continue; 
          
          if (pSim.HasBaseSignal(totalComp))
            continue;

          // no total PMT base signal exists yet
          pSim.MakeBaseSignal(totalComp);

          for (unsigned int i = 0; i != nFundComponents; ++i) {
            const StationConstants::SignalComponent sigComp = ((pIt->GetType() == sdet::PMTConstants::eWaterCherenkovLarge) ? fFundamentalComponents[i] : fFundamentalComponents_ASCII[i]);
            if (!pSim.HasPETimeDistribution(sigComp) ) continue;                
            if (!pSim.HasBaseSignal(sigComp)) pSim.MakeBaseSignal(sigComp);
            TimeDistributionD& baseSignal = pSim.GetBaseSignal(sigComp);
            ConvertPEToBaseSignal(pSim.GetPETimeDistribution(sigComp), baseSignal,PMTGain);
          }

          // Make total or add the fundamental components
          TimeDistributionD& totalBaseSignal = pSim.GetBaseSignal(totalComp);
          if ( nFundComponents == 0) {
            ConvertPEToBaseSignal(pSim.GetPETimeDistribution(totalComp), totalBaseSignal, PMTGain);
          }
          else {
            for (unsigned int i = 0; i != nFundComponents; ++i) {
              const StationConstants::SignalComponent sigComp = ((pIt->GetType() == sdet::PMTConstants::eWaterCherenkovLarge) ? fFundamentalComponents[i] : fFundamentalComponents_ASCII[i]);
              if (!pSim.HasPETimeDistribution(sigComp) ) continue;
              TimeDistributionD& componentSignal = pSim.GetBaseSignal(sigComp);
              for (TimeDistributionD::SparseIterator peIt = componentSignal.SparseBegin(), end = componentSignal.SparseEnd();
                   peIt != end; ++peIt) {
                totalBaseSignal[peIt->get<0>()] += peIt->get<1>();
              }
            }
          }
          if (fSaturationFunction)
            SimulateSaturation(totalBaseSignal, totalBaseSignal);

    }//Loop over PMTs
  
    ++fStationsProcessed;

  }//Loop over stations

  return eSuccess;
}


VModule::ResultFlag
SdPMTSimulator::Finish()
{
  if (fWantRootHistos) {
    TFile outFile(fRootHistoFilename.c_str(), "RECREATE");

    fChargeHisto->Write();
    fAnodePulseHisto->Write();

    //fPEHistos->Write();
    outFile.Write();
    outFile.Close();

    delete fChargeHisto;
    fChargeHisto = 0;
    delete fAnodePulseHisto;
    fAnodePulseHisto = 0;
  }

  // some of this root objects get owned by TFile
  delete fAveragePulseShape;
  fAveragePulseShape = 0;
  delete fChargeDist;
  fChargeDist = 0;
  delete fSaturationFunction;
  fSaturationFunction = 0;

  return eSuccess;
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: