Deprecated/UpgradeASCIITests/SdBaselineSimulatorASCII/SdBaselineSimulator.cc

Go to the documentation of this file.

// SdBaselineSimulator - take fadc trace and add pedestal
// 25 Feb 2004 T. Ohnuki UCLA
// revised Jun 2004  T. Paul

#include <sstream>

#include "SdBaselineSimulator.h"

#include <evt/Event.h>

#include <det/Detector.h>

#include <sdet/SDetector.h>
#include <sdet/Station.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 <CLHEP/Random/RandGauss.h>

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

using CLHEP::RandGauss;


SdBaselineSimulator::SdBaselineSimulator() :
  fPrefixBins(0),
  fPostfixBins(0),
  fRandomEngine(0),
  fStoreNonSaturated(0)
{
}


VModule::ResultFlag
SdBaselineSimulator::Init()
{
  Branch topBranch =
    CentralConfig::GetInstance()->GetTopBranch("SdBaselineSimulatorASCII");

  topBranch.GetChild("NumberPrefixBins").GetData(fPrefixBins);
  topBranch.GetChild("NumberPostfixBins").GetData(fPostfixBins);

  topBranch.GetChild("StoreNonSaturated").GetData(satString);

  fStoreNonSaturated = (satString == "yes");

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

  return eSuccess;
}


VModule::ResultFlag
SdBaselineSimulator::Run(evt::Event& event)
{
  SEvent& sEvent = event.GetSEvent();

  for (SEvent::StationIterator sIt = sEvent.StationsBegin();
       sIt != sEvent.StationsEnd(); ++sIt) {

    // skip station if no particles
    if (!sIt->HasSimData())
      continue;

    const sdet::Station& dStation = det::Detector::GetInstance().GetSDetector().GetStation(*sIt);

    string simModuleName = sIt->GetSimData().GetSimulatorSignature();


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

      const sdet::PMT& dPMT = dStation.GetPMT(*pmtIt);
      const unsigned int saturationValue = dStation.GetSaturationValue();

      // Only add baseline to the total trace
            if (pmtIt->HasSimData()) {
              
        PMTSimData& pmtSim = pmtIt->GetSimData();

        if (pmtSim.HasFADCTrace()) {

          sevt::StationConstants::SignalComponent component=(StationConstants::eTotal); 
          sevt::StationConstants::SignalComponent componentNoSat=(StationConstants::eTotalNoSaturation);

          double pedHG = dPMT.GetBaseline(sdet::PMTConstants::eHighGain);
          double noiseHG = dPMT.GetBaselineRMS(sdet::PMTConstants::eHighGain);
          double pedLG = dPMT.GetBaseline(sdet::PMTConstants::eLowGain);
          double noiseLG = dPMT.GetBaselineRMS(sdet::PMTConstants::eLowGain);

          TimeDistributionI& hgTrace = pmtSim.GetFADCTrace(PMTConstants::eHighGain, component);
          TimeDistributionI& lgTrace = pmtSim.GetFADCTrace(PMTConstants::eLowGain,  component);
          TimeDistributionD& hgTrace_d = pmtSim.GetFADCTraceD(PMTConstants::eHighGain, component);
          TimeDistributionD& lgTrace_d = pmtSim.GetFADCTraceD(PMTConstants::eLowGain, component);

          if (fStoreNonSaturated) {
            pmtSim.MakeFADCTrace(componentNoSat);
            TimeDistributionD& hgUnSatTrace_d = pmtSim.GetFADCTraceD(PMTConstants::eHighGain, componentNoSat);
            TimeDistributionD& lgUnSatTrace_d = pmtSim.GetFADCTraceD(PMTConstants::eLowGain, componentNoSat);
            AddPedestal(hgTrace, hgTrace_d, hgUnSatTrace_d, saturationValue, pedHG, noiseHG);
            AddPedestal(lgTrace, lgTrace_d, lgUnSatTrace_d, saturationValue, pedLG, noiseLG);
          } else {
            TimeDistributionD dummy(hgTrace.GetBinning());
            AddPedestal(hgTrace, hgTrace_d, dummy, saturationValue, pedHG, noiseHG);  // don't bother writing unsaturated to event
            AddPedestal(lgTrace, lgTrace_d, dummy, saturationValue, pedLG, noiseLG);
        }
            }// Has FADC trace    
          }// Has Sim data
        }// PMT loop
    
  }

  return eSuccess;
}


VModule::ResultFlag
SdBaselineSimulator::Finish()
{
  return eSuccess;
}


int
SdBaselineSimulator::AddPedestal(TimeDistributionI& trace,
                                 TimeDistributionD& doubleTrace,
                                 TimeDistributionD& unsaturatedTrace,
                                 const unsigned int saturationValue,
                                 const double ped, const double noise)
{
  const int firstBin = trace.GetStart() - fPrefixBins;
  const int lastBin = trace.GetStop() + fPostfixBins;
  RandomEngine::RandomEngineType* const randomEngine =
    &fRandomEngine->GetEngine();

  for (int i = firstBin; i <= lastBin; ++i) {
    const double diff = RandGauss::shoot(randomEngine, ped, noise);
    const int idiff = int(round(diff));
    
    unsaturatedTrace.SetTime(i, max(doubleTrace.At(i) + idiff, 0.));
    doubleTrace.SetTime(i, min(unsaturatedTrace.At(i), double(saturationValue)));

    const int value = max(trace.At(i) + idiff, 0);
    const int realValue = min(value, int(saturationValue));
    trace.SetTime(i, realValue);
  }

  return 0;
}