GGSIntHitSD.cpp

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

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

#include <boost/unordered_map.hpp>

#include <numeric>

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

RegisterSD(GGSIntHitSD);

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

GGSIntHitSD::GGSIntHitSD(G4String name) :
    G4VSensitiveDetector(name), _intHitCollection(NULL), _timeBins(0), _intHitClass("GGSIntHit"), _partHitClass(
        "GGSPartHit"), _posHitClass("GGSPosHit"), _storePartHits(false), _storePosHits(false), _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;
      }
    }
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

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

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

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]);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

G4bool GGSIntHitSD::ProcessHits(G4Step *aStep, G4TouchableHistory*) {

  // 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;

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

  // Build the absolute touchable ID
  // Definition: the multiplicity of a physical volume is defined as the number of that volumes inside a given mother logical volume.
  // Definition: given a history, the multiplicity at depth i is the multiplicity of the volume in history at depth i
  // Definition: given a history, the copy number at depth i is the copy number of the volume in history at depth i
  // Being N_i and M_i the copy number and multiplicity at depth i, the absolute ID of the touchable at history depth 0 is defined as:
  //   ID = N_0 + M_0*N_1 + M_0*M_1*N_2 + ...
  // At each step, the product M_0*M_1*N_2*... is stored in the variable totMult.

  int historyDepth = touchable->GetHistoryDepth();
  int totMult = 1;
  int touchableID = 0;
  // Map to store the multiplicity of each physical volume
  static boost::unordered_map<G4VPhysicalVolume*, int> multMap;

  for (int iDepth = 0; iDepth < historyDepth; iDepth++) { // This loop will go up in history up to the daughter of the world volume (i.e. iDepth == historyDepth - 1)
    touchableID += totMult * touchable->GetCopyNumber(iDepth);
    // Compute total multiplicity for next step
    if (iDepth != historyDepth - 1) { // No need to compute the total multiplicity on last step
      int currVolMult = 0;
      if (!(touchable->GetVolume(iDepth)->IsReplicated())) {
        // See if the multiplicity for the current non-replicated volume has been computed before
        boost::unordered_map<G4VPhysicalVolume*, int>::iterator storedMult = multMap.find(touchable->GetVolume(iDepth));
        if (storedMult == multMap.end()) {
          // Compute the multiplicity for non-replicated volumes placed inside the current mother logical
          G4LogicalVolume *currLogVol = touchable->GetVolume(iDepth)->GetLogicalVolume();
          G4LogicalVolume *motherLogVol = touchable->GetVolume(iDepth)->GetMotherLogical();
          int nDaughters = motherLogVol->GetNoDaughters();
          for (int iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
            if (motherLogVol->GetDaughter(iDaughter)->GetLogicalVolume() == currLogVol)
              currVolMult++;
          }
          // Store the computed multiplicity
          multMap[touchable->GetVolume(iDepth)] = currVolMult;
        }
        else {
          // Use the stored multiplicity
          currVolMult = storedMult->second;
        }
      }
      else {
        // Use the multiplicity of the replicated volume
        EAxis axis;
        G4int nReplicas;
        G4double width, offset;
        G4bool consuming;
        touchable->GetVolume(iDepth)->GetReplicationData(axis, nReplicas, width, offset, consuming);
        currVolMult = nReplicas;
      }
      totMult *= currVolMult;

    }
  }

  /* *** Add step to integrated hit *** */
  GGSIntHit *intHit = nullptr;
  bool intFound = false;
  // Search for existing integrated hit for this volume
  for (int iHit = 0; iHit < _intHitCollection->entries(); 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);
    // 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();
      // Set global hit properties
      _currPartHit->SetPosHitsStorage(_storePosHits);
    }
    // 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);
    }
  }

  return true;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void GGSIntHitSD::EndOfEvent(G4HCofThisEvent *hitCollection) {

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

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....