GGSHitsAction.cpp

Go to the documentation of this file.

/*
 * GGSHitsAction.cpp
 *
 *  Created on: Oct 8, 2010
 *      Author: Elena Vannuccini and Nicola Mori
 */

#include "montecarlo/useractions/GGSHitsAction.h"
#include "montecarlo/pluginmanagers/GGSMCPluginMacros.h"
#include "montecarlo/scoring/GGSIntHitSD.h"
#include "montecarlo/scoring/GGSMultiSensitiveDetector.h"
#include "montecarlo/services/GGSRootFileService.h"
#include "utils/GGSNameDecoder.h"
#include "utils/GGSSmartLog.h"
#include "utils/GGSStringUtils.h"

#include "G4SDManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UnitsTable.hh"

#include "TArrayF.h"

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

RegisterUA(GGSHitsAction);

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

GGSHitsAction::GGSHitsAction()
    : GGSUserAction(), _messenger(this, "/GGS/userActions/hitsAction/", "Commands for hits action"), _outFileBase(""),
      _outFile(nullptr), _outTree(nullptr), _outTreeName(""), _privateMessenger(this) {

  // Setup the messenger
  _messenger.DeclareProperty("fileBase", _outFileBase)
      .SetGuidance("Sets the base name for ROOT output file")
      .SetGuidance("   Can be with or without extension (.root will be used automatically)")
      .SetStates(G4State_PreInit, G4State_Idle);
  _messenger.DeclareProperty("treeName", _outTreeName)
      .SetGuidance("Set the name of the TTree object in the ROOT output file.")
      .SetStates(G4State_PreInit, G4State_Idle);

  _messenger.DeclareMethod("outputIntHitClass", &GGSHitsAction::SetOutputIntHitClass)
      .SetGuidance("Class name of output integrated hits")
      .SetStates(G4State_PreInit, G4State_Idle);
  _messenger.DeclareMethod("outputPartHitClass", &GGSHitsAction::SetOutputPartHitClass)
      .SetGuidance("Class name of output particle hits")
      .SetStates(G4State_PreInit, G4State_Idle);
  _messenger.DeclareMethod("outputPosHitClass", &GGSHitsAction::SetOutputPosHitClass)
      .SetGuidance("Class name of output position hits")
      .SetStates(G4State_PreInit, G4State_Idle);
}

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

GGSHitsAction::~GGSHitsAction() {
  // Let the TObjArry disown their TClonesArrays to avoid a double deletion since they are owned by the
  // GGSTClonesArrayService.
  for (auto &handleElem : _outputHandles) {
    OutputHandle &handle = handleElem.second;
    if (handle.partHitArray) {
      handle.partHitArray->SetOwner(false);
    }
    if (handle.posHitArray) {
      handle.posHitArray->SetOwner(false);
    }
  }
}

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

void GGSHitsAction::BeginOfRunAction(const G4Run *run) {
  const std::string routineName("GGSHitsAction::BeginOfRunAction");

  // -----------------
  // Open the run file
  // -----------------
  _outFile = GGSRootFileService::GetInstance().GetFileForThisRun(_outFileBase, run);
  if (!_outFile || (_outFile && _outFile->IsZombie())) {
    GGSCOUT(ERROR) << "Cannot open ROOT file " << _outFileBase << "for run" << run->GetRunID() << GGSENDL;
    throw -1;
  }

  // -----------------
  // Create a new tree
  // -----------------
  if (_outTreeName == "") {
    _outTree = GGSRootFileService::GetInstance().GetDefaultTree(_outFile);
    _outTree->SetTitle(TString(_outTree->GetTitle()) + "Hits ");
  } else
    _outTree = new TTree(_outTreeName, "GGS hits");

  // Loop over logical volumes...
  G4LogicalVolumeStore *logicalVS = G4LogicalVolumeStore::GetInstance();
  if (!logicalVS)
    return;

  std::vector<G4VSensitiveDetector *> processedSDs;
  for (unsigned int i = 0; i < logicalVS->size(); i++) {
    G4LogicalVolume *logicalV = logicalVS->at(i);

    // ...and select sensitive volumes with GGS hits SD
    std::vector<G4VSensitiveDetector *> sDets(0);
    G4VSensitiveDetector *sd = logicalV->GetSensitiveDetector();
    if (sd) {
      GGSMultiSensitiveDetector *msd = dynamic_cast<GGSMultiSensitiveDetector *>(sd);
      if (msd) {
        sDets = msd->GetListOfSensitiveDetectors();
      } else {
        sDets.push_back(sd);
      }
    }

    for (auto &sdPtr : sDets) {
      if (std::find(processedSDs.begin(), processedSDs.end(), sdPtr) != processedSDs.end()) {
        continue;
      }
      auto sdName = sdPtr->GetName();
      auto ggssdNamePos = sdName.find(".GGSIntHitSD", 0);
      if (ggssdNamePos != std::string::npos) {
        GGSIntHitSD *intHitSD = dynamic_cast<GGSIntHitSD *>(sdPtr);
        if (intHitSD && intHitSD->isActive()) {
          std::string sdNameReduced = sdName;
          sdNameReduced.erase(ggssdNamePos, 12); // Remove ".GGSIntHitSD"

          std::unordered_map<std::string, OutputHandle>::iterator handleIter = _outputHandles.find(sdName);
          if (handleIter == _outputHandles.end()) {
            // The handle for this detector does not exist yet, so create a standard one

            auto insertResult = _outputHandles.insert(std::pair<std::string, OutputHandle>{sdName, OutputHandle()});
            if (insertResult.second) {
              handleIter = insertResult.first;
            }
          }

          // Create arrays et al.
          auto &newHandle = handleIter->second;
          newHandle.intHitArray = std::unique_ptr<TClonesArray>(new TClonesArray(
              newHandle.intHitClassName.c_str())); // Create it here and Clear("C") it at beginning of event
          // Create a branch to store hits in the output file
          _outTree->Branch((sdNameReduced + "-GGSIntHits").data(), "TClonesArray", newHandle.intHitArray.get(), 64000);

          // Set output threshold
          {
            auto thresholdIter = std::find_if(
                _thresholds.begin(), _thresholds.end(), [&sdNameReduced](const ThresholdStruct &threshold) {
                  return (threshold.hitType == 0) && (sdNameReduced == threshold.detectorName);
                });
            if (thresholdIter != _thresholds.end()) {
              newHandle.intHitThreshold = thresholdIter->value;
            }
          }
          if (intHitSD->GetPartHitsStorage()) {
            newHandle.partHitArray = std::unique_ptr<TObjArray>(new TObjArray());
            // Set ownership to avoid memory leaks when reading back the TObjArray.
            // http://root.cern.ch/root/roottalk/roottalk02/0444.html
            newHandle.partHitArray->SetOwner(true);
            _outTree->Branch((sdNameReduced + "-GGSPartHits").data(), "TObjArray", newHandle.partHitArray.get(), 64000);
            newHandle.partHitCAService =
                std::unique_ptr<GGSTClonesArrayService>(new GGSTClonesArrayService(newHandle.partHitClassName.c_str()));
            // Set output threshold
            auto thresholdIter = std::find_if(
                _thresholds.begin(), _thresholds.end(), [&sdNameReduced](const ThresholdStruct &threshold) {
                  return (threshold.hitType == 1) && (sdNameReduced == threshold.detectorName);
                });
            if (thresholdIter != _thresholds.end()) {
              newHandle.partHitThreshold = thresholdIter->value;
            }
          }
          if (intHitSD->GetPosHitsStorage()) {
            newHandle.posHitArray = std::unique_ptr<TObjArray>(new TObjArray());
            newHandle.posHitArray->SetOwner(true);
            _outTree->Branch((sdNameReduced + "-GGSPosHits").data(), "TObjArray", newHandle.posHitArray.get(), 64000);
            newHandle.posHitCAService =
                std::unique_ptr<GGSTClonesArrayService>(new GGSTClonesArrayService(newHandle.posHitClassName.c_str()));
            // Set output threshold
            auto thresholdIter = std::find_if(
                _thresholds.begin(), _thresholds.end(), [&sdNameReduced](const ThresholdStruct &threshold) {
                  return (threshold.hitType == 2) && (sdNameReduced == threshold.detectorName);
                });
            if (thresholdIter != _thresholds.end()) {
              newHandle.posHitThreshold = thresholdIter->value;
            }
          }

          // Save the time bins as run products
          if (GGSRootFileService::GetInstance().GetOutputMode() == GGSRootFileService::OutputMode::MULTIFILE ||
              G4Threading::G4GetThreadId() <= 0) {
            unsigned int nBins = intHitSD->GetTimeBins().size();
            TArrayF timeBins(nBins);
            for (unsigned int iBin = 0; iBin < nBins; iBin++)
              timeBins[iBin] = intHitSD->GetTimeBins()[iBin];
            _outFile->WriteObject(&timeBins, (sdNameReduced + "-GGSIntHits-TimeBins").data());

            // Save the hit thresholds
            std::vector<float> detThresholds{(float)(newHandle.intHitThreshold / GeV),
                                             (float)(newHandle.partHitThreshold / GeV),
                                             (float)(newHandle.posHitThreshold / GeV)};
            _outFile->WriteObject(&detThresholds, (sdNameReduced + "-GGSIntHits-Thresholds").data());
          }

        } // Check on active volume

      } // Check on sensitive volume name
      processedSDs.push_back(sdPtr);
    } // End loop on sensitive detectors

  } // End loop on logical volumes
}

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

void GGSHitsAction::EndOfRunAction(const G4Run *run) {

  // -----------------
  // close output file
  // -----------------

  if (_outFile) {
    _outFile->cd();
    if (_outTreeName != "") {
      // Write the tree if we are not using the default one
      _outTree->Write();
    }

    GGSRootFileService::GetInstance().CloseFileForThisRun(_outFileBase,
                                                          run); // Will automatically save detector informations
    _outFile = nullptr;
    _outTree = nullptr; // The tree is delete by the file's destructor.
  } else {
    throw -1; // TODO handle errors
  }

  // -----------------
  // delete detectors
  // -----------------

  // Prepare for an eventual next run
  _GiveBackAllArrays();
  for (auto &handleElem : _outputHandles) {
    if (handleElem.second.posHitCAService) {
      handleElem.second.posHitCAService->DeleteAll();
    }
    if (handleElem.second.partHitCAService) {
      handleElem.second.partHitCAService->DeleteAll();
    }
  }
}

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

void GGSHitsAction::BeginOfEventAction(const G4Event *) {

  // Give back arrays to the TCLonesArray services.
  _GiveBackAllArrays();

  for (auto &handleElem : _outputHandles) {
    OutputHandle &handle = handleElem.second;
    if (handle.partHitArray) {

      // Clear the TObjarray persistence structures. These contains TClonesArrays which
      // must not be deleted since they are owned and managed by the GGSTClonesArrayService...

      // ... so set the TObjArray to not own them...
      handle.partHitArray->SetOwner(false);
      // ... then clean to make size=0 without deleting the contained objects...
      handle.partHitArray->Clear();
      // ... and then restore ownership to avoid memory leak on readout of ROOT file.
      // (see http://root.cern.ch/root/roottalk/roottalk02/0444.html)
      handle.partHitArray->SetOwner(true);
    }

    if (handle.posHitArray) {
      handle.posHitArray->SetOwner(false);
      handle.posHitArray->Clear();
      handle.posHitArray->SetOwner(true);
    }
  }
}

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

void GGSHitsAction::EndOfEventAction(const G4Event *event) {

  //-------------------------------
  // Set hits info
  //-------------------------------

  // Retrieve SD manager
  G4SDManager *SDman = G4SDManager::GetSDMpointer();

  // Retrieve event hits
  G4HCofThisEvent *evHitColl = event->GetHCofThisEvent();

  // Loop on detectors
  for (auto &handleElem : _outputHandles) {

    const std::string &detectorName = handleElem.first;
    OutputHandle &handle = handleElem.second;
    // Get hits of current detector
    G4int hitsCollID = SDman->GetCollectionID(detectorName);

    // Convert hits to persistent data format
    handle.intHitArray->Clear(); // Set entries to 0 but do not clear elements (will be cleared in _Convert)
    _Convert(*((GGSIntHitsCollection *)(evHitColl->GetHC(hitsCollID))), *(handle.intHitArray), handle, detectorName);

    // Link particle and position hit collections to persistence array
    for (int iIntHit = 0; iIntHit < handle.intHitArray->GetEntries(); ++iIntHit) {
      GGSTIntHitBase *tIntHit = (GGSTIntHitBase *)(handle.intHitArray->At(iIntHit));
      if (tIntHit->_partHits) {
        handle.partHitArray->AddAtAndExpand(tIntHit->_partHits, handle.partHitArray->GetEntries());
        if (tIntHit->_partHits->GetEntries() > 0 && ((GGSTPartHitBase *)(tIntHit->_partHits->At(0)))->_posHits) {
          static TClonesArray *partHits = nullptr;
          static Int_t nPartHits = 0;
          partHits = tIntHit->_partHits;
          nPartHits = partHits->GetEntries();
          for (Int_t iPartHit = 0; iPartHit < nPartHits; iPartHit++) {
            GGSTPartHitBase *partHit = (GGSTPartHitBase *)(tIntHit->_partHits->At(iPartHit));
            partHit->_posHitIndex = handle.posHitArray->GetEntries(); // Save reference to pos hit array
            handle.posHitArray->AddAtAndExpand(partHit->_posHits, handle.posHitArray->GetEntries());
          }
        }
      }
    }
  }

  //-----------
  // Fill tree
  //-----------
  if (_outTreeName != "") {
    // Call fill only if not using the default tree
    _outFile->cd();
    _outTree->Fill();
  }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void GGSHitsAction::_Convert(const GGSPosHitsCollection &posHits, TClonesArray &tPosHits, OutputHandle &handle) {

  if (tPosHits.GetEntries() != 0) {
    throw std::runtime_error("GGSHitsAction::_Convert: dirty TClonesArray for position hits");
  }
  for (size_t iHit = 0, currTPosHit = 0; iHit < posHits.GetSize(); ++iHit) {
    auto &posHit = *((GGSPosHit *)(posHits.GetHit(iHit)));
    if (posHit.GetEnergyDeposit() > handle.posHitThreshold) {
      GGSTPosHitBase &tPosHit = *((GGSTPosHitBase *)(tPosHits.ConstructedAt(currTPosHit, "C"))); // Get a cleared hit
      tPosHit.eDep = posHit.GetEnergyDeposit() / GeV;
      tPosHit.UserConversion(posHit);
      // Increase counter for output position hits
      ++currTPosHit;
    }
  }
}

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

void GGSHitsAction::_Convert(const GGSPartHitsCollection &partHits, TClonesArray &tPartHits, OutputHandle &handle) {

  if (tPartHits.GetEntries() != 0) {
    throw std::runtime_error("GGSHitsAction::_Convert: dirty TClonesArray for particle hits");
  }
  for (size_t iHit = 0, currTPartHit = 0; iHit < partHits.GetSize(); ++iHit) {
    auto &partHit = *((GGSPartHit *)(partHits.GetHit(iHit)));
    const GGSPosHitsCollection *posHitColl = partHit.GetPosHits();
    auto &tPartHit = *((GGSTPartHitBase *)(tPartHits.ConstructedAt(currTPartHit, "C"))); // Get a cleared hit
    if (posHitColl) {
      // Convert the position hits
      if (!(tPartHit._posHits)) {
        tPartHit._posHits = handle.posHitCAService->ProvideArray(
            posHitColl->GetSize()); // Get a hit array with 0 entries and non-cleared elements
      } else {
        // The particle hit should have given back its position hits array at BeginOfEvent, so this code
        // should never be executed. Throw an exception.
        throw std::logic_error(
            "GGSHitsAction::_Convert: particle hit already has position hits, this should not happen.");
      }
      _Convert(*posHitColl, *(tPartHit._posHits), handle);
    }
    // Save the particle hit if its energy deposit is above threshold or if there are position hits above threshold
    if (partHit.GetEnergyDeposit() > handle.partHitThreshold || (posHitColl && tPartHit._posHits->GetEntries() != 0)) {
      tPartHit.eDep = partHit.GetEnergyDeposit() / GeV;
      tPartHit.UserConversion(partHit);
      // Increase counter for output particle hits
      ++currTPartHit;
    } else {                   // Otherwise remove the particle hit
      if (tPartHit._posHits) { // Before deleting the particle hit we need to give back the array of position hits
        handle.posHitCAService->TakeBackArray(tPartHit._posHits);
      }
      tPartHits.RemoveLast();
    }
  }
}

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

void GGSHitsAction::_Convert(const GGSIntHitsCollection &intHits, TClonesArray &tIntHits, OutputHandle &handle,
                             const std::string &detectorName) {

  if (tIntHits.GetEntries() != 0) {
    throw std::runtime_error("GGSHitsAction::_Convert: dirty TClonesArray for particle hits");
  }
  for (size_t iHit = 0, currTIntHit = 0; iHit < intHits.GetSize(); ++iHit) {
    auto &intHit = *((GGSIntHit *)(intHits.GetHit(iHit)));
    const GGSPartHitsCollection *partHitColl = intHit.GetPartHits();
    auto &tIntHit = *((GGSTIntHitBase *)(tIntHits.ConstructedAt(currTIntHit, "C"))); // Get a cleared hit
    if (partHitColl) {                                                               // Convert the particle hits
      if (!tIntHit._partHits) {
        tIntHit._partHits = handle.partHitCAService->ProvideArray(
            partHitColl->GetSize()); // Get a hit array with 0 entries and non-cleared elements
      } else {
        // The integrated hit should have given back its particle hits array at BeginOfEvent, so this code
        // should never be executed. Throw an exception.
        throw std::logic_error(
            "GGSHitsAction::_Convert: integrated hit already has particle hits, this should not happen.");
      }
      _Convert(*partHitColl, *(tIntHit._partHits), handle);
    }

    // Save the integrated hit if its energy deposit is above threshold or if there are particle hits
    auto totEDep = intHit.GetTotalEnergyDep();
    if (totEDep > handle.intHitThreshold || (partHitColl && tIntHit._partHits->GetEntries() != 0)) {
      tIntHit.eDep = totEDep / GeV;
      if (tIntHit.eDepTimeBin.GetSize() != (int)(intHit.GetEnergyDep().size()))
        tIntHit.eDepTimeBin.Set(intHit.GetEnergyDep().size());
      for (unsigned int i = 0; i < intHit.GetEnergyDep().size(); i++) {
        tIntHit.eDepTimeBin[i] = intHit.GetEnergyDep()[i] / GeV;
      }
      tIntHit.UserConversion(intHit);

      // ----- Save volume informations of this hit ------
      tIntHit._volumeIndex = GGSRootFileService::GetInstance().StoreVolume(_outFile, detectorName, intHit.GetVolume(),
                                                                           intHit.GetAbsTranslation(), intHit.GetID());
      // Increase counter for output integrated hits
      ++currTIntHit;
    } else {                   // Otherwise remove the integrated hit
      if (tIntHit._partHits) { // Before deleting the integrated hit we need to give back the array of particle hits
        handle.partHitCAService->TakeBackArray(tIntHit._partHits);
        // No need to give back arrays for position hits: this code portion is executed only when no particle hits are
        // present, so there are no position hits
      }
      tIntHits.RemoveLast();
    }
  }
}

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

void GGSHitsAction::SetOutputIntHitClass(std::string detectorName, const std::string &className) {
  const std::string routineName("[GGSHitsAction::SetOutputIntHitClass] ");

  // Check if the given class has a dictionary
  TClass *intHitClass = TClass::GetClass(className.c_str());
  if (!intHitClass) {
    throw std::runtime_error(routineName + "Dictionary for class " + className + " not found");
  }

  // Check if the given detector already has a specified output class
  detectorName.insert(detectorName.find_first_of('.'), ".GGSIntHitSD");
  auto handleIter = _outputHandles.find(detectorName);
  if (handleIter == _outputHandles.end()) {
    handleIter = _outputHandles.insert(std::pair<std::string, OutputHandle>{detectorName, OutputHandle()}).first;
  }

  // Set output facilities
  handleIter->second.intHitClassName = className;
}

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

void GGSHitsAction::SetOutputPartHitClass(std::string detectorName, const std::string &className) {
  const std::string routineName("[GGSHitsAction::SetOutputPartHitClass] ");

  // Check if the given class has a dictionary
  TClass *partHitClass = TClass::GetClass(className.c_str());
  if (!partHitClass) {
    throw std::runtime_error(routineName + "Dictionary for class " + className + " not found");
  }

  // Check if the given detector already has a specified output class
  detectorName.insert(detectorName.find_first_of('.'), ".GGSIntHitSD");
  auto handleIter = _outputHandles.find(detectorName);
  if (handleIter == _outputHandles.end()) {
    handleIter = _outputHandles.insert(std::pair<std::string, OutputHandle>{detectorName, OutputHandle()}).first;
  }

  // Set output facilities
  handleIter->second.partHitClassName = className;
}

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

void GGSHitsAction::SetOutputPosHitClass(std::string detectorName, const std::string &className) {
  const std::string routineName("[GGSHitsAction::SetOutputPosHitClass] ");

  // Check if the given class has a dictionary
  TClass *posHitClass = TClass::GetClass(className.c_str());
  if (!posHitClass) {
    throw std::runtime_error(routineName + "Dictionary for class " + className + " not found");
  }

  // Check if the given detector already has a specified output class
  detectorName.insert(detectorName.find_first_of('.'), ".GGSIntHitSD");
  auto handleIter = _outputHandles.find(detectorName);
  if (handleIter == _outputHandles.end()) {
    handleIter = _outputHandles.insert(std::pair<std::string, OutputHandle>{detectorName, OutputHandle()}).first;
  }

  // Set output facilities
  handleIter->second.posHitClassName = className;
}

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

void GGSHitsAction::SetHitThreshold(const std::string &detectorName, const std::string &hitType,
                                    const std::string &valueStr, const std::string &unit) {

  const std::string routineName("GGSHitsAction::SetHitThreshold");

  int hitTypeNo = -1;
  if (hitType == "integrated") {
    hitTypeNo = 0;
  } else if (hitType == "particle") {
    hitTypeNo = 1;
  } else if (hitType == "position") {
    hitTypeNo = 2;
  } else {
    throw std::runtime_error(routineName + ": Unknown hit type " + hitType + " for detector " + detectorName);
  }

  float value = 0.;
  try {
    value = std::stof(valueStr);
  } catch (std::invalid_argument &exc) {
    GGSCOUT(ERROR) << "Invalid threshold for " << hitType << " hits of detector " << detectorName << GGSENDL;
    throw exc;
  } catch (std::out_of_range &exc) {
    GGSCOUT(ERROR) << "Invalid value for threshold of " << hitType << " hits of detector " << detectorName << GGSENDL;
    throw exc;
  }

  auto &unitsTable = G4UnitDefinition::GetUnitsTable();
  auto unitCatIter = std::find_if(unitsTable.begin(), unitsTable.end(),
                                  [](const G4UnitsCategory *unitCat) { return (unitCat->GetName() == "Energy"); });
  if (unitCatIter == unitsTable.end()) {
    throw std::runtime_error(routineName + ": energy unit category not found");
  }

  auto &unitsList = (*unitCatIter)->GetUnitsList();
  auto unitIter = std::find_if(unitsList.begin(), unitsList.end(), [&unit](const G4UnitDefinition *unitDef) {
    return unit == unitDef->GetSymbol().c_str();
  });
  if (unitIter == unitsList.end()) {
    throw std::runtime_error(routineName + ": unit " + unit + " not found");
  }

  auto thresholdIter = std::find_if(
      _thresholds.begin(), _thresholds.end(), [&detectorName, hitTypeNo](const ThresholdStruct &threshold) {
        return (threshold.hitType == hitTypeNo) && (threshold.detectorName == detectorName);
      });
  if (thresholdIter == _thresholds.end()) {
    _thresholds.push_back(ThresholdStruct{detectorName, hitTypeNo, -1.});
    thresholdIter = _thresholds.end() - 1;
  }

  auto *sd =
      dynamic_cast<GGSIntHitSD *>(G4SDManager::GetSDMpointer()->FindSensitiveDetector(detectorName + ".GGSIntHitSD"));
  if (!sd) {
    throw std::invalid_argument(routineName + ": detector " + detectorName + " not found");
  }

  switch (hitTypeNo) {
  case 0:
    thresholdIter->value = value * (*unitIter)->GetValue();
    sd->SetIntHitThreshold(thresholdIter->value);
    break;
  case 1:
    thresholdIter->value = value * (*unitIter)->GetValue();
    sd->SetPartHitThreshold(thresholdIter->value);
    break;
  case 2:
    thresholdIter->value = value * (*unitIter)->GetValue();
    sd->SetPosHitThreshold(thresholdIter->value);
    break;
  }
}

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

GGSHitsAction::GGSHitsActionPrivateMessenger::GGSHitsActionPrivateMessenger(GGSHitsAction *hitsAction)
    : _hitsAction(hitsAction) {
  _setThreshCmd = std::make_unique<G4UIcmdWithAString>("/GGS/userActions/hitsAction/setHitThreshold", this);
  _setThreshCmd->SetGuidance("Set energy deposit threshold for hit output");
  _setThreshCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
}

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

void GGSHitsAction::GGSHitsActionPrivateMessenger::SetNewValue(G4UIcommand *command, G4String newValue) {
  const std::string routineName("GGSHitsAction::GGSHitsActionPrivateMessenger::SetNewValue");
  if (command == _setThreshCmd.get()) {
    auto tokens = GGSStringUtils::Tokenize(newValue);
    if (tokens.size() != 4) {
      throw std::runtime_error(routineName + ": setHitTreshold expects 4 parameters, but " +
                               std::to_string(tokens.size()) + " were given.");
    }

    _hitsAction->SetHitThreshold(tokens[0], tokens[1], tokens[2], tokens[3]);
  }
}

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

void GGSHitsAction::_GiveBackAllArrays() {
  for (auto &handleElem : _outputHandles) {
    OutputHandle &handle = handleElem.second;
    int nIntHits = handle.intHitArray->GetEntries();
    for (int iIntHit = 0; iIntHit < nIntHits; iIntHit++) {
      GGSTIntHitBase *intHit = (GGSTIntHitBase *)(handle.intHitArray->At(iIntHit));
      if (intHit->_partHits) {
        int nPartHits = intHit->GetNPartHits();
        for (int iPartHit = 0; iPartHit < nPartHits; iPartHit++) {
          GGSTPartHitBase *partHit = intHit->GetPartHit(iPartHit);
          if (partHit->_posHits) {
            handle.posHitCAService->TakeBackArray(partHit->_posHits);
          }
        }
        handle.partHitCAService->TakeBackArray(intHit->_partHits);
      }
    }
  }
}