CachedShowerRegeneratorOG/CachedShowerRegenerator.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "CachedShowerRegenerator.h"
#include "LogGaussSmearing.h"

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

#include <evt/Event.h>
#include <evt/ShowerSimData.h>

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

#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <sdet/Scintillator.h>

#include <mdet/MDetector.h>

#include <cdet/CDetector.h>
#include <cdet/Station.h>

#include <utl/ErrorLogger.h>
#include <utl/GeometryUtilities.h>
#include <utl/MathConstants.h>
#include <utl/Particle.h>
#include <utl/ParticleCases.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/RandomEngine.h>
#include <utl/TimeStamp.h>
#include <utl/Math.h>
#include <utl/Reader.h>
#include <utl/TabularStream.h>
#include <utl/TabulatedFunction.h>

#include <CLHEP/Random/RandFlat.h>
#include <CLHEP/Random/RandPoisson.h>

#include <cmath>
#include <sstream>
#include <vector>
#include <limits>

using namespace CachedShowerRegeneratorOG;
using namespace fwk;
using namespace evt;
using namespace sevt;
using namespace utl;
using namespace std;

using CLHEP::RandFlat;
using CLHEP::RandPoisson;


//#define DUMP_PARTICLES


namespace CachedShowerRegeneratorOG {

#ifdef DUMP_PARTICLES
  ShadowPtr<ofstream> gParticleDump;


  inline
  double
  RadiusXY(const Vector& v, const CoordinateSystemPtr& cs)
  {
    return sqrt(Sqr(v.GetX(cs)) + Sqr(v.GetY(cs)));
  }


#  define DUMP(_pStationDistance_, _pStationDistanceXY_, _timeShift_, _injectionZ_) \
    *gParticleDump \
      << sId << ' ' \
      << Point(0,0,0, dStationCS).GetX(showerCS) << ' ' \
      << Point(0,0,0, dStationCS).GetY(showerCS) << ' ' \
      << RadiusXY(corePosition - Point(0,0,0, dStationCS), showerCS) / meter << ' ' \
      << particle.GetType() << ' ' \
      << particle.GetKineticEnergy() / GeV << ' ' \
      << pWeight << ' ' \
      << avgN << ' ' \
      << n << ' ' \
      << pPosition.GetX(showerCS) << ' ' \
      << pPosition.GetY(showerCS) << ' ' \
      << sqrt(exp(pLnR2)) / meter << ' ' \
      << (_pStationDistance_) / meter << ' ' \
      << (_pStationDistanceXY_) / meter << ' ' \
      << (_timeShift_) / nanosecond << ' ' \
      << (_injectionZ_) / meter \
      << '\n'
#  define DUMP_REJECT DUMP(0, 0, 0, 0)

#else

#  define DUMP(_pStationDistance_, _pStationDistanceXY_, _timeShift_, _injectionZ_)
#  define DUMP_REJECT

#endif


  template<typename Map, typename T>
  inline
  void
  InsertValue(Map& map, const int sId, const T& value)
  {
    const auto sIt = map.find(sId);
    if (sIt != map.end())
      sIt->second(value);
    else
      map.insert(make_pair(sId, typename Map::mapped_type(value)));
  }


  inline
  double
  Round(const double div, const double val)
  {
    return round(div * val) / div;
  }



  inline
  double
  PlaneFrontTime(const CoordinateSystemPtr& showerCS, const Point& corePosition,
                 const Point& position)
  {
    // assuming t_c iz zero
    return (corePosition.GetZ(showerCS) - position.GetZ(showerCS)) / kSpeedOfLight;
  }


  inline
  bool
  IsMuonic(const Particle& p)
  {
    return
      p.GetType() == utl::Particle::eMuon ||
      p.GetType() == utl::Particle::eAntiMuon ||
      p.GetSource() == utl::Particle::eShowerFromMuonDecay ||
      p.GetSource() == utl::Particle::eShowerFromLocalHadron;
  }


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

    const unsigned int maxUI = numeric_limits<unsigned int>::max();

    topB.GetChild("LimitParticlesPerCycle").GetData(fLimitParticlesPerCycle);
    fParticlesPerCycle = maxUI;
    if (fLimitParticlesPerCycle)
      topB.GetChild("ParticlesPerCycle").GetData(fParticlesPerCycle);

    topB.GetChild("UseWeightedStationSimulation").GetData(fUseWeightedStationSimulation);
    fWeightedStationSimulationParticleLimit = maxUI;
    if (fUseWeightedStationSimulation) {
      topB.GetChild("WeightedStationSimulationParticleLimit").GetData(fWeightedStationSimulationParticleLimit);
      Branch wtFactorB = topB.GetChild("WeightedStationSimulationThinningFactor");
      if (wtFactorB)
        wtFactorB.GetData(fWeightedStationSimulationThinningFactor);
    }

    if (fParticlesPerCycle != maxUI && fWeightedStationSimulationParticleLimit != maxUI) {
      ERROR("The <LimitParticlesPerCycle> and <UseWeightedStationSimulation> options are incompatible, "
            "only one should be enabled.");
      return eFailure;
    }

    Branch distanceCutsB = topB.GetChild("DistanceCuts");
    fInnerRadiusCut = 0;
    const AttributeMap useAtt({{ "use", "yes" }});
    Branch ircB = distanceCutsB.GetChild("InnerRadiusCut", useAtt);
    if (ircB)
      ircB.GetData(fInnerRadiusCut);

    fOuterRadiusCut = 1e6*km;
    Branch orcB = distanceCutsB.GetChild("OuterRadiusCut", useAtt);
    if (orcB)
      orcB.GetData(fOuterRadiusCut);

    Branch energyCutsB = topB.GetChild("EnergyCuts");
    energyCutsB.GetChild("ElectronEnergyCut").GetData(fElectronEnergyCut);
    energyCutsB.GetChild("MuonEnergyCut").GetData(fMuonEnergyCut);
    energyCutsB.GetChild("PhotonEnergyCut").GetData(fPhotonEnergyCut);
    energyCutsB.GetChild("HadronEnergyCut").GetData(fHadronEnergyCut);
    energyCutsB.GetChild("MesonEnergyCut").GetData(fMesonEnergyCut);

    Branch algoB = topB.GetChild("AlgorithmParameters");
    algoB.GetChild("DeltaROverR").GetData(fDeltaROverR);
    algoB.GetChild("DeltaPhi").GetData(fDeltaPhi);
    algoB.GetChild("HorizontalParticleCosThetaCut").GetData(fHorizontalParticleCut);
    fHorizontalParticleCut = std::abs(fHorizontalParticleCut);
    algoB.GetChild("UseStationPositionMatrix").GetData(fUseStationPositionMatrix);
    algoB.GetChild("PhiGranularity").GetData(fPhiGranularity);
    algoB.GetChild("RGranularity").GetData(fRGranularity);
    algoB.GetChild("LogGaussSmearingWidth").GetData(fLogGaussSmearingWidth);
    algoB.GetChild("UseWeightDependentResamplingArea").GetData(fUseWeightDependentResamplingArea);

    fMuonWeightScale = 1;
    Branch muonWeightScaleB = topB.GetChild("MuonWeightScale");
    if (muonWeightScaleB)
      muonWeightScaleB.GetData(fMuonWeightScale);

    fSimulateUMD = false;
    Branch umdB = topB.GetChild("UMD");
    if (umdB) {
      umdB.GetChild("Simulate").GetData(fSimulateUMD);
      umdB.GetChild("TightRadius").GetData(fUMDTightRadius);
      umdB.GetChild("MaxRadius").GetData(fUMDMaxRadius);
      if (fSimulateUMD) {
        INFO("UMD tank+ground G4 simulation has been activated");
        if (fUseWeightedStationSimulation)
          WARNING("UMD does not handle weighted particles. The LimitParticlesPerCycle option should be used instead.");
      }
    }

    fSimulateMARTA = false;
    Branch martaB = topB.GetChild("MARTA");
    if (martaB) {
      martaB.GetChild("Simulate").GetData(fSimulateMARTA);
      martaB.GetChild("Radius").GetData(fMARTARadius);
      if (fSimulateMARTA)
        INFO("MARTA tank+RPC G4 simulation has been activated");
    }

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

    if (!fUseStationPositionMatrix)
      WARNING("Not using StationPositionMatrix. "
              "All stations will be searched for injection.");

  #ifdef DUMP_PARTICLES
    gParticleDump =
      new ofstream("CachedShowerRegenerator-particle_dump-"
                   "wdr" + to_string(int(fUseWeightDependentResamplingArea)) +
                   ".dat");
  #endif

    return eSuccess;
  }


  inline
  bool
  CachedShowerRegenerator::IsParticleEnergyLow(const int type,
                                               const double energy)
    const
  {
    switch (type) {
    case OFFLINE_PHOTON:
      return energy < fPhotonEnergyCut;
    case OFFLINE_ELECTRONS:
      return energy < fElectronEnergyCut;
    case OFFLINE_MUONS:
      return energy < fMuonEnergyCut;
    case OFFLINE_BARYONS:
      return energy < fHadronEnergyCut;
    case OFFLINE_MESONS:
      return energy < fMesonEnergyCut;
    default:
      return true;
    }
  }


  void
  CachedShowerRegenerator::InitNewShower(Event& event)
  {
    fShowerData = new ShowerData(StationPositionMatrix(fPhiGranularity, fRGranularity));

    fShowerData->fWeightCounterMap.clear();

    if (fLogGaussSmearingWidth)
      fShowerData->fLogGauss =
        new LogGaussSmearing(fLogGaussSmearingWidth, fRandomEngine);

    auto& simShower = event.GetSimShower();
    simShower.SetMuonWeightScale(fMuonWeightScale);

    if (simShower.GetMinRadiusCut() < fInnerRadiusCut) {
      ostringstream msg;
      msg << "Setting new value for fMinRadiusCut: " << fInnerRadiusCut;
      INFO(msg);
      simShower.SetMinRadiusCut(fInnerRadiusCut);
    }

    fShowerData->fParticleIt = simShower.GroundParticlesBegin();
    fShowerData->fParticlesEnd = simShower.GroundParticlesEnd();

    const auto localCS = simShower.GetLocalCoordinateSystem();
    const double showerZenith = (-simShower.GetDirection()).GetTheta(localCS);
    // area = pi (r1^2 - r2^2) (2 dPhi) / (2 pi). the 2 dPhi since its +- dPhi
    const double samplingAreaFactor = fDeltaPhi / cos(showerZenith);

    if (!event.HasSEvent())
      event.MakeSEvent();
    auto& sEvent = event.GetSEvent();

    const auto showerCS = simShower.GetShowerCoordinateSystem();
    const auto& corePosition = simShower.GetPosition();
    const auto& coreTime = simShower.GetTimeStamp();
    const auto& sDetector = det::Detector::GetInstance().GetSDetector();

    int nStations = 0;
    int nInner = 0;
    int nOuter = 0;

    for (const auto& sd : sDetector.StationsRange()) {

      ++nStations;

      const int sId = sd.GetId();

      if (!sEvent.HasStation(sId))
        sEvent.MakeStation(sId);

      auto& station = sEvent.GetStation(sId);

      // clear particle lists for the station. The stuff above
      // presumably only needs to be called once per shower, but
      // particle clearing has to happen each time through.
      if (!station.HasSimData())
        station.MakeSimData();
      auto& sSim = station.GetSimData();
      if (station.HasScintillator() && !station.GetScintillator().HasSimData())
        station.GetScintillator().MakeSimData();

      sSim.SetMaxNParticles(fWeightedStationSimulationParticleLimit);
      sSim.SetThinningFactor(fWeightedStationSimulationThinningFactor);

      const auto& sPosition = sd.GetPosition();

      const TimeStamp planeTime =
        coreTime + TimeInterval(PlaneFrontTime(showerCS, corePosition, sPosition));
      sSim.SetPlaneFrontTime(planeTime);

      const double sR = sPosition.GetRho(showerCS);

      if (sR < fInnerRadiusCut) {
        ++nInner;
        sSim.SetIsInsideMinRadius();
        continue;
      }

      if (sR > fOuterRadiusCut) {
        ++nOuter;
        continue;
      }

      // it is best to create station matrix based on log(r^2), see Reference Manual
      const double r1 = max(sR*(1 - fDeltaROverR), fInnerRadiusCut);
      const double r2 = min(sR*(1 + fDeltaROverR), fOuterRadiusCut);
      const double lnSqrR1 = 2 * log(r1/meter);
      const double lnSqrR2 = 2 * log(r2/meter);

      // DV: Particle weights are calculated by the shower simulation according
      // to the particle densities as measured in the ground-particle plane.
      // Thus all reweighting area ratios calculated in the resampler have to
      // have area values expressed in the corresponding ground-particle plane
      // projections, ie if there is a small missalignement of the tank vertical
      // and the ground-particle plane vertical, it has to be taken into account
      // with appropriate cosine of this small angle.
      // First calculate the ground-particle area of the resampling pie
      // section of a station:
      const double samplingArea = samplingAreaFactor * (Sqr(r2) - Sqr(r1));

      fShowerData->fStationMatrix.PushBack(sd, station,
                                           sPosition.GetPhi(showerCS), fDeltaPhi,
                                           sR, lnSqrR1, lnSqrR2,
                                           samplingArea);

    }

    fShowerData->fStationMatrix.CreateMatrix(fUseStationPositionMatrix);
    fShowerData->fMinSqrR = exp(fShowerData->fStationMatrix.GetMinLnSqrR())*meter2;

    ostringstream info;
    info << "Out of " << nStations << " stations: inner cut " << nInner << ", outer " << nOuter;
    INFO(info);
  }


  VModule::ResultFlag
  CachedShowerRegenerator::Run(Event& event)
  {
    if (!event.HasSimShower()) {
      ERROR("Current event does not have a simulated shower.");
      return eFailure;
    }

    const auto& simShower = event.GetSimShower();

    // sim event finished ?
    if (fShowerData && fShowerData->fParticleIt == fShowerData->fParticlesEnd) {
      fShowerData = nullptr;
      // break inner loop
      return eBreakLoop;
    }

    if (!fShowerData)
      InitNewShower(event);

    const unsigned int maxUI = numeric_limits<unsigned int>::max();

    if (fParticlesPerCycle < maxUI) {
      ostringstream info;
      info << "Regenerating batch of " << fParticlesPerCycle << " particles.";
      INFO(info);
    }

    if (fWeightedStationSimulationParticleLimit < maxUI) {
      ostringstream info;
      info << "Will use weighted simulation for more than " << fWeightedStationSimulationParticleLimit
           << " particles per station.";
      INFO(info);
    }

    const auto& sDetector = det::Detector::GetInstance().GetSDetector();
    auto& sEvent = event.GetSEvent();

    // clear particle list for the stations
    for (const auto& sd : sDetector.StationsRange())
      sEvent.GetStation(sd.GetId()).GetSimData().ClearParticleList();

    const auto showerCS = simShower.GetShowerCoordinateSystem();
    const auto& corePosition = simShower.GetPosition();
    const auto grParticleCS = simShower.GetLocalCoordinateSystem();

    // Based on the cuts read in for the outer and inner radius in Init(), we
    // build the list of stations that lie within the region of interest in
    // the shower plane. We then loop only over these stations, since anything
    // outside the limits defined by the outer and inner radius cuts will never
    // receive any regenerated particles anyway.

    // keep count to avoid excessive memory consumption
    unsigned int nParticlesInjected = 0;

    auto& particleIt = fShowerData->fParticleIt;
    for (auto pEnd = fShowerData->fParticlesEnd; particleIt != pEnd; ++particleIt) {

      // over the cache limit?
      if (fParticlesPerCycle && nParticlesInjected >= fParticlesPerCycle) {
        ostringstream info;
        info << nParticlesInjected << " particles processed.";
        INFO(info);
        return eSuccess;
      }

      const auto& particle = *particleIt;

      const auto& pPosition = particle.GetPosition();
      const double pR2ShowerFrame = pPosition.GetRho2(showerCS);

      if (pR2ShowerFrame <= fShowerData->fMinSqrR)
        continue;

      if (IsParticleEnergyLow(particle.GetType(), particle.GetKineticEnergy()))
        continue;

      const double pPhi = pPosition.GetPhi(showerCS);
      const double pLnSqrR = log(pR2ShowerFrame);

      const auto& sList =
        fShowerData->fStationMatrix.GetInjectionStations(pPhi, pLnSqrR);

      if (sList.empty())
        continue;

      // particle time relative to the core time
      const double pTime = particle.GetTime().GetInterval();

      // due to a bug in CORSIKA (version 6.900) when simulating upward going particles (option UPWARD)
      // a few particles per shower get negative delays w.r.t. to the plane front, they behave like tachyons.
      // see http://www-ik.fzk.de/~maurel/corsika_negative_particle_delays for a list
      // of "superluminal" particles in a library simulated using the UPWARD option.
      // Negative delays on the order of 0.01 ns are possible due to roundoff errors,
      // only particles with a delay smaller than -0.1 ns are rejected here.
      const double pPlaneFrontDelay = pTime - PlaneFrontTime(showerCS, corePosition, pPosition);
      if (pPlaneFrontDelay < -0.1*ns)
        continue;

      Particle newParticle = particle;
      newParticle.SetWeight(1);

      for (const auto sInfo : sList) {

        const auto& dStation = sInfo->GetDetStation();
        const auto stationCS = dStation.GetLocalCoordinateSystem();
        const unsigned int sId = dStation.GetId();
        const auto& pDirection = newParticle.GetDirection();

        // positive for down-going particles
        const double pStationCosine = -pDirection.GetZ(stationCS);
        // skip upward-going particles
        if (pStationCosine <= 0) {
          ++fShowerData->fUpwardSideParticles[sId];
          continue;
        }

        // positive for down-going particles
        const double pPlaneCosine = -pDirection.GetZ(grParticleCS);
        // skip horizontal particles
        if (pPlaneCosine < fHorizontalParticleCut) {
          ++fShowerData->fNHorizontalParticles;
          continue;
        }

        const double sThickness = dStation.GetThickness();
        const double sRadius = dStation.GetRadius() + sThickness;
        const double sHeight = dStation.GetHeight() + 2*sThickness;

        double injectionRadius = sRadius;
        double injectionHeight = sHeight;

        // if MARTA simulation is enabled the tank is raised by the precast total height and so should the injection height
        if (fSimulateMARTA) {
          const auto& cStation =
            det::Detector::GetInstance().GetCDetector().GetStation(sId);
          injectionHeight +=
            cStation.GetTankSupportTopSlabDimensions()[2] +
            cStation.GetTankSupportCentralFootDimensions()[2] +
            cStation.GetTankSupportCentralFootBaseDimensions()[2];
          injectionRadius = max(injectionRadius, fMARTARadius);
        }

        if (dStation.HasScintillator()) {
          // Check if scintillator requires larger injection cylinder
          const double scintRadius = dStation.GetScintillator().GetMaxRadius();
          const double scintHeight = dStation.GetScintillator().GetMaxHeight();
          injectionRadius = max(injectionRadius, scintRadius);
          injectionHeight = max(injectionHeight, scintHeight);
        }

        if (fSimulateUMD && dStation.ExistsAssociatedCounter()) {
          const auto& mDet = det::Detector::GetInstance().GetMDetector();
          const auto& dCounter = mDet.GetCounter(dStation.GetAssociatedCounterId());
          const auto& mod = *dCounter.ModulesBegin();
          const auto& modPos = mod.GetPosition();
          const double depth = -modPos.GetZ(stationCS);
          const double pPlaneSine = pDirection.GetRho(grParticleCS);
          // increase radius for inclined particles so that shadow always hits UMD + depth * tan(theta)
          const double umdRadius = fUMDTightRadius + depth * pPlaneSine/pPlaneCosine;
          injectionRadius = min(fUMDMaxRadius, max(injectionRadius, umdRadius));
        }

        const double pStationSine = pDirection.GetRho(stationCS);

        // "eff" here stands for effective areas as seen from the particle standpoint
        const double effAreaTop = kPi * Sqr(injectionRadius) * pStationCosine;
        const double effAreaSide = 2 * injectionRadius * injectionHeight * pStationSine;
        const double effArea = effAreaTop + effAreaSide;  // total

        const double pWeight = IsMuonic(particle) ?
          fMuonWeightScale * particle.GetWeight() : particle.GetWeight();

        const double effResArea = sInfo->GetResamplingArea() * pPlaneCosine;

        // all particle positions within the resampling area are equivalent, define
        // "weight density", ie the probability to find a particle with such properties
        // within some area
        const double avgN = pWeight * effArea / effResArea;

        unsigned int n = 0;
        if (fUseWeightDependentResamplingArea && avgN < 1) {
          // see GAP-2015-086 for details
          if (!sInfo->IsInScaledArea(avgN, pPhi, pLnSqrR)) {
            DUMP_REJECT;
            continue;
          }
          // direct injection
          InsertValue(fShowerData->fWeightStat, sId, 1);
          n = 1;
        } else {
          // cloning
          InsertValue(fShowerData->fWeightStat, sId, avgN);
          n = RandPoisson::shoot(fRandomEngine, avgN);
        }

        auto& station = sInfo->GetEvtStation();

        for (unsigned int i = 0; i < n; ++i) {

          if (RandFlat::shoot(fRandomEngine, 0, effArea) < effAreaTop) {
            // top injection
            const double r = injectionRadius * sqrt(RandFlat::shoot(fRandomEngine, 0, 1));
            const double phi = RandFlat::shoot(fRandomEngine, 0, kTwoPi);
            newParticle.SetPosition(Point(r, phi, injectionHeight, stationCS, Point::kCylindrical));
            ++fShowerData->fTopParticles[sId];
          } else {
            // side injection
            const double pPhi = pDirection.GetPhi(stationCS) + kPiOnTwo;
            // For the X and Y component we need to find the angle made between
            // the injection cylinder's X-Y axis and the incoming particle direction,
            // then add an angle sampled from a cosine distribution.
            const double phi = pPhi + acos(RandFlat::shoot(fRandomEngine, -1, 1));
            const double z = RandFlat::shoot(fRandomEngine, 0, injectionHeight);
            newParticle.SetPosition(Point(injectionRadius, phi, z, stationCS, Point::kCylindrical));
            ++fShowerData->fSideParticles[sId];
          }

          // Shift the time to account for differences in position between the
          // hit tank and the position the particle was placed at the ground
          // by the shower simulation. Also, include a shift for multiple
          // particle entries from regeneration of one weighted particle. So far,
          // no one has provided any real justification for what distribution
          // should be used for the energy shifting.
          // see P. Billoir, Astroparticle Physics 30 (2008) 270-285
          const auto& pShiftedPosition = newParticle.GetPosition();
          double pShiftedTime = pTime + PlaneFrontTime(showerCS, pPosition, pShiftedPosition);
          if (i && fShowerData->fLogGauss) {
            const double planeFrontTime = PlaneFrontTime(showerCS, corePosition, pShiftedPosition);
            pShiftedTime = fShowerData->fLogGauss->GetSmearedTime(planeFrontTime, pShiftedTime);
          }
          newParticle.SetTime(pShiftedTime);

          InsertValue(fShowerData->fTimeStat, sId, pShiftedTime);
          station.AddParticle(newParticle);

          DUMP(
            (pPosition - Point(0,0,0, dStationCS)).GetMag(),
            RadiusXY(pPosition - Point(0,0,0, dStationCS), showerCS),
            pShiftedTime - pTime,
            newParticle.GetPosition().GetZ(dStationCS)
          );

        }

        nParticlesInjected += n;

      }

    }

    ostringstream info;
    info << nParticlesInjected << " particles processed.";
    INFO(info);

    if (fShowerData->fParticleIt == fShowerData->fParticlesEnd)
      OutputStats(event);

  #ifdef DUMP_PARTICLES
    *gParticleDump << "\n\n\n---\n\n\n" << endl;
  #endif

    return eSuccess;
  }


  void
  CachedShowerRegenerator::OutputStats(Event& event)
  {
    const auto& sList = fShowerData->fStationMatrix.GetStationList();

    int nStations = 0;

    {
      INFO("\nStation particle statistics:");

      TabularStream tab("r|.|r|r|r|.|.|.|l");
      tab              << endc << "Rho"  << endc           << endc              << endc         << endc
          << "wght"    << endc << "wght" << endc << "wght" << endc << "station"
          << endr
          << "station" << endc << "[km]" << endc << "top"  << endc << "side"    << endc << "up" << endc
          << "min"     << endc << "avg"  << endc << "max"  << endc << "thinning"
          << endr << hline;

      for (const auto& s : sList) {

        const unsigned sId = s.GetDetStation().GetId();

        if (fShowerData->fTopParticles[sId] ||
            fShowerData->fSideParticles[sId] ||
            fShowerData->fUpwardSideParticles[sId]) {

          const auto& w = *fShowerData->fWeightStat.find(sId);
          tab << sId << ' ' << endc
              << setprecision(4) << s.GetR()/km << endc
              << ' ' << fShowerData->fTopParticles[sId] << ' ' << endc
              << ' ' << fShowerData->fSideParticles[sId] << ' ' << endc
              << ' ' << fShowerData->fUpwardSideParticles[sId] << ' ' << endc
              << ' ' << setprecision(4) << w.second.GetMin() << endc
              << ' ' << setprecision(4) << w.second.GetAverage() << endc
              << ' ' << setprecision(4) << w.second.GetMax() << endc
              << s.GetEvtStation().GetSimData().GetThinning()
              << endr;
          ++nStations;
        }

      }

      tab << delr;

      if (nStations)
        cerr << tab << endl;
      else
        cerr << "No stations were hit." << endl;
    }

    {
      INFO("\nStation timing statistics:\n"
           "abs = absolute time difference to core time\n"
           "rel = relative time to station plane-front arrival");
      TabularStream tab("r|.|.|.|.|.");
      tab <<              endc << "Rho"  << endc << "min abs"   << endc << "max abs"   << endc
          << "min rel"   << endc << "max rel"   << endr

          << "station" << endc << "[km]" << endc << "time [ns]" << endc << "time [ns]" << endc
          << "time [ns]" << endc << "time [ns]" << endr << hline;

      const auto& simShower = event.GetSimShower();
      const auto& coreTime = simShower.GetTimeStamp();

      for (const auto& s : sList) {

        const unsigned sId = s.GetDetStation().GetId();

        if (fShowerData->fTopParticles[sId] ||
            fShowerData->fSideParticles[sId] ||
            fShowerData->fUpwardSideParticles[sId]) {

          const TimeInterval sTime =
            s.GetEvtStation().GetSimData().GetPlaneFrontTime() - coreTime;

          const auto& mm = *fShowerData->fTimeStat.find(sId);
          tab << sId << ' ' << endc
              << setprecision(4) << s.GetR()/km << endc
              << ' ' << Round(10, mm.second.GetMin().GetInterval()/ns) << endc
              << ' ' << Round(10, mm.second.GetMax().GetInterval()/ns) << endc
              << ' ' << Round(10, (mm.second.GetMin() - sTime).GetInterval()/ns) << endc
              << ' ' << Round(10, (mm.second.GetMax() - sTime).GetInterval()/ns) << endr;

        }

      }

      tab << delr;

      if (nStations)
        cerr << tab << endl;
    }

    if (nStations) {
      ostringstream info;
      info << fShowerData->fNHorizontalParticles << " horizontal particles rejected.";
      INFO(info);
    }
  }

}