TabulatedTankSimulator.cc

Go to the documentation of this file.

#include "TabulatedTankSimulator.h"
#include <utl/ErrorLogger.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerUnits.h>
#include <utl/ParticleCases.h>

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

#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>
#include <det/Detector.h>
#include <sdet/SDetector.h>

#include <CLHEP/Random/RandFlat.h>


//using namespace sdet;
using namespace std;
using namespace sevt;
using namespace evt;
using namespace fwk;
using namespace utl;
using CLHEP::RandFlat;
using namespace TabulatedTankSimulatorNS;


namespace TabulatedTankSimulatorNS {
  const double pulse_Tau = 60.; //ns
  const double pulse_tMax = 25.; //ns
  const double pulse_cdfMax = 0.3; //Probability (time <= pulse_tMax)
}


TabulatedTankSimulator::TabulatedTankSimulator() :
  fTauMuon(kMuonLifetime/ns)
{ }


TabulatedTankSimulator::~TabulatedTankSimulator()
{
  delete fTankResponse;
}


VModule::ResultFlag
TabulatedTankSimulator::Init()
{
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("TabulatedTankSimulator");

  string sepMode;
  topB.GetChild("signalSeparationMode").GetData(sepMode);
  if (sepMode == "Standard")
    fSignalSeparationMode = eStandard;
  else if (sepMode == "Universality")
    fSignalSeparationMode = eUniversality;

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

  delete fTankResponse;

  bool useG4 = false;
  // 0: WCD, 1: ASCII, 2: AMIGA (this module can be used only for WCD)
  int DetectorType = 0;
  bool UseEnlargeArea = false;
  fTankResponse = new DetectorResponse(useG4, DetectorType, UseEnlargeArea);

  return eSuccess;
}


VModule::ResultFlag
TabulatedTankSimulator::Run(evt::Event& event)
{
  INFO("TabulatedTankSimulator::Run()");

  SEvent& sEvent = event.GetSEvent();
  for (SEvent::StationIterator sIt = sEvent.StationsBegin();
       sIt != sEvent.StationsEnd(); ++sIt)
    SimulateStation(*sIt);

  return eSuccess;
}


VModule::ResultFlag
TabulatedTankSimulator::Finish()
{
  INFO("TabulatedTankSimulator::Finish()");
  return eSuccess;
}


int
TabulatedTankSimulator::PDGToCDF(const int pid)
{
  switch (pid) {
  case utl::Particle::eMuon: return 0;
  case utl::Particle::eAntiMuon: return 0;
  case utl::Particle::eElectron: return 1;
  case utl::Particle::ePositron: return 1;
  case utl::Particle::ePhoton: return 2;
  case utl::Particle::eProton: return -5000;
  case utl::Particle::eAntiProton: return -5000;
  case utl::Particle::eNeutron: return -5000;
  case utl::Particle::eAntiNeutron: return -5000;
  default: return -5000;
  }
}


double
TabulatedTankSimulator::GetTimePE()
{
  double time = 0;
  double rnd = RandFlat::shoot(fRandomEngine, 0., 1.);

  if (rnd < pulse_cdfMax)
    time = rnd * pulse_tMax / pulse_cdfMax; //Simple linear interpolation
  else {
    rnd = (rnd - pulse_cdfMax) / (1 - pulse_cdfMax);
    time = pulse_tMax - pulse_Tau * log(1 - rnd);
  }

  if (!std::isnan(time) && !std::isinf(time))
    return time;

  return 25.;
}


void
TabulatedTankSimulator::SimulateStation(sevt::Station& station)
{
  if (!station.HasSimData())
    return;

  StationSimData& simData = station.GetSimData();
  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

  simData.SetSimulatorSignature("TabulatedTankSimulatorKG");

  const unsigned long numParticles =
    simData.ParticlesEnd() - simData.ParticlesBegin();

  if (!numParticles)
    return;

  /*
   * Counting the number of particles that are actually simulated (post
   * station-level thinning). Due to resampling simulation option that
   * uses cycles with maximum particle number per cycle, we always check
   * if particles have already been simulated.
   */
  const unsigned int nParticlesAlreadySimulated = simData.GetTotalSimParticleCount();
  if (!nParticlesAlreadySimulated)
    simData.SetTotalSimParticleCount(numParticles);
  else
    simData.SetTotalSimParticleCount(nParticlesAlreadySimulated + numParticles);

  const sdet::Station& detStation = sDetector.GetStation(station);
  for (sevt::Station::PMTIterator pmtIt = station.PMTsBegin();
       pmtIt != station.PMTsEnd(); ++pmtIt) {
    if (!pmtIt->HasSimData()) {
      pmtIt->MakeSimData();
      pmtIt->GetSimData().MakePETimeDistribution();
    }
  }

  for (StationSimData::ParticleIterator pIt = simData.ParticlesBegin();
       pIt != simData.ParticlesEnd(); ++pIt) {

    const double pke = pIt->GetKineticEnergy();
    const double ptheta = kPi - pIt->GetDirection().GetTheta(detStation.GetLocalCoordinateSystem());

    const int itype = PDGToCDF(pIt->GetType());
    if (itype == -5000)
      continue;

    const double rnd = RandFlat::shoot(fRandomEngine, 0, 1);
    const int nPe = int(fTankResponse->GetPe(pke/GeV, ptheta, itype, rnd));

    int nPe_decay = 0;

    if (itype == 0) {
      const double rnd_d = RandFlat::shoot(fRandomEngine, 0, 1);
      nPe_decay = int(fTankResponse->GetPe(pke/GeV, ptheta, 4, rnd_d));  // type 4: decay electrons
    }

    if (nPe + nPe_decay == 0)
      continue;

    ostringstream msg;
    msg << "Theta " << ptheta*180/3.14 << "|   KE  GeV " << pke/GeV << " |    Id " << itype
        << " | NPE=  " << nPe << " - "  << nPe_decay;
    //INFO(msg);

    double tdecay = 0;
    if (nPe_decay > 0)
      tdecay = -fTauMuon * log(1 - RandFlat::shoot(fRandomEngine, 0, 1));

    for (int kk = 0; kk < (nPe + nPe_decay); ++kk) {

      TimeInterval pe_time(pIt->GetTime());
      pe_time += GetTimePE();
      if (kk >= nPe)
        pe_time += tdecay;

      const int pmtId = int(RandFlat::shoot(fRandomEngine, 0, 3)) + 1;
      if (!station.HasPMT(pmtId))
        continue;

      if (!station.GetPMT(pmtId).HasSimData())
        station.GetPMT(pmtId).MakeSimData();
      PMTSimData& pmtSim = station.GetPMT(pmtId).GetSimData();
      AddPhoton(*pIt, pmtSim, pe_time.GetInterval());
    }
  }
}


void
TabulatedTankSimulator::AddPhoton(const utl::Particle& particle, sevt::PMTSimData& pmtSim, const double peTime)
  const
{
  const int type = particle.GetType();
  const int source = particle.GetSource();

  sevt::StationConstants::SignalComponent component = sevt::StationConstants::eTotal;
  sevt::StationConstants::SignalComponent extraComponent = sevt::StationConstants::eTotal;

  switch (fSignalSeparationMode) {
  case eStandard:
    switch (type) {
    case OFFLINE_ELECTRONS:
      component = sevt::StationConstants::eElectron;
      if (source == utl::Particle::eShowerFromMuonDecay)
        extraComponent = sevt::StationConstants::eShowerMuonDecayElectron;
      break;
    case OFFLINE_PHOTON:
      component = sevt::StationConstants::ePhoton;
      if (source == utl::Particle::eShowerFromMuonDecay)
        extraComponent = sevt::StationConstants::eShowerMuonDecayPhoton;
      break;
    case OFFLINE_MUONS:
      component = sevt::StationConstants::eMuon;
      break;
    case OFFLINE_HADRONS:
      component = sevt::StationConstants::eHadron;
      break;
    default:
      cerr << "Unknown particle type" << endl;
      exit(-1);
    }
    break;
  case eUniversality:
    switch (type) {
    case OFFLINE_ELECTRONS:
      switch (source) {
      case utl::Particle::eShowerFromLocalHadron:
        component = sevt::StationConstants::eShowerLocalHadronElectron;
        break;
      case utl::Particle::eShowerFromMuonDecay:
        component = sevt::StationConstants::eShowerMuonDecayElectron;
        break;
      default:
        component = sevt::StationConstants::eElectron;
        break;
      }
      break;
    case OFFLINE_PHOTON:
      switch (source) {
      case utl::Particle::eShowerFromLocalHadron:
        component = sevt::StationConstants::eShowerLocalHadronPhoton;
        break;
      case utl::Particle::eShowerFromMuonDecay:
        component = sevt::StationConstants::eShowerMuonDecayPhoton;
        break;
      default:
        component = sevt::StationConstants::ePhoton;
        break;
      }
      break;
    case OFFLINE_MUONS:
      component = sevt::StationConstants::eMuon;
      break;
    case OFFLINE_HADRONS:
      component = sevt::StationConstants::eHadron;
      break;
    default:
      cerr << "Unknown particle type" << endl;
      exit(-1);
    }
    break;
  }

  if (!pmtSim.HasPETimeDistribution())
    pmtSim.MakePETimeDistribution();
  pmtSim.GetPETimeDistribution().AddTime(peTime);

  if (!pmtSim.HasPETimeDistribution(component))
    pmtSim.MakePETimeDistribution(component);
  pmtSim.GetPETimeDistribution(component).AddTime(peTime);

  if (extraComponent != sevt::StationConstants::eTotal) {
    if (!pmtSim.HasPETimeDistribution(extraComponent))
      pmtSim.MakePETimeDistribution(extraComponent);
    pmtSim.GetPETimeDistribution(extraComponent).AddTime(peTime);
  }

}