SdAccidentalInjector.cc

Go to the documentation of this file.

#include "SdAccidentalInjector.h"

#include <det/Detector.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>

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

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>

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

#include <utl/Math.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/Particle.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/ErrorLogger.h>
#include <utl/Vector.h>
#include <utl/Cache.h>
#include <utl/CoordinateSystemPtr.h>

#include <TH2D.h>
#include <TFile.h>

#include <cstddef>
#include <sstream>
#include <vector>

using namespace SdAccidentalInjectorKG;
using namespace fwk;
using namespace utl;
using namespace std;


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

  // Beware: We assume that all stations are identical
  const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();
  const double radius = sDetector.StationsBegin()->GetRadius();
  const double height = sDetector.StationsBegin()->GetHeight();
  const double hull = sDetector.StationsBegin()->GetThickness();
  fCyl = CylinderSurfaceSampler(radius, height, hull);

  topBranch.GetChild("muonRatePerStation").GetData(fMuonRatePerStation);
  // these times are relative to core time
  topBranch.GetChild("timeBefore").GetData(fTimeBefore);
  topBranch.GetChild("timeAfter").GetData(fTimeAfter);
  topBranch.GetChild("maxDistance").GetData(fMaxDistance2);
  fMaxDistance2 *= fMaxDistance2;

  Branch bFlux = topBranch.GetChild("fluxHistogram");
  const auto nCosThetaColumns = bFlux.GetChild("nCosThetaColumns").Get<unsigned int>();
  const auto nLog10EGeVRows = bFlux.GetChild("nLog10EGeVRows").Get<unsigned int>();
  const auto rangeCosTheta = bFlux.GetChild("rangeCosTheta").Get<vector<double>>();
  const auto rangeLog10EGeV = bFlux.GetChild("rangeLog10EGeV").Get<vector<double>>();
  const auto table = bFlux.GetChild("table").Get<vector<double>>();

  const double minimumEnergy = topBranch.GetChild("minimumEnergy").Get<double>();
  const double lgEmin = log10(minimumEnergy / GeV);
  const int iemin = max(0., floor((lgEmin - rangeLog10EGeV.at(0)) / (rangeLog10EGeV.at(1) - rangeLog10EGeV.at(0)) * nLog10EGeVRows));
  if (iemin >= int(nLog10EGeVRows)) {
    ERROR("minimumEnergy cuts the whole histogram away!");
    return eFailure;
  }
  const double lgEminEdge = rangeLog10EGeV[0] + iemin * (rangeLog10EGeV[1] - rangeLog10EGeV[0]) / nLog10EGeVRows;

  // we only need histogram above minimum energy
  TH2D h("SdAccidentalInjector::fHist", "",
         nLog10EGeVRows - iemin, lgEminEdge, rangeLog10EGeV[1],       // x
         nCosThetaColumns, rangeCosTheta.at(0), rangeCosTheta.at(1)); // y

  // muon rate per station =
  //   integral (muon_flux * effective_station_area dlogE dCosTheta dPhi)
  // only relative intensities are important for sampling,
  // so constant factors are skipped
  int i = 0;
  for (unsigned int ie = 1; ie <= nLog10EGeVRows; ++ie) {
    //const double lgE = h.GetXaxis().GetBinCenter(ie);
    for (unsigned int ic = 1; ic <= nCosThetaColumns; ++ic) {
      const double cosTh = h.GetYaxis()->GetBinCenter(ic);
      const double area = StationArea(radius, height, cosTh);
      const double flux = table.at(i++);
      if (int(ie) >= iemin)
        h.SetBinContent(ie - iemin, ic, flux * area);
    }
  }

  fHist = HistogramSampler(h);

  ostringstream msg;
  msg << "Configuration:\n"
         "Muon rate per station [Hz]             : " << fMuonRatePerStation/(1./s) << "\n"
         "Minimum energy [GeV]                   : " << minimumEnergy/GeV << "\n"
         "Start injection [ns, rel. to core time]: " << -fTimeBefore/ns << "\n"
         "Stop  injection [ns, rel. to core time]: " << fTimeAfter/ns;
  INFO(msg);

  return eSuccess;
}


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

  const Point& core = event.GetSimShower().GetPosition();

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

  ostringstream msg;
  msg << "Injecting accidental muons:";
  int nTotal = 0;
  for (sdet::SDetector::StationIterator sdIt = sDetector.StationsBegin();
       sdIt != sDetector.StationsEnd(); ++sdIt) {

    const sdet::Station& dStation = *sdIt;

    if ((dStation.GetPosition() - core).GetMag2() > fMaxDistance2)
      continue;

    const int nMuons = fPoisson(fMuonRatePerStation * (fTimeAfter + fTimeBefore));

    if (!nMuons)
      continue;

    if (!sEvent.HasStation(dStation.GetId()))
      sEvent.MakeStation(dStation.GetId());
    sevt::Station& station = sEvent.GetStation(dStation.GetId());

    if (!station.HasSimData())
      station.MakeSimData();

    const CoordinateSystemPtr dStationCS = dStation.GetLocalCoordinateSystem();

    msg << "\n"
           "station " << dStation.GetId() << " gets " << nMuons << " muon" << (nMuons == 1 ? "" : "s") << ':';

    for (int iMu = 0; iMu < nMuons; ++iMu) {

      const auto ec = fHist(); // returns: lg(E/GeV), cos(th)
      const double energy = pow(10, ec.first) * GeV;
      const double theta = acos(ec.second);
      const double phi = fFlat(2 * kPi);

      const double totalMomentum = sqrt(energy*energy - kMuonMass*kMuonMass);
      const Vector pMomentum(-totalMomentum, theta, phi, dStationCS, Vector::kSpherical);
      const TimeInterval pTime = fFlat(-fTimeBefore, fTimeAfter);

      const tuple<double, double, double> xyz = fCyl(theta, phi);
      const Point pPosition(get<0>(xyz), get<1>(xyz), get<2>(xyz), dStationCS);

      const Particle particle(Particle::eMuon, Particle::eBackground, pPosition, pMomentum, pTime, 1);

      station.AddParticle(particle);

      msg << "\n"
             "  energy = " << (energy/GeV) << " GeV, "
             "theta = " << theta/degree << " deg, "
             "time = " << pTime.GetInterval() / nanosecond << " ns";

    }

    nTotal += nMuons;

  }

  msg << "\n"
         "Total number of accidental muons: " << nTotal;
  INFO(msg);

  return eSuccess;
}