G4StationSteppingAction.cc

Go to the documentation of this file.

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

#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <sevt/PMT.h>
#include <sevt/PMTSimData.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>
#include <utl/config.h>
#include <utl/ErrorLogger.h>
#include <utl/Particle.h>
#include <utl/Trace.h>
#include <utl/AugerCoordinateSystem.h>
#include <fwk/RunController.h>

#include <G4Step.hh>
#include <G4VProcess.hh>
#include <G4SteppingManager.hh>
#include <G4Track.hh>
#include <G4TrackVector.hh>
#include <G4Trajectory.hh>
#include <G4TrajectoryContainer.hh>

#include <sstream>
#include <vector>

using namespace std;


namespace G4StationSimulatorOG {

  G4StationSteppingAction::RPCTrackIDListMap G4StationSteppingAction::fgRPCTrackIDListMap;
  G4StationSteppingAction::RPCParticleListMap G4StationSteppingAction::fgRPCParticleListMap;

  int G4StationSteppingAction::fgParticleId = 0;


  void
  G4StationSteppingAction::UserSteppingAction(const G4Step* const step)
  {
    const auto& track = *step->GetTrack();
    fgParticleId = track.GetDefinition()->GetPDGEncoding();

    if (!fgParticleId && track.GetNextVolume()) {
      const auto& next = track.GetNextVolume()->GetName();
      if (next == "top" || next == "side" || next == "bottom")
        ++G4StationTrackingAction::fNumBounces;
    }

    if (!fMARTAEnabled)
      return;

    const auto& currentVolume = track.GetVolume()->GetName();
    if (currentVolume != "rpc_gas")
      return;

    const auto& mother = *track.GetTouchable()->GetVolume(1);
    const auto& grandmother = *track.GetTouchable()->GetVolume(2);

    if (mother.GetName() != "rpc_glass" || grandmother.GetName() != "rpc") {
      ostringstream err;
      err << "RPC gas not found where expected : mother is "
          << mother.GetName() << " (should be rpc_glass) and grandmother is "
          << grandmother.GetName() << " (should be rpc)";
      throw utl::InvalidConfigurationException(err.str());
    }

    const G4int rpcId = grandmother.GetCopyNo();

    const int particleID = track.GetDefinition()->GetPDGEncoding();
    const int trackID = track.GetTrackID();

    const auto& prePoint = *step->GetPreStepPoint();
    const double ke1 = prePoint.GetKineticEnergy() / CLHEP::MeV * utl::MeV;
    const double dedx = step->GetTotalEnergyDeposit() / CLHEP::MeV * utl::MeV;
    const double enDepNonIoni = step->GetNonIonizingEnergyDeposit() / CLHEP::MeV * utl::MeV;
    const double enDepIoni = dedx - enDepNonIoni;

    const double w = track.GetWeight();
    const bool isSec = w < 0;

    auto& secondary = *G4UserSteppingAction::fpSteppingManager->GetfSecondary();

    for (const auto sp : secondary) {
      auto& s = *sp;
      if (s.GetVolume()->GetName() != "rpc_gas")
        continue;
      if (!isSec)
        s.SetWeight(-trackID);
    }

    bool trackSave = true;
    auto& trackId = fgRPCTrackIDListMap[rpcId];
    auto& particleList = fgRPCParticleListMap[rpcId];

    if (isSec) {  // Particle it is a secondary, add dEdX to primary

      for (unsigned int t = 0, n = trackId.size(); t < n; ++t) {
        if (trackId[t] != -w)  // find trackID of the primary using W of secondary
          continue;
        if (enDepIoni > 0) {
          auto& p = particleList[t];
          const double enLoss = p.GetWeight();
          p.SetWeight(enDepIoni + enLoss);
        }
        trackSave = false;
      }

    } else {

      for (unsigned int t = 0, n = trackId.size(); t < n; ++t) {
        if (trackId[t] != trackID)  // find trackID of the primary
          continue;
        if (enDepIoni > 0) {
          auto& p = particleList[t];
          const double enLoss = p.GetWeight();
          p.SetWeight(enDepIoni + enLoss);
        }
        trackSave = false;
      }

    }

    // Save a new primary track which is not secondary!
    if (trackSave && !isSec) {

      trackId.push_back(trackID);

      const sdet::Station& station = G4StationSimulator::fgCurrent.GetDetectorStation();
      const auto cs = station.GetLocalCoordinateSystem();

      // Position in Tank CS (centered in center of tank)
      const G4ThreeVector pos = prePoint.GetPosition() - G4StationConstruction::fgTankCenter;  // <-- this is not right
      const double time = prePoint.GetGlobalTime() / CLHEP::ns * utl::ns;
      const auto position = To<utl::Point>(pos, cs, utl::meter/CLHEP::meter);
      const auto direction = To<utl::Vector>(prePoint.GetMomentumDirection(), cs, 1);
      const auto particleType = utl::Particle::Type(particleID);

      utl::Particle& currentParticle = G4StationSimulator::fgCurrent.GetParticle();

      const utl::Particle::Source source = currentParticle.GetSource();
      utl::Particle newParticle(particleType, source, position, direction, time, enDepIoni, ke1);

      // Set parent particle (i.e. the particle entering the tank)

      if (track.GetParentID())
        newParticle.SetParent(currentParticle);

      // Copy production point
      if (currentParticle.HasProductionPoint()) {
        const auto& pProd = currentParticle.GetProductionPoint();
        newParticle.SetProductionPoint(pProd);
      }

      particleList.push_back(newParticle);

    }
  }

}