GGSIntHitSD.cpp

#include "montecarlo/scoring/GGSIntHitSD.h"
#include "montecarlo/pluginmanagers/GGSMCPluginMacros.h"
#include "montecarlo/scoring/GGSIntHit.h"
#include "montecarlo/scoring/GGSIntHitSDMessenger.h"
#include "montecarlo/scoring/GGSTouchableIDComputer.h"
#include "utils/GGSFactory.h"
#include "utils/GGSStringUtils.h"

#include "G4PVReplica.hh"
#include "G4SDManager.hh"
#include "G4Step.hh"
#include "G4TouchableHandle.hh"
#include "G4TransportationManager.hh"
#include "G4VProcess.hh"
#include "G4ios.hh"

#include <boost/unordered_map.hpp>

#include <numeric>

RegisterSD(GGSIntHitSD);

GGSIntHitSD::GGSIntHitSD(G4String name)
    : G4VSensitiveDetector(name), _intHitCollection(NULL), _timeBins(0), _intHitClass("GGSIntHit"),
      _partHitClass("GGSPartHit"), _posHitClass("GGSPosHit"), _storePartHits(false), _storePosHits(false),
      _idComputer(nullptr), _intHitThreshold{0.}, _partHitThreshold{0.}, _posHitThreshold{0.}, _currPartHit{nullptr} {

  collectionName.insert(name);
  _messenger = new GGSIntHitSDMessenger(this);

  auto params = name.substr(name.find_first_of('.') + 1); // Remove volume name
  auto secondPointPos = params.find_first_of('.');
  if (secondPointPos != std::string::npos) {
    // Parameters have been specified so let's handle them
    params = params.substr(secondPointPos + 1); // Remove GGSIntHitSD
    if (params.size() > 0) {
      // Retrieve hit classes
      std::string intHitClassName, partHitClassName, posHitClassName;
      auto start = params.find_first_of('[');
      auto end = params.find_first_of(']');
      if (start < std::string::npos && end != std::string::npos) {
        intHitClassName = params.substr(start + 1, end - start - 1);
      } else {
        throw std::runtime_error(
            std::string("[GGSIntHitSD::GGSIntHitSD] Missing definition of custom integrated hit class in \"") + params +
            ("\""));
      }

      start = params.find_first_of('[', end + 1);
      end = params.find_first_of(']', start + 1);
      if (start != std::string::npos && end != std::string::npos) {
        partHitClassName = params.substr(start + 1, end - start - 1);
      } else {
        throw std::runtime_error(
            std::string("[GGSIntHitSD::GGSIntHitSD] Missing definition of custom particle hit class in \"") + params +
            ("\""));
      }

      start = params.find_first_of('[', end + 1);
      end = params.find_first_of(']', start + 1);
      if (start != std::string::npos && end != std::string::npos) {
        posHitClassName = params.substr(start + 1, end - start - 1);
      } else {
        throw std::runtime_error(
            std::string("[GGSIntHitSD::GGSIntHitSD] Missing definition of custom position hit class in \"") + params +
            ("\""));
      }

      if (intHitClassName.size() > 0) {
        auto &buildersList = GGSIntHitFactory::GetInstance().GetListOfRegisteredBuilders();
        if (std::find(buildersList.begin(), buildersList.end(), intHitClassName) == buildersList.end()) {
          throw std::runtime_error(std::string("[GGSIntHitSD::GGSIntHitSD] Integrated hit class \"") + intHitClassName +
                                   "\" not found.");
        }
        _intHitClass = intHitClassName;
      }
      if (partHitClassName.size() > 0) {
        auto &buildersList = GGSPartHitFactory::GetInstance().GetListOfRegisteredBuilders();
        if (std::find(buildersList.begin(), buildersList.end(), partHitClassName) == buildersList.end()) {
          throw std::runtime_error(std::string("[GGSIntHitSD::GGSIntHitSD] Particle hit class \"") + partHitClassName +
                                   "\" not found.");
        }
        _partHitClass = partHitClassName;
      }
      if (posHitClassName.size() > 0) {
        auto &buildersList = GGSPosHitFactory::GetInstance().GetListOfRegisteredBuilders();
        if (std::find(buildersList.begin(), buildersList.end(), posHitClassName) == buildersList.end()) {
          throw std::runtime_error(std::string("[GGSIntHitSD::GGSIntHitSD] Position hit class \"") + posHitClassName +
                                   "\" not found.");
        }
        _posHitClass = posHitClassName;
      }
    }
  }

  // Create the default ID computer
  SetTouchableIDComputer("GGSUniqueTouchableIDComputer");
}

GGSIntHitSD::~GGSIntHitSD() { delete _messenger; }

void GGSIntHitSD::SetTouchableIDComputer(const std::string &tidcClassName) {
  _idComputer = GGSFactory<GGSTouchableIDComputer>::GetInstance().CreateObject(tidcClassName);
  if (!_idComputer) {
    throw std::runtime_error("[GGSIntHitSD::GGSIntHitSD] Touchable ID computer class " + tidcClassName + " not found");
  }
}

void GGSIntHitSD::Initialize(G4HCofThisEvent *) {

  // The hits collection is added to G4HCofThisEvent in EndOfEvent.
  // It is automatically deleted at the end of the event so it must be recreated every time.
  // However, G4Allocator should provide fast memory allocation for this purpose.
  _intHitCollection = new GGSIntHitsCollection(SensitiveDetectorName, collectionName[0]);
}

G4bool GGSIntHitSD::ProcessHits(G4Step *aStep, G4TouchableHistory *) {
  static const std::string routineName("GGSIntHitSD::ProcessHits");

  // The G4TouchableHistory argument of this method is non-null only if a readout geometry has been
  // defined. Since this feature is not supported, it is currently always null.
  // TODO: investigate the possibility to define a readout geometry

  if (aStep->GetTotalEnergyDeposit() == 0. && !_storePartHits)
    return false;

  int touchableID = _idComputer->ComputeTouchableID(aStep);
  G4TouchableHandle touchable = aStep->GetPreStepPoint()->GetTouchableHandle();
  G4ThreeVector volPos = touchable->GetTranslation(); // Absolute volume translation
  G4VPhysicalVolume *physVol = touchable->GetVolume();

  /* *** Add step to integrated hit *** */
  GGSIntHit *intHit = nullptr;
  bool intFound = false;
  // Search for existing integrated hit for this volume
  auto nIntHits = _intHitCollection->entries();
  for (int iHit = nIntHits - 1; iHit >= 0; iHit--) {
    intHit = (*(_intHitCollection))[iHit];
    if (physVol == intHit->GetVolume() && touchableID == intHit->GetID() && volPos == intHit->GetAbsTranslation()) {
      intFound = true;
      break;
    }
  }
  if (!intFound) {
    // If we reach this point, then this is a new hit
    intHit = GGSIntHitFactory::GetInstance().CreateObject(_intHitClass, _timeBins).release();
    if (!intHit) {
      throw std::runtime_error(
          std::string("[GGSIntHitSD::ProcessHits] Can't create integrated hit of class ").append(_intHitClass));
    }
    // Set global hit properties
    intHit->SetAbsTranslation(volPos);
    intHit->SetID(touchableID);
    intHit->SetVolume(physVol);
    intHit->SetPartHitsStorage(_storePartHits);
    intHit->SetPosHitsStorage(_storePosHits);
    // Call the user initialization routine
    intHit->UserInit(aStep);
    // Insert the new integrated hit into the hit collection
    _intHitCollection->insert(intHit);
  }
  // Add the current step to the hit
  intHit->AddStep(*aStep);

  if (_storePartHits) {
    /* *** Add step to particle hit *** */
    // Create the particle hit if it does not exists, i.e. on first step inside the sensitive volume.
    // If a particle entered in the volume in previous step and immediately got stopped but not killed then
    // IsFirstStepInVolume for the current step will be true, so check also that kinetic energy is non-zero for pre-step
    if (aStep->IsFirstStepInVolume() && aStep->GetPreStepPoint()->GetKineticEnergy() != 0) {
      if (_currPartHit != nullptr) {
        // On first step there must be no current hit
        throw std::runtime_error("GGSIntHitSD::ProcessHit: existing particle hit on first step in sensitive volume");
      }
      _currPartHit = GGSPartHitFactory::GetInstance().CreateObject(_partHitClass).release();
      if (!_currPartHit) {
        throw std::runtime_error(
            std::string("[GGSIntHitSD::ProcessHits] Can't create particle hit of class ").append(_partHitClass));
      }
      // Set global hit properties
      _currPartHit->SetPosHitsStorage(_storePosHits);
    }

    if (_currPartHit != nullptr) {
      // Add the current step to the hit
      _currPartHit->AddStep(*aStep);

      // Check if this is the last step in the volume
      // Do not use aStep->IsLastStepInVolume(); see https://geant4-forum.web.cern.ch/t/last-step-in-a-volume/856
      if ((aStep->GetPostStepPoint()->GetStepStatus() == G4StepStatus::fGeomBoundary &&
           aStep->GetTrack()->GetTrackStatus() != G4TrackStatus::fStopButAlive) ||
          aStep->GetTrack()->GetTrackStatus() == G4TrackStatus::fStopAndKill ||
          aStep->GetTrack()->GetTrackStatus() == G4TrackStatus::fKillTrackAndSecondaries) {
        if (_currPartHit->GetEnergyDeposit() > _partHitThreshold) {
          // Hit above threshold: insert it into the hit collection
          intHit->GetPartHits()->insert(_currPartHit);
        } else {
          // Hit below threshold: delete it
          delete _currPartHit;
        }
        // Do not point at the hit anymore (at this point, it's either stored in the collection or deleted)
        _currPartHit = nullptr;
      }

      if (_storePosHits && aStep->GetTotalEnergyDeposit() > _posHitThreshold) {
        /* *** Create position hit for current step *** */
        GGSPosHit *posHit = GGSPosHitFactory::GetInstance().CreateObject(_posHitClass).release();
        if (!posHit) {
          throw std::runtime_error(
              std::string("[GGSIntHitSD::ProcessHits] Can't create position hit of class ").append(_posHitClass));
        }
        // Insert the new position hit into the hit collection
        auto *partHit =
            (_currPartHit ? _currPartHit : (*(intHit->GetPartHits()))[intHit->GetPartHits()->GetSize() - 1]);
        partHit->GetPosHits()->insert(posHit);
        // Set members of the position hit
        posHit->SetStep(*aStep);
      }
    } else {
      std::vector<std::string> statusNames{
          "fAlive", "fStopButAlive", "fStopAndKill", "fKillTrackAndSecondaries", "fSuspend", "fPostponeToNextEvent"};
      COUT(WARNING) << "No particle hit has been created for particle:\n";
      CCOUT(WARNING) << " - type: "
                     << aStep->GetTrack()->GetDynamicParticle()->GetParticleDefinition()->GetParticleName() << "\n";
      CCOUT(WARNING) << " - creator process: " << aStep->GetTrack()->GetCreatorProcess()->GetProcessName() << "\n";
      CCOUT(WARNING) << " - trackID: " << aStep->GetTrack()->GetTrackID() << "\n";
      CCOUT(WARNING) << " - status: " << statusNames[static_cast<int>(aStep->GetTrack()->GetTrackStatus())] << "\n";
      CCOUT(WARNING) << " - eDep: " << aStep->GetTotalEnergyDeposit() / CLHEP::MeV << " MeV \n";
      CCOUT(WARNING) << " - kin. energy (pre) : " << aStep->GetPreStepPoint()->GetKineticEnergy() / CLHEP::MeV
                     << " MeV \n";
      CCOUT(WARNING) << " - kin. energy (post): " << aStep->GetPostStepPoint()->GetKineticEnergy() / CLHEP::MeV
                     << " MeV \n";
      CCOUT(WARNING) << " - process (pre) : "
                     << (aStep->GetPreStepPoint()->GetProcessDefinedStep()
                             ? aStep->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName()
                             : "none")
                     << "\n";
      CCOUT(WARNING) << " - process (post): "
                     << (aStep->GetPostStepPoint()->GetProcessDefinedStep()
                             ? aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName()
                             : "none")
                     << "\n";
    }
  }

  return true;
}

void GGSIntHitSD::EndOfEvent(G4HCofThisEvent *hitCollection) {

  G4int IHID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
  hitCollection->AddHitsCollection(IHID, _intHitCollection);
}