G4StationSimulator/G4TankPMTAction.cc

Go to the documentation of this file.

#include "G4TankPMTAction.h"
#include "G4StationEventAction.h"
#include "G4StationSimulator.h"
#include "G4StationSteppingAction.h"

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

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

#include <utl/ErrorLogger.h>
#include <utl/Particle.h>
#include <utl/ParticleCases.h>
#include <utl/Trace.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 <iostream>

using namespace std;
using namespace sevt;
using namespace sdet;
using namespace det;
using namespace utl;
using namespace fwk;


namespace G4StationSimulatorOG {

  G4TankPMTAction::G4TankPMTAction(const G4String& name, const int pmtIndex) :
    G4VSensitiveDetector(name),
    fPMTIndex(pmtIndex),
    fG4StationSimulator(
      dynamic_cast<G4StationSimulator&>(RunController::GetInstance().GetModule("G4StationSimulatorOG"))
    )
  {
    if (pmtIndex < 1 || pmtIndex > 5) {
      FATAL("G4PMTAction() : invalid PMT index.");
      exit(EXIT_FAILURE);
    }
  }


  G4bool
  G4TankPMTAction::ProcessHits(G4Step* const step, G4TouchableHistory* const /*rOHist*/)
  {
    // use to reject hits registered in sensitive volume in case the particle is not a photon
    const int particle = G4StationSteppingAction::GetCurrentParticleId();
    if (particle)
      return true;  // kill

    const double time = step->GetPreStepPoint()->GetGlobalTime() * (utl::second / CLHEP::second);
    if (time >= 1*CLHEP::second)
      return true;  // kill

    // sometimes particles (PDG ID=0, photons?) with crazy energy
    const double energy = step->GetPreStepPoint()->GetKineticEnergy() * (utl::MeV / CLHEP::MeV);

    auto& station = G4StationSimulator::fgCurrent.GetEventStation();
    const auto& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();

    if (station.GetSimData().GetSimulatorSignature() == "G4StationSimulatorFullTrackOG") {

      // if SimulatorSignature is G4StationSimulatorFullTrackOG randomly kill
      // photons with correct quntum and collection efficiency for each PMT type

      const auto& pmt = dStation.GetPMT(fPMTIndex);
      const auto& qeff = pmt.GetQuantumEfficiency();
      if (energy < qeff.XFront() || energy > qeff.XBack())
        return true;  // kill

      const double eff = pmt.GetCollectionEfficiency() * qeff.InterpolateY(energy, 1);
      const auto rand = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector).GetEngine();
      if (CLHEP::RandFlat::shoot(rand, 0, 1) > eff)
        return true;  // kill

    } else {

      // if not random killing is done already in StationStackingAction.cc to
      // speed up simulation (standard way); here only the difference between
      // the maximal efficiency envelope and the individual PMTs is considered

      if (dStation.HasSmallPMT()) {

        const auto& pmt = dStation.GetPMT(fPMTIndex);
        const auto& qeff = pmt.GetQuantumEfficiency();
        if (energy < qeff.XFront() || energy > qeff.XBack())
          return true;  // kill
        const double eff = pmt.GetCollectionEfficiency() * qeff.InterpolateY(energy, 1);

        const auto& pmtLarge = dStation.GetPMT(1);
        const double effLarge =
          pmtLarge.GetCollectionEfficiency() * pmtLarge.GetQuantumEfficiency().InterpolateY(energy, 1);
        const auto& pmtSmall = dStation.GetPMT(4);
        const double effSmall =
          pmtSmall.GetCollectionEfficiency() * pmtSmall.GetQuantumEfficiency().InterpolateY(energy, 1);
        const double effMax = std::max(effSmall, effLarge);
        if (effMax < 1e-5)
          return true;  // kill

        const auto rand = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector).GetEngine();
        if (effMax * CLHEP::RandFlat::shoot(rand, 0, 1) > eff)
          return true;  // kill
      }

    }

    fG4StationSimulator.AddPhoton(fPMTIndex, time);  // count the photon
    // note that in the current simplified simulation, the photons hitting the pmt are always absorbed:
    return true;  // kill
  }

}