G4UMDScintStripAction.cc

Go to the documentation of this file.

#include "G4UMDScintStripAction.h"
#include "G4StationSimulator.h"
#include "G4StationTrackingAction.h"
#include "G4Utils.h"

#include <utl/Particle.h>
#include <utl/ParticleCases.h>

#include <utl/Segment.h>
#include <utl/Plane.h>
#include <utl/Point.h>
#include <utl/GeometryUtilities.h>
#include <utl/PhysicalConstants.h>

#include <mdet/MDetector.h>
#include <mdet/Counter.h>

#include <G4Step.hh>
#include <G4TouchableHistory.hh>
#include <G4ios.hh>
#include <G4SystemOfUnits.hh>
#include <G4RunManager.hh>
#include <G4OpticalPhoton.hh>


namespace G4StationSimulatorOG {

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


  G4bool
  G4UMDScintStripAction::ProcessHits(G4Step* const step, G4TouchableHistory* const /*rOhist*/)
  {
    // If it is an opticalphoton, continue
    if (step->GetTrack()->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition())
      return true;

    // Retrieve current physical volumen and its name (i.e. the current physical strip)
    G4VPhysicalVolume* const currentPhysVol = step->GetPreStepPoint()->GetPhysicalVolume();
    // Retrieve mother volume of the current physical volumen and its name (i.e. the current physical module)
    G4VPhysicalVolume* const motherPhysVol = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume(1);
    //const G4String& moduleName = motherPhysVol->GetName();

    const G4int stripId = currentPhysVol->GetCopyNo();
    const G4int moduleId = motherPhysVol->GetCopyNo();

    // Particle
    const int particleId = step->GetTrack()->GetDefinition()->GetPDGEncoding();
    const double stepLen = step->GetStepLength() * (utl::m / CLHEP::m);

    const mdet::Counter& dCounter = G4StationSimulator::fgCurrent.GetDetectorUMDCounter();

    const utl::CoordinateSystemPtr cs = dCounter.GetLocalCoordinateSystem();

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

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

    // Only first step inside scintillator (PreStepPoint() == fGeomBoundary on:
    //
    // vol1  |         vol2 (Scint)               |      vol3
    //       |                                    |
    //       |     CurrStep                       |
    //       |o----------------o                  |
    //       |                                    |
    //       |PrePoint         PostPoint          |
    //       |                                    |
    //       |                                    |
    //       |                                    |
    //
    // GeomBoundary Vol1-Vol2               GeomBoundary Vol2-Vol3

    if (step->GetPreStepPoint()->GetStepStatus() == fGeomBoundary) {

      G4StationTrackingAction* const trackingA = (G4StationTrackingAction*)G4RunManager::GetRunManager()->GetUserTrackingAction();
      const auto type = trackingA->GetFirstStepInVolume() ? utl::Particle::eShower : utl::Particle::eCornerClipping;

      const double kineticEnergy = step->GetTrack()->GetKineticEnergy() * (utl::MeV / CLHEP::MeV);
      const auto direction = To<utl::Vector>(step->GetPreStepPoint()->GetMomentumDirection(), cs, 1);
      const utl::Particle injectedParticle(particleId, type, position, direction, time, 1, kineticEnergy);

      // Skip non muon particles if desired (useful to study corner-clipping impact by particles other than muons)
       if (fG4StationSimulator.OnlyMuonsInUMD() &&
           particleId != utl::Particle::eMuon &&
           particleId != utl::Particle::eAntiMuon)
         return false;

      fG4StationSimulator.AddInjectedParticle(moduleId, stripId, injectedParticle);
      // Allows to keep memory of the first time the current particle enters in the current strip
      trackingA->SetFirstStepInVolume(false);

    }

    if (fG4StationSimulator.fUMDFastMode) {

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

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

      const mdet::Module& detMod = dCounter.GetModule(moduleId);
      const mdet::Scintillator& detSci = detMod.GetScintillator(stripId);
      const mdet::Fiber& fiber = detMod.GetFiberFor(detSci);

      auto& counter = G4StationSimulator::fgCurrent.GetEventUMDCounter();

      if (!counter.HasModule(moduleId))
        counter.MakeModule(moduleId);

      mevt::Module& evtMod = counter.GetModule(moduleId);

      if (!evtMod.HasScintillator(stripId))
        evtMod.MakeScintillator(stripId);

      mevt::Scintillator& evtScint = evtMod.GetScintillator(stripId);

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

      mevt::ScintillatorSimData& simData = evtScint.GetSimData();

      if (!simData.HasAnalogicTrace())
        simData.MakeAnalogicTrace();

      simData.AddEnergyDeposit(depE);
      if (particleId == utl::Particle::eMuon || particleId == utl::Particle::eAntiMuon)
        simData.AddEnergyDepositMuons(depE);

      // Compute a plane perpendicular to the fiber.
      const auto& fcs = fiber.GetLocalCoordinateSystem();
      const utl::Plane targetPlane(utl::Point(0, 0, 0, fcs), utl::Vector(1, 0, 0, fcs));

      // The on-scintillator distance is the distance from the current point to the plane
      // Notice this distance is signed; can be negative value (depending on the normal-to-the-plane orientation).
      const double onScint = std::abs(utl::Distance(position, targetPlane));

      // Extra fiber length to arrive to the optical sensor (either PMT or SiPM).
      const double onManifold = fiber.GetOnManifoldLength();

      // Total distance.
      const double pathLen = onScint + onManifold;

      // Compute the number of photo-electron from attenuation curve
      unsigned int nSPE = fiber.ComputeSPENumber(pathLen, stepLen, depE);

      // These are the same for every spe:
      const double pathDelay = pathLen * fiber.GetRefractionIndex() / utl::kSpeedOfLight;

      // Now the arrival time
      const double particleDelay = time.GetInterval();

      for (unsigned int n = 0; n < nSPE; ++n) {
        const double delay =
          pathDelay +
          fiber.ComputeDecayDelay() +
          detSci.ComputeDecayDelay() +
          particleDelay;
        fG4StationSimulator.AddUMDPhoton(moduleId, stripId, delay);
      }
    }

    return true;
  }

}