G4VUserMoleculeReactionCounter.hh

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// Author: Christian Velten (2025)
 
#ifndef G4VUSERMOLECULEREACTIONCOUNTER_HH
#define G4VUSERMOLECULEREACTIONCOUNTER_HH 1
 
#include "G4DNAChemistryManager.hh"
#include "G4MoleculeCounterTemplates.hh"
#include "G4Scheduler.hh"
#include "G4UnitsTable.hh"
#include "G4VMoleculeReactionCounter.hh"
 
//------------------------------------------------------------------------------
 
template<class TIndex>
class G4VUserMoleculeReactionCounter : public G4VMoleculeReactionCounter
{
    static_assert(std::is_base_of<G4VMoleculeReactionCounter::G4VMoleculeReactionCounterIndex, TIndex>::value,
      "TIndex must be derived from G4VMoleculeReactionCounter::G4VMoleculeReactionCounterIndex! "
      "No forward declaration is allowed.");
 
  protected:
    struct Search;
 
  public:
    G4VUserMoleculeReactionCounter();
    G4VUserMoleculeReactionCounter(const G4String&,
                                   MoleculeReactionCounterType = MoleculeReactionCounterType::Basic);
    ~G4VUserMoleculeReactionCounter() override = default;
 
  public:
    void Initialize() final;
    void InitializeUser() override = 0;
    void ResetCounter() override;
    void Dump() const override;
    void DumpCounterMapIndices() const override;
 
    void AbsorbCounter(const G4VMoleculeCounterInternalBase*) override;
 
    std::unique_ptr<G4VMoleculeReactionCounterIndex> BuildSimpleIndex(const G4DNAMolecularReactionData*) const override = 0;
 
    void RecordReaction(std::unique_ptr<G4VMoleculeReactionCounterIndex>, G4double, G4int = 1) override;
 
    std::set<const G4DNAMolecularReactionData*> GetRecordedReactions() const override;
    std::set<G4double> GetRecordedTimes() const override;
 
  protected:
    std::map<TIndex, InnerCounterMapType> fCounterMap{};
 
  public:
    const std::map<TIndex, InnerCounterMapType>& GetCounterMap() const { return fCounterMap; }
    std::vector<TIndex> GetMapIndices() const;
 
    virtual G4int GetNbReactionsAtTime(const TIndex&, G4double) const;
    virtual G4int GetNbReactionsAtTime(Search&, const TIndex&, G4double) const;
    virtual std::vector<G4int> GetNbReactionsAtTimes(const TIndex&,
                                                     const std::vector<G4double>&) const;
 
    //-SEARCH-----------------------------------------------------------------------
  protected:
    struct Search
    {
        Search() : fLowerBoundSet(false) {}
        typename std::map<TIndex, InnerCounterMapType>::const_iterator fLastIndexSearched;
        InnerCounterMapType::const_iterator fLowerBoundTime;
        G4bool fLowerBoundSet;
    };
    G4bool SearchIndexUpdated(Search&, const TIndex&) const;
    G4int SearchUpperBoundTime(Search&, G4double, G4bool) const;
};
 
//------------------------------------------------------------------------------
 
// #include "G4VUserMoleculeReactionCounter.icc"
 
//------------------------------------------------------------------------------
 
template<typename T>
G4VUserMoleculeReactionCounter<T>::G4VUserMoleculeReactionCounter() : G4VMoleculeReactionCounter()
{}
 
//------------------------------------------------------------------------------
 
template<typename T>
G4VUserMoleculeReactionCounter<T>::G4VUserMoleculeReactionCounter(const G4String& name,
                                                                  MoleculeReactionCounterType type)
  : G4VMoleculeReactionCounter(name, type)
{}
 
//------------------------------------------------------------------------------
 
template<typename TIndex>
void G4VUserMoleculeReactionCounter<TIndex>::Initialize()
{
  InitializeUser();
  fIsInitialized = true;
}
 
//------------------------------------------------------------------------------
 
template<typename TIndex>
G4int G4VUserMoleculeReactionCounter<TIndex>::GetNbReactionsAtTime(const TIndex& index, G4double time) const
{
  Search search = {};
  return GetNbReactionsAtTime(search, index, time);
}
 
//------------------------------------------------------------------------------
 
template<typename TIndex>
G4int G4VUserMoleculeReactionCounter<TIndex>::GetNbReactionsAtTime(Search& search,
                                                                   const TIndex& index,
                                                                   G4double time) const
{
  G4bool sameIndex = !SearchIndexUpdated(search, index);
  return SearchUpperBoundTime(search, time, sameIndex);
}
 
//------------------------------------------------------------------------------
 
template<typename TIndex>
std::vector<G4int> G4VUserMoleculeReactionCounter<TIndex>::GetNbReactionsAtTimes(
  const TIndex& index, const std::vector<G4double>& times) const
{
  Search search = {};
  std::vector<G4int> counts = {};
  for (auto time : times)
    counts.push_back(GetNbReactionsAtTime(search, index, time));
  return counts;
}
 
//------------------------------------------------------------------------------
 
template<typename TIndex>
void G4VUserMoleculeReactionCounter<TIndex>::RecordReaction(
  std::unique_ptr<G4VMoleculeReactionCounter::G4VMoleculeReactionCounterIndex> pIndex,
  G4double time, G4int number)
{
  const TIndex* mapIndex = dynamic_cast<TIndex*>(pIndex.get());
 
  if (IsTimeAboveUpperBound(time)) {
    if (fVerbose > 3) {
      G4cout << "G4VUserMoleculeReactionCounter<"
             << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">(" << GetName()
             << ")::RecordReaction : " << mapIndex->GetReactionData()->GetReactionID()
             << " at time : " << G4BestUnit(time, "Time") << G4endl;
      G4cout << ":: [IsTimeAboveUpperBound] Skipping since IsTimeAboveUpperBound == true" << G4endl;
    }
    return;
  }
  else if (IsTimeBelowLowerBound(time)) {
    if (fVerbose > 3) {
      G4cout << "G4VUserMoleculeReactionCounter<"
             << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">(" << GetName()
             << ")::RecordReaction : " << mapIndex->GetReactionData()->GetReactionID()
             << " at time : " << G4BestUnit(time, "Time") << G4endl;
      G4cout << ":: [IsTimeBelowLowerBound] Skipping since IsTimeBelowLowerBound == true" << G4endl;
    }
    return;
  }
 
  if (fVerbose > 2) {
    G4cout << "G4VUserMoleculeReactionCounter<"
           << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">(" << GetName()
           << ")::RecordReaction : " << mapIndex->GetReactionData()->GetReactionID()
           << " at time : " << G4BestUnit(time, "Time") << G4endl;
  }
 
  auto [it, indexIsNew] = fCounterMap.emplace(*mapIndex, InnerCounterMapType{fTimeComparer});
 
  if (indexIsNew || it->second.empty()) {
    it->second.emplace(time, number);
    // it->second[time] = number;
  }
  else {
    if (G4MoleculeCounterManager::Instance()->GetResetCountersBeforeEvent())
    // can only do consistency check if the counters are cleared before each event
    {
      auto end = it->second.rbegin();
 
      if ((end->first <= time || std::fabs(end->first - time) <= fTimeComparer.GetPrecisionAtTime(time)))
      // Case 1 = new time comes after last recorded data
      // Case 2 = new time is about the same as the last recorded one
      {
        // it->second[time] = end->second + number;
        auto [it_time, _] = it->second.emplace(time, end->second);
        it_time->second += number;
      }
      else {
        G4ExceptionDescription errMsg;
        errMsg << "Time of reaction " << mapIndex->GetReactionData()->GetReactionID() << " is "
               << G4BestUnit(time, "Time") << "while the global time is "
               << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
               << "(last counter time: " << G4BestUnit(end->first, "Time") << ")" << G4endl;
        G4Exception(G4String("G4VUserMoleculeReactionCounter<"
                             + G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
                             + ">::RecordReaction"),
                    "TIME_DONT_MATCH", FatalException, errMsg);
      }
    }
    else {
      // since counters are not cleared after chemical run (i.e., after event)
      // there will already be numbers in the map, so...
      // (1) find the closest time
      // (2) emplace entry using closest value as init + number
      // (3) add number to all "future" entries as well
      auto it_closest = G4::MoleculeCounter::FindClosestEntryForKey(it->second, time);
      auto [it_new, _] = it->second.emplace(time, it_closest->second);
      do {
        it_new->second += number;
      } while (++it_new != it->second.end());
    }
  }
}
 
//------------------------------------------------------------------------------
 
template<typename TIndex>
std::vector<TIndex> G4VUserMoleculeReactionCounter<TIndex>::GetMapIndices() const
{
  if (fVerbose > 2) {
    G4cout << "Entering in G4VUserMoleculeReactionCounter::GetMapIndices" << G4endl;
  }
  return G4::MoleculeCounter::GetMapIndices(fCounterMap);
}
 
//------------------------------------------------------------------------------
 
template<typename T>
std::set<const G4DNAMolecularReactionData*> G4VUserMoleculeReactionCounter<T>::GetRecordedReactions() const
{
  if (fVerbose > 2) {
    G4cout << "Entering in G4VUserMoleculeReactionCounter::GetRecordedReactions" << G4endl;
  }
  std::set<const G4DNAMolecularReactionData*> output{};
  for (const auto& it : fCounterMap) {
    output.insert(it.first.GetReactionData());
  }
  return output;
}
 
//------------------------------------------------------------------------------
 
template<typename T>
std::set<G4double> G4VUserMoleculeReactionCounter<T>::GetRecordedTimes() const
{
  return G4::MoleculeCounter::GetRecordedTimes<T>(fCounterMap);
}
 
//------------------------------------------------------------------------------
 
template<typename T>
void G4VUserMoleculeReactionCounter<T>::Dump() const
{
  DumpCounterMapIndices();
  G4::MoleculeCounter::DumpCounterMapContents<T>(fCounterMap);
}
 
template<typename T>
void G4VUserMoleculeReactionCounter<T>::DumpCounterMapIndices() const
{
  G4::MoleculeCounter::DumpCounterMapIndices<T>(fCounterMap);
}
 
//------------------------------------------------------------------------------
 
template<typename T>
void G4VUserMoleculeReactionCounter<T>::ResetCounter()
{
  if (fVerbose > 1) {
    G4cout << "G4VUserMoleculeReactionCounter<" << G4::MoleculeCounter::GetTemplateTypeName<T>()
           << ">(" << GetName() << ")::ResetCounter" << G4endl;
  }
  fCounterMap.clear();
}
 
//------------------------------------------------------------------------------
 
template<typename TIndex>
G4bool G4VUserMoleculeReactionCounter<TIndex>::SearchIndexUpdated(Search& search, const TIndex& index) const
{
  if (search.fLowerBoundSet && !(search.fLastIndexSearched->first < index)
      && !(index < search.fLastIndexSearched->first))
  {
    return true;
  }
 
  auto mol_it = fCounterMap.find(index);
  search.fLastIndexSearched = mol_it;
 
  if (mol_it != fCounterMap.end()) {
    search.fLowerBoundTime = search.fLastIndexSearched->second.end();
    search.fLowerBoundSet = true;
  }
  else {
    search.fLowerBoundSet = false;
  }
 
  return false;
}
 
//------------------------------------------------------------------------------
 
template<typename T>
G4int G4VUserMoleculeReactionCounter<T>::SearchUpperBoundTime(Search& search, G4double time, G4bool sameIndex) const
{
  auto mol_it = search.fLastIndexSearched;
  if (mol_it == fCounterMap.end()) {
    return 0;
  }
 
  InnerCounterMapType const& timeMap = mol_it->second;
  if (timeMap.empty()) {
    return 0;
  }
 
  if (sameIndex) {
    if (search.fLowerBoundSet && search.fLowerBoundTime != timeMap.end()) {
      if (search.fLowerBoundTime->first < time) {
        auto upperToLast = search.fLowerBoundTime;
        upperToLast++;
 
        if (upperToLast == timeMap.end()) {
          return search.fLowerBoundTime->second;
        }
 
        if (upperToLast->first > time) {
          return search.fLowerBoundTime->second;
        }
      }
    }
  }
 
  auto up_time_it = timeMap.upper_bound(time);
 
  if (up_time_it == timeMap.end()) {
    auto last_time = timeMap.rbegin();
    return last_time->second;
  }
  if (up_time_it == timeMap.begin()) {
    return 0;
  }
 
  up_time_it--;
 
  search.fLowerBoundTime = up_time_it;
  search.fLowerBoundSet = true;
 
  return search.fLowerBoundTime->second;
}
 
//------------------------------------------------------------------------------
 
template<typename TIndex>
void G4VUserMoleculeReactionCounter<TIndex>::AbsorbCounter(
  const G4VMoleculeCounterInternalBase* pCounterBase)
{
  if (pCounterBase == nullptr) {
    G4ExceptionDescription errMsg;
    errMsg << "Could not cast the pointer to type G4VUserMoleculeReactionCounter<"
           << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">!\n"
           << "Because the pointer is nullptr!" << G4endl;
    G4Exception(G4String("G4VUserMoleculeReactionCounter<"
                         + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AbsorbCounter"),
                "BAD_REFERENCE", FatalException, errMsg);
  }
 
  auto pCounter = dynamic_cast<G4VUserMoleculeReactionCounter<TIndex> const*>(pCounterBase);
 
  if (pCounter == nullptr) {
    G4ExceptionDescription errMsg;
    errMsg << "Could not cast the pointer to type G4VUserMoleculeReactionCounter<"
           << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">!\n"
           << "Because the objects aren't of the same type!" << G4endl;
    G4Exception(G4String("G4VUserMoleculeReactionCounter<"
                         + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AbsorbCounter"),
                "BAD_REFERENCE", FatalException, errMsg);
  }
 
  if (pCounter->GetType() != GetType()) {
    G4ExceptionDescription errMsg;
    errMsg << "You are trying to absorb a counter with different type!" << G4endl;
    G4Exception(G4String("G4VUserMoleculeReactionCounter<"
                         + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AbsorbCounter"),
                "TYPE_DIFF", JustWarning, errMsg);
  }
 
  for (auto const& worker_it : pCounter->GetCounterMap()) {
    auto [master_it, indexIsNew] =
      fCounterMap.emplace(worker_it.first, InnerCounterMapType{fTimeComparer});
 
    G4int currentNumber = 0, previousNumber = 0;
    for (auto const& [time, number] : worker_it.second) {
      currentNumber = number - previousNumber;
      previousNumber = number;
 
      if (master_it->second.empty()) {
        master_it->second.emplace(time, currentNumber);
      }
      else {  // at least one element exists, so we can try to find the closest key
        auto it_closest = G4::MoleculeCounter::FindClosestEntryForKey(master_it->second, time);
        auto [it, _] = master_it->second.emplace(time, it_closest->second);
        do {
          it->second += currentNumber;
        } while (++it != master_it->second.end());
      }
    }
  }
}
 
//------------------------------------------------------------------------------
 
#endif