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/services/GGSRootFileService.h"
#include "utils/GGSNameDecoder.h"
#include "utils/GGSSmartLog.h"
#include "montecarlo/scoring/GGSIntHitSD.h"
#include "montecarlo/scoring/manager/GGSMultiSensitiveDetector.h"
#include "montecarlo/pluginmanagers/GGSMCPluginMacros.h"
#include "utils/GGSStringUtils.h"

#include "G4SystemOfUnits.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(
    NULL), _outTree(NULL), _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) {
  static const std::string routineName("GSHitsAction::BeginOfRunAction");

  // -----------------
  // Open the run file
  // -----------------
  _outFile = GGSRootFileService::GetInstance().GetFileForThisRun(_outFileBase, run);
  if (!_outFile || (_outFile && _outFile->IsZombie())) {
    COUT(ERROR) << "Cannot open ROOT file " << _outFileBase << "for run" << run->GetRunID() << ENDL;
    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;

  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) {
      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
          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

    } // End loop on sensitive detectors

  } // End loop on logical volumes

}

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

void GGSHitsAction::EndOfRunAction(const G4Run *) {

// -----------------
// 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); // Will automatically save detector informations
    _outFile = NULL;
    _outTree = NULL; // 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 *) {

  _GiveBackAllArrays();

  // Empty the array and call Clear(Option_t*) of the objects
  // This will clear also the eventual position and particle hit arrays
  for (auto &handleElem : _outputHandles) {
    OutputHandle &handle = handleElem.second;
    handle.intHitArray->Clear("C");
    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
    _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) {
      GGSTIntHit *tIntHit = (GGSTIntHit*) (handle.intHitArray->At(iIntHit));
      if (tIntHit->_partHits) {
        handle.partHitArray->AddAtAndExpand(tIntHit->_partHits, handle.partHitArray->GetEntries());
        if (tIntHit->_partHits->GetEntries() > 0 && ((GGSTPartHit*) (tIntHit->_partHits->At(0)))->_posHits) {
          static TClonesArray *partHits = NULL;
          static Int_t nPartHits = 0;
          partHits = tIntHit->_partHits;
          nPartHits = partHits->GetEntries();
          for (Int_t iPartHit = 0; iPartHit < nPartHits; iPartHit++) {
            GGSTPartHit *partHit = (GGSTPartHit*) (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) {

  for (size_t iHit = 0; iHit < posHits.GetSize(); ++iHit) {
    auto &posHit = *((GGSPosHit*) (posHits.GetHit(iHit)));
    if (posHit.GetEnergyDeposit() > handle.posHitThreshold) {
      GGSTPosHit &tPosHit = *((GGSTPosHit*) (tPosHits.ConstructedAt(tPosHits.GetEntries(), "")));
      tPosHit.eDep = posHit.GetEnergyDeposit() / GeV;
      tPosHit.time = posHit.GetTime() / ns;
      tPosHit.pathLength = posHit.GetPathLength() / cm;
      tPosHit.startPoint[0] = posHit.GetStartPoint()[0] / cm;
      tPosHit.startPoint[1] = posHit.GetStartPoint()[1] / cm;
      tPosHit.startPoint[2] = posHit.GetStartPoint()[2] / cm;
      tPosHit.endPoint[0] = posHit.GetEndPoint()[0] / cm;
      tPosHit.endPoint[1] = posHit.GetEndPoint()[1] / cm;
      tPosHit.endPoint[2] = posHit.GetEndPoint()[2] / cm;
      tPosHit.startMomentum[0] = posHit.GetStartMomentum()[0] / GeV;
      tPosHit.startMomentum[1] = posHit.GetStartMomentum()[1] / GeV;
      tPosHit.startMomentum[2] = posHit.GetStartMomentum()[2] / GeV;
      tPosHit.startEnergy = posHit.GetStartEnergy() / GeV;
      tPosHit.UserConversion(posHit);
    }
  }
}

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

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

  for (size_t iHit = 0; iHit < partHits.GetSize(); ++iHit) {
    auto &partHit = *((GGSPartHit*) (partHits.GetHit(iHit)));
    const GGSPosHitsCollection *posHitColl = partHit.GetPosHits();
    auto &tPartHit = *((GGSTPartHit*) (tPartHits.ConstructedAt(tPartHits.GetEntries(), "")));
    if (posHitColl) {
      // Convert the position hits
      if ((tPartHit._posHits)) {
        tPartHit._posHits->Clear("C");
      }
      else {
        tPartHit._posHits = handle.posHitCAService->ProvideArray(posHitColl->GetSize());
      }
      _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.time = partHit.GetTime() / ns;
      tPartHit.pathLength = partHit.GetPathLength() / cm;
      tPartHit.entrancePoint[0] = partHit.GetEntrancePoint()[0] / cm;
      tPartHit.entrancePoint[1] = partHit.GetEntrancePoint()[1] / cm;
      tPartHit.entrancePoint[2] = partHit.GetEntrancePoint()[2] / cm;
      tPartHit.exitPoint[0] = partHit.GetExitPoint()[0] / cm;
      tPartHit.exitPoint[1] = partHit.GetExitPoint()[1] / cm;
      tPartHit.exitPoint[2] = partHit.GetExitPoint()[2] / cm;
      tPartHit.entranceMomentum[0] = partHit.GetEntranceMomentum()[0] / GeV;
      tPartHit.entranceMomentum[1] = partHit.GetEntranceMomentum()[1] / GeV;
      tPartHit.entranceMomentum[2] = partHit.GetEntranceMomentum()[2] / GeV;
      tPartHit.entranceEnergy = partHit.GetEntranceEnergy() / GeV;
      tPartHit.trackID = partHit.GetTrackID();
      tPartHit.parentID = partHit.GetParentID();
      tPartHit.particlePdg = partHit.GetParticlePdg();
      tPartHit.UserConversion(partHit);
    }
    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) {

  for (size_t iHit = 0; iHit < intHits.GetSize(); ++iHit) {
    auto &intHit = *((GGSIntHit*) (intHits.GetHit(iHit)));
    const GGSPartHitsCollection *partHitColl = intHit.GetPartHits();
    auto &tIntHit = *((GGSTIntHit*) (tIntHits.ConstructedAt(tIntHits.GetEntries(), "")));
    if (partHitColl) { // Convert the particle hits
      if (tIntHit._partHits) {
        tIntHit._partHits->Clear("C");
      }
      else {
        tIntHit._partHits = handle.partHitCAService->ProvideArray(partHitColl->GetSize());
      }
      _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.time = intHit.GetTime() / ns;
      tIntHit.UserConversion(intHit);

      // ----- Save volume informations of this hit ------
      tIntHit._volumeIndex = GGSRootFileService::GetInstance().StoreVolume(_outFile, detectorName, intHit.GetVolume(),
          intHit.GetAbsTranslation(), intHit.GetID());
    }
    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) {
    COUT(ERROR) << "Invalid threshold for " << hitType << " hits of detector " << detectorName << ENDL;
    throw exc;
  } catch (std::out_of_range &exc) {
    COUT(ERROR) << "Invalid value for threshold of " << hitType << " hits of detector " << detectorName << ENDL;
    throw exc;
  }

  auto &unitsTable = G4UnitDefinition::GetUnitsTable();
  auto unitCatIter = std::find_if(unitsTable.begin(), unitsTable.end(),
      [&unit](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;
  }
  switch (hitTypeNo) {
  case 0:
    thresholdIter->value = value * (*unitIter)->GetValue();
    break;
  case 1:
    thresholdIter->value = value * (*unitIter)->GetValue();
    break;
  case 2:
    thresholdIter->value = value * (*unitIter)->GetValue();
    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++) {
      GGSTIntHit *intHit = (GGSTIntHit*) (handle.intHitArray->At(iIntHit));
      if (intHit->_partHits) {
        int nPartHits = intHit->GetNPartHits();
        for (int iPartHit = 0; iPartHit < nPartHits; iPartHit++) {
          GGSTPartHit *partHit = intHit->GetPartHit(iPartHit);
          if (partHit->_posHits) {
            handle.posHitCAService->TakeBackArray(partHit->_posHits);
          }
        }
        handle.partHitCAService->TakeBackArray(intHit->_partHits);
      }
    }
  }
}