CachedXShowerRegenerator.cc

Go to the documentation of this file.

#include <utl/config.h>

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

#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 <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 <CLHEP/Random/RandGauss.h>

#include "CachedXShowerRegenerator.h"
#include "LogXGaussSmearing.h"

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

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

namespace CachedXShowerRegeneratorAG {

  template<typename Map, typename T>
  inline
  void
  InsertValue(Map& map, const int sId, const T& value)
  {
    const typename Map::iterator 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)
  {
    return (corePosition.GetZ(showerCS) - position.GetZ(showerCS)) / kSpeedOfLight;
  }

} // namespace CachedXShowerRegeneratorAG


CachedXShowerRegenerator::CachedXShowerRegenerator() :
  fMaxParticles(numeric_limits<unsigned int>::max()),
  fInnerRadiusCut(0.),
  fOuterRadiusCut(1.e6*km),
  fElectronEnergyCut(0),
  fMuonEnergyCut(0),
  fPhotonEnergyCut(0),
  fHadronEnergyCut(0),
  fMesonEnergyCut(0),
  fDeltaROverR(0),
  fDeltaPhi(0),
  fHorizontalParticleCut(0),
  fUseStationPositionMatrix(true),
  fPhiGranularity(0),
  fRGranularity(0),
  fMuToVem(0),
  fEToVem(0),
  fUseWeightLimiting(false),
  fWeightThreshold(0),
  fAccumulatedWeightLimit(0),
  fMuonWeightScale(1),
  fRandomEngine(0),
  fTimeSlotSize(1.*ns),
  fNOutOfRangeWeight(0),
  fWeightLimit(0),
  fSimulateAmiga(false)
{
}


VModule::ResultFlag
CachedXShowerRegenerator::Init()
{

  unityWeightCtr = 0; // temp
  largeWeightCtr = 0; // temp

  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("CachedXShowerRegenerator");

  AttributeMap useAtt;
  useAtt["use"] = string("yes");

  Branch limitB = topB.GetChild("LimitParticlesPerCycle", useAtt);
  if (limitB)
    limitB.GetData(fMaxParticles);

  Branch distanceCutB = topB.GetChild("DistanceCuts");

  Branch inB = distanceCutB.GetChild("InnerRadiusCut", useAtt);
  if (inB)
    inB.GetData(fInnerRadiusCut);

  Branch outB = distanceCutB.GetChild("OuterRadiusCut", useAtt);
  if (outB)
    outB.GetData(fOuterRadiusCut);

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

  Branch algoB = topB.GetChild("AlgorithmParameters");
  algoB.GetChild("DeltaROverR").GetData(fDeltaROverR);
  algoB.GetChild("DeltaPhi").GetData(fDeltaPhi);
  algoB.GetChild("HorizontalParticleCut").GetData(fHorizontalParticleCut);
  algoB.GetChild("UseStationPositionMatrix").GetData(fUseStationPositionMatrix);
  algoB.GetChild("PhiGranularity").GetData(fPhiGranularity);
  algoB.GetChild("RGranularity").GetData(fRGranularity);
  algoB.GetChild("LogXGaussSmearingWidth").GetData(fLogXGaussSmearingWidth);

  Branch limitSwitchB = topB.GetChild("ResamplingWeightLimiting", useAtt);
  if (limitSwitchB) {
    fUseWeightLimiting = true;
    limitSwitchB.GetChild("WeightThreshold").GetData(fWeightThreshold);
    limitSwitchB.GetChild("AccumulatedWeightLimit").GetData(fAccumulatedWeightLimit);
  }

  if (topB.GetChild("MuonWeightScale"))
    topB.GetChild("MuonWeightScale").GetData(fMuonWeightScale);

  Branch amigaSimB = topB.GetChild("AmigaSimulation", useAtt);
  if ( amigaSimB ) {
    INFO( "AMIGA tank+ground G4 simulation has been activated" );
    fSimulateAmiga = true;
    amigaSimB.GetChild("AmigaRadius").GetData(fAmigaRadius);
  }

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

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

  return eSuccess;
}


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


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

  fShowerData->fWeightCounterMap.clear();

  if (fLogXGaussSmearingWidth)
    fShowerData->fLogXGauss =
      new LogXGaussSmearing(fLogXGaussSmearingWidth, fRandomEngine);

  ShowerSimData& 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 CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
  const double showerZenith = (-simShower.GetDirection()).GetTheta(localCS); 
  const double samplingAreaFactor = 4 * fDeltaPhi * fDeltaROverR / cos(showerZenith);

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

  const CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  const Point& corePosition = simShower.GetPosition();
  const TimeStamp& coreTime = simShower.GetTimeStamp();
  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

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

  for (sdet::SDetector::StationIterator sdIt = sDetector.StationsBegin();
       sdIt != sDetector.StationsEnd(); ++sdIt) {

    ++nStations;

    const int sId = sdIt->GetId();

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

    Station& 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();

    StationSimData& sSim = station.GetSimData();

    const Point& sPosition = sdIt->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;
    }

    // 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(sR);

    // it is best to create station matrix based on log(r^2)
    const double r1 = 2 * log(sR * (1 - fDeltaROverR));
    const double r2 = 2 * log(sR * (1 + fDeltaROverR));

    //cout << "Pushing back station " << sdIt->GetId() << endl;
    fShowerData->fStationMatrix.PushBack(*sdIt, station,
                                         sPosition.GetPhi(showerCS), fDeltaPhi,
                                         sR, r1, r2,
                                         samplingArea);

  }

  fShowerData->fStationMatrix.CreateMatrix(fUseStationPositionMatrix);
  fShowerData->fMinR2 = exp(fShowerData->fStationMatrix.GetMinR());

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


VModule::ResultFlag
CachedXShowerRegenerator::Run(Event& event)
{
  INFO("Starting G4XTank simulation...");

  if (!event.HasSimShower()) {
    ERROR("Current event does not have a simulated shower.");
    return eFailure;
  }

  const ShowerSimData& simShower = event.GetSimShower();

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

  if (!fShowerData)
    InitNewShower(event);

  {
    ostringstream info;
    info << "Regenerating batch of "
         << fMaxParticles << " particles.";
    INFO(info);
  }

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

  // clear particle list for the stations
  for (sdet::SDetector::StationIterator sdIt = sDetector.StationsBegin();
       sdIt != sDetector.StationsEnd(); ++sdIt)
     sEvent.GetStation(sdIt->GetId()).GetSimData().ClearParticleList();

  const CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  const Point& corePosition = simShower.GetPosition();
  const CoordinateSystemPtr 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;

  ShowerParticleIterator& particleIt = fShowerData->fParticleIt;

  for (; particleIt != fShowerData->fParticlesEnd; ++particleIt) {

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

    const Point& pPosition = particleIt->GetPosition();
    const double pR2ShowerFrame = pPosition.GetRho2(showerCS);

    // particle time relative to the core time
    const double pTime = particleIt->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.
    double pPlaneFrontDelay = pTime - PlaneFrontTime(showerCS, corePosition, pPosition);
    if (pPlaneFrontDelay < -0.1*ns )
      continue;

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

    const double pLogR2ShowerFrame = log(pR2ShowerFrame);

    const double pPhiShowerFrame = pPosition.GetPhi(showerCS);

    const StationPositionMatrix::StationInfoPtrList& sList =
      fShowerData->fStationMatrix.GetStationList(pPhiShowerFrame, pLogR2ShowerFrame);

    if (sList.empty())
      continue;

    if (IsParticleEnergyLow(particleIt->GetType(), particleIt->GetKineticEnergy()))
      continue;

    for (StationPositionMatrix::StationInfoPtrList::const_iterator sIt = sList.begin();
         sIt != sList.end(); ++sIt) {

      const StationInfo& sInfo = **sIt;

      if (!sInfo.IsIn(pPhiShowerFrame, pLogR2ShowerFrame))
        continue;

      const double pWeight =
        (particleIt->GetType() == utl::Particle::eMuon ||
         particleIt->GetType() == utl::Particle::eAntiMuon ||
         particleIt->GetSource() == utl::Particle::eShowerFromMuonDecay ?
         fMuonWeightScale*particleIt->GetWeight():
         particleIt->GetWeight());

      Particle newParticle(*particleIt);
      //Set Parent (same particle at this level, keep this for underground propagation backtracking)
      newParticle.SetParent( newParticle );

      // This decision should be taken elsewhere, after
      // calculating AvgN
      //
      //original      newParticle.SetWeight(1);

      // 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 pWeightDensity = pWeight / sInfo.GetResamplingArea();

      const sdet::Station& dStation = sInfo.GetDetStation();
      const double sThickness = dStation.GetThickness();
      const double sRadius = dStation.GetRadius() + sThickness;
      const double sHeight = dStation.GetHeight() + sThickness;
      const CoordinateSystemPtr dStationCS =
        dStation.GetLocalCoordinateSystem();
      Station& station = sInfo.GetEvtStation();

      const unsigned int sId = dStation.GetId();

      const Vector& pDirection = newParticle.GetDirection();
      const double pStationCosine = pDirection.GetZ(dStationCS);
      const double pPlaneCosine = pDirection.GetZ(grParticleCS);
      const double pPlaneAbsCosine = abs(pPlaneCosine);

      // top entry

      // tank top surface seen by the particle?
      if (pPlaneAbsCosine < fHorizontalParticleCut) {
        ++fShowerData->fNHorizontalParticles;
      } else if (pStationCosine < 0 && pPlaneCosine < fHorizontalParticleCut) {

        // projected top area of a tank
        const double effArea = kPi*Sqr((fSimulateAmiga?fAmigaRadius:sRadius)) * pStationCosine/pPlaneCosine;

        const double avgN = pWeightDensity * effArea;

        InsertValue(fShowerData->fWeightStat, sId, avgN);

        unsigned int n;
        double weight;

        if (fUseWeightLimiting) {
          boost::tie(n, weight) = ParticleNumberAndWeight(avgN, sId);

          if (int(weight) == 1)
            ++unityWeightCtr;
          else
            largeWeightCtr += weight;

        } else {
          n = RandPoisson::shoot(fRandomEngine, avgN);
          weight = 1;
        }

        newParticle.SetWeight(weight);

        for (unsigned int i = 0; i < n; ++i) {
          //Different definition of the radius when simulating for AMIGA
          //with the tank and the ground beneath simulated simultaneously
          double r = (fSimulateAmiga?fAmigaRadius:sRadius) * sqrt(RandFlat::shoot(fRandomEngine, 0, 1));
          const double phi = RandFlat::shoot(fRandomEngine, 0, kTwoPi);

          const Point pShiftedPosition(r, phi, sHeight, dStationCS, Point::kCylindrical);
          newParticle.SetPosition(pShiftedPosition);

          // 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
          double pShiftedTime =
            pTime + PlaneFrontTime(showerCS, pPosition, pShiftedPosition);

          if (i && fShowerData->fLogXGauss) {
            const double planeFrontTime =
              PlaneFrontTime(showerCS, corePosition, pShiftedPosition);
            pShiftedTime =
              fShowerData->fLogXGauss->GetSmearedTime(planeFrontTime, pShiftedTime);
          }

          newParticle.SetTime(pShiftedTime);

          InsertValue(fShowerData->fTimeStat, sId, pShiftedTime);

          station.AddParticle(newParticle);

        } // for n

        fShowerData->fTopParticles[sId] += n;
        nParticlesInjected += n;

      } // end of for top entries

      //In the case of AMIGA simulation, side entries have no importance
      //particles are injected onthe top of a "big" cylinder containing both the
      //tank and the ground beneath
      if(!fSimulateAmiga){
        // side entry
        const double pStationSine = pDirection.GetRho(dStationCS);

        if (pPlaneAbsCosine < fHorizontalParticleCut) {

          ++fShowerData->fNHorizontalParticles;

        } else if (pStationSine) {

          // projected side area of a tank
          const double effArea = 2*sRadius*sHeight * pStationSine/pPlaneAbsCosine;
          const double avgN = pWeightDensity * effArea;
          //original      const unsigned int n = RandPoisson::shoot(fRandomEngine, avgN);
          const double pPhi = pDirection.GetPhi(dStationCS) + kPiOnTwo;

          unsigned int n;
          double weight;

          if (fUseWeightLimiting) {
            boost::tie(n, weight) = ParticleNumberAndWeight(avgN, sId);
          } else {
            n = RandPoisson::shoot(fRandomEngine, avgN);
            weight = 1;
          }

          newParticle.SetWeight(weight);

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

            // For the X and Y component we need to find the angle made between
            // the (tank) 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, -sThickness, sHeight);
            const Point pShiftedPosition(sRadius, phi, z, dStationCS, Point::kCylindrical);
            newParticle.SetPosition(pShiftedPosition);

            // Shift the time (read comment above).
            double pShiftedTime =
              pTime + PlaneFrontTime(showerCS, pPosition, pShiftedPosition);
            if (i && fShowerData->fLogXGauss) {
              const double planeFrontTime =
                PlaneFrontTime(showerCS, corePosition, pShiftedPosition);
              pShiftedTime =
                fShowerData->fLogXGauss->GetSmearedTime(planeFrontTime, pShiftedTime);
            }
            newParticle.SetTime(pShiftedTime);

            InsertValue(fShowerData->fTimeStat, sId, pShiftedTime);

            station.AddParticle(newParticle);
          }
          fShowerData->fSideParticles[sId] += n;
          nParticlesInjected += n;
          if (pStationCosine >= 0)
            fShowerData->fUpwardSideParticles[sId] += n;
        }  // for side entries
      }  // for side entries
    } // for station info
  }

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

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

  return eSuccess;
}


void
CachedXShowerRegenerator::OutputStats(Event& event)
{
  const StationPositionMatrix::StationInfoList& sList =
    fShowerData->fStationMatrix.GetStationList();

  int nStations = 0;

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

    TabularStream tab("r|.|r|r|r|.|.|.");
    tab <<              endc << "Rho"  << endc <<          endc <<           endc <<         endc
        << "wght" << endc << "wght" << endc << "wght" << endr

        << "station" << endc << "[km]" << endc << "top" << endc << "side" << endc << "up" << endc
        << "min"  << endc << "avg"  << endc << "max"  << endr << hline;

    for (StationPositionMatrix::StationInfoList::const_iterator sIt = sList.begin();
         sIt != sList.end(); ++sIt) {

      const unsigned sId = sIt->GetDetStation().GetId();

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

        ShowerData::WeightStatMap::const_iterator wIt = fShowerData->fWeightStat.find(sId);
        // ShowerData::TimeStatMap::const_iterator mmIt = fShowerData->fTimeStat.find(sId);
        tab << sId << ' ' << endc
            << setprecision(4) << sIt->GetR()/km << endc
            << ' ' << fShowerData->fTopParticles[sId] << ' ' << endc
            << ' ' << fShowerData->fSideParticles[sId] << ' ' << endc
            << ' ' << fShowerData->fUpwardSideParticles[sId] << ' ' << endc
            << ' ' << setprecision(4) << wIt->second.GetMin() << endc
            << ' ' << setprecision(4) << wIt->second.GetAverage() << endc
            << ' ' << setprecision(4) << wIt->second.GetMax() << 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 ShowerSimData& simShower = event.GetSimShower();
    const TimeStamp& coreTime = simShower.GetTimeStamp();

    for (StationPositionMatrix::StationInfoList::const_iterator sIt = sList.begin();
         sIt != sList.end(); ++sIt) {

      const unsigned sId = sIt->GetDetStation().GetId();

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

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

        // ShowerData::WeightStatMap::const_iterator wIt = fShowerData->fWeightStat.find(sId);
        ShowerData::TimeStatMap::const_iterator mmIt = fShowerData->fTimeStat.find(sId);
        tab << sId << ' ' << endc
            << setprecision(4) << sIt->GetR()/km << endc
            << ' ' << Round(10, mmIt->second.GetMin().GetInterval()/ns) << endc
            << ' ' << Round(10, mmIt->second.GetMax().GetInterval()/ns) << endc
            << ' ' << Round(10, (mmIt->second.GetMin() - sTime).GetInterval()/ns) << endc
            << ' ' << Round(10, (mmIt->second.GetMax() - sTime).GetInterval()/ns) << endr;

      }

    }

    tab << delr;

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

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


VModule::ResultFlag
CachedXShowerRegenerator::Finish()
{
  if (fNOutOfRangeWeight) {
    ostringstream warn;
    warn << "During shower regeneration, a total of " << fNOutOfRangeWeight << " "
            "weighted particle(s) resulted in generation of unity-weight particles exceeding "
         << fWeightLimit << " in number.";
    WARNING(warn);
  }

//   cout << "unityWeightCtr = " << unityWeightCtr << endl;
//   cout << "largeWeightCtr = " << largeWeightCtr << endl;

  return eSuccess;
}


  boost::tuple<unsigned int, double>
  CachedXShowerRegenerator::ParticleNumberAndWeight(const double& avgN, const int& sId)
  {
    unsigned int n;
    double weight;

    if (avgN > fWeightThreshold) {

      //      fShowerData->fWeightCounterMap[sId] += avgN;

      if (fShowerData->fWeightCounterMap[sId] <= fAccumulatedWeightLimit) {

//         if (sId == 5669)
//           cout << "a (avgN = " << avgN << ") ";

        n = RandPoisson::shoot(fRandomEngine, avgN);
        weight = 1;
      } else  {

//         if (sId == 5669)
//           cout << "b (avgN = " << avgN << ") ";

        n = 1;
        weight = avgN;
      }
    }
    else {
      n = RandPoisson::shoot(fRandomEngine, avgN);
      weight = 1;
    }

    fShowerData->fWeightCounterMap[sId] += n;
    return boost::make_tuple(n, weight);
  }

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. CachedXShowerRegenerator/CachedXShowerRegenerator.o -k"
// End: