G4ScintillatorAction.cc

Go to the documentation of this file.

#include "G4ScintillatorAction.h"
#include "G4StationEventAction.h"
#include "G4StationSimulator.h"
#include "G4StationSteppingAction.h"
#include "G4StationTrackingAction.h"
#include "G4Utils.h"

#include <utl/ErrorLogger.h>
#include <utl/RandomEngine.h>

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

#include <G4Event.hh>
#include <G4HCofThisEvent.hh>
#include <G4SDManager.hh>
#include <G4Step.hh>
#include <G4TouchableHistory.hh>
#include <G4Track.hh>
#include <G4Trajectory.hh>
#include <G4TrajectoryContainer.hh>
#include <G4VHitsCollection.hh>

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

#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <G4SystemOfUnits.hh>
#include <G4ThreeVector.hh>

#include <iostream>
#include <fstream>

using namespace fwk;
using namespace std;

using CLHEP::RandPoisson;
using CLHEP::RandExponential;


namespace G4StationSimulatorOG {

  G4ScintillatorAction::G4ScintillatorAction(const G4String& name) :
    G4VSensitiveDetector(name),
    fG4StationSimulator(
      dynamic_cast<G4StationSimulator&>(RunController::GetInstance().GetModule("G4StationSimulatorOG"))
    )
  { }


  G4bool
  G4ScintillatorAction::ProcessHits(G4Step* const step, G4TouchableHistory* const /*rOHist*/)
  {
    // this code uses Offline units, converted only when passed to G4

    const double eDep = step->GetTotalEnergyDeposit() * (utl::eV / CLHEP::eV);

    // if no energy is deposited, no photons (and thus no photoelectrons) will be produced
    if (!eDep)
      return true;  // drop

    const sdet::Station& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();

    /* Unfortunately at the time of writing this comment, the geant detector volumes are only
     * constructed once and re-used for the simulation of different stations. Because of this, it is
     * possible for scintillator volumes to exist in the geant world volume above the water-cherenkov
     * detector even if the station in question does not have a scintillator (meaning there is no
     * sdet/sevt Scintillator). If the detector volumes were to be properly updated between the
     * simulation of each station with Geant, the check below would be unnecessary and could be
     * removed. */
    if (!dStation.HasScintillator())
      return true;  // drop

    const double time = step->GetPreStepPoint()->GetGlobalTime() * (utl::ns / CLHEP::ns);

    // The coordinate system of the Geant4 world volume is identical to the Station coordinate system.
    const utl::CoordinateSystemPtr cs = dStation.GetLocalCoordinateSystem();

    const utl::Point position = To<utl::Point>(step->GetTrack()->GetPosition(), cs, utl::meter/CLHEP::meter);

    const sdet::Scintillator& scin = dStation.GetScintillator();

    // number of PEs for a given reference energy
    const double refPENumber = scin.GetReferencePENumber();
    const double refEnergy = scin.GetReferenceEnergy();

    // fiber attenuation length and parameters describing photon leakage at bar ends
    const double attLength = scin.GetAttenuationLength();
    const double barLeakExpCoeff = scin.GetBarLeakageExpCoefficient();
    const double barLeakMaxRatio = scin.GetBarLeakageMaxRatio();

    // properties governing PE arrival time distribution
    const double effRefractiveIndex = scin.GetEffectiveRefractiveIndex();
    const double fiberExpDecayConst = scin.GetFiberExpDecayConstant();
    const double barExpDecayConst = scin.GetBarExpDecayConstant();

    // distances from point of traversal to PMT along fiber
    const std::vector<double>& fiberDistances = scin.GetDistancesToPMT(position);
    // distances to both ends of the bar traversed (first element is end closer to PMT)
    const std::pair<double, double>& barDistances = scin.GetDistancesToBarEnds(position);

    // for both legs of signal
    for (double d : fiberDistances) {
      // attenuation along fiber for distance to PMT
      const double att = exp(-d / attLength);
      // photon loss at end of bar closer to PMT
      const double leakInner = 1 - barLeakMaxRatio * exp(-barDistances.first / barLeakExpCoeff);
      // photon loss at end of bar further from PMT
      const double leakOuter = 1 - barLeakMaxRatio * exp(-barDistances.second / barLeakExpCoeff);

      const double npeMean = 0.5 * eDep / refEnergy * refPENumber * att * leakInner * leakOuter;
      const auto rand = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector).GetEngine();
      const unsigned int npe = CLHEP::RandPoisson::shoot(rand, npeMean);
      for (unsigned int i = 0; i < npe; ++i) {
        const double scinDelay = CLHEP::RandExponential::shoot(rand, barExpDecayConst);
        const double fiberDelay = CLHEP::RandExponential::shoot(rand, fiberExpDecayConst);
        const double pathDelay = d * effRefractiveIndex / utl::kSpeedOfLight;
        const double peTime = time + scinDelay + fiberDelay + pathDelay;
        // needs to be a better way than hardcoding scint PMT id while maintaining flexibility,
        // but static constants for ids in sdet station are private (and should stay that way)
        fG4StationSimulator.AddPhoton(5, peTime);  // keep
      }
    }

    return true;
  }

}