G4StationStackingAction.cc

Go to the documentation of this file.

#include "G4StationStackingAction.h"
#include "G4StationSimulator.h"

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

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

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

#include <utl/TabulatedFunction.h>

#include <G4Track.hh>
#include <G4TrackingManager.hh>
#include <G4TrackStatus.hh>
#include <G4OpticalPhoton.hh>

using namespace sdet;
using namespace fwk;


namespace G4StationSimulatorOG {

  G4ClassificationOfNewTrack
  G4StationStackingAction::ClassifyNewTrack(const G4Track* const track)
  {
    // NOTE: for now, we assume that each pmt has the same collection
    // and quantum efficiencies

    if (track->GetDefinition() != G4OpticalPhoton::OpticalPhoton())
      return fWaiting;  // let G4 handle everything but photons

    // optical photons at this point

    const sevt::Station& station = G4StationSimulator::fgCurrent.GetEventStation();

    // if SimulatorSignature is G4StationSimulatorFullTrackOG, the quantum and
    // collection efficiencies are taken into account in G4TankPMTAction.cc
    // (for each PMT type accordingly), ie all photons are tracked by G4 until
    // they hit one of the PMTs (or get absorbed)

    const auto signature = station.GetSimData().GetSimulatorSignature();
    if (signature == "G4StationSimulatorFullTrackOG")
      return fWaiting;  // let G4 handle all the photons

    // if SimulatorSignature is G4StationSimulatorFullOG randomly kill photons
    // according to a maximal envelope of the quantum and collection
    // efficiencies of PMTs to speed up simulation, currently only two
    // different PMTs: large (id=1) and small (id=4) WCD PMT

    const double energy = track->GetKineticEnergy() * (utl::MeV / CLHEP::MeV);

    // large PMT
    const auto& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();
    const auto& pmt = dStation.GetPMT(1);
    const double eff =
      pmt.GetCollectionEfficiency() * pmt.GetQuantumEfficiency().InterpolateY(energy, 1);

    const auto rand = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector).GetEngine();
    const double keep = CLHEP::RandFlat::shoot(rand, 0, 1);

    if (!dStation.HasSmallPMT())
      return (keep > eff) ? fKill : fWaiting;

    // small PMT
    const auto& pmtSmall = dStation.GetPMT(4);
    const double effSmall =
      pmtSmall.GetCollectionEfficiency() * pmtSmall.GetQuantumEfficiency().InterpolateY(energy, 1);

    return (keep > std::max(eff, effSmall)) ? fKill : fWaiting;
  }

}