include/Geant4/G4VUserMoleculeReactionCounter.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Author: Christian Velten (2025)
0027
0028 #ifndef G4VUSERMOLECULEREACTIONCOUNTER_HH
0029 #define G4VUSERMOLECULEREACTIONCOUNTER_HH 1
0030
0031 #include "G4DNAChemistryManager.hh"
0032 #include "G4MoleculeCounterTemplates.hh"
0033 #include "G4Scheduler.hh"
0034 #include "G4UnitsTable.hh"
0035 #include "G4VMoleculeReactionCounter.hh"
0036
0037 //------------------------------------------------------------------------------
0038
0039 template<class TIndex>
0040 class G4VUserMoleculeReactionCounter : public G4VMoleculeReactionCounter
0041 {
0042     static_assert(std::is_base_of<G4VMoleculeReactionCounter::G4VMoleculeReactionCounterIndex, TIndex>::value,
0043       "TIndex must be derived from G4VMoleculeReactionCounter::G4VMoleculeReactionCounterIndex! "
0044       "No forward declaration is allowed.");
0045
0046   protected:
0047     struct Search;
0048
0049   public:
0050     G4VUserMoleculeReactionCounter();
0051     G4VUserMoleculeReactionCounter(const G4String&,
0052                                    MoleculeReactionCounterType = MoleculeReactionCounterType::Basic);
0053     ~G4VUserMoleculeReactionCounter() override = default;
0054
0055   public:
0056     void Initialize() final;
0057     void InitializeUser() override = 0;
0058     void ResetCounter() override;
0059     void Dump() const override;
0060     void DumpCounterMapIndices() const override;
0061
0062     void AbsorbCounter(const G4VMoleculeCounterInternalBase*) override;
0063
0064     std::unique_ptr<G4VMoleculeReactionCounterIndex> BuildSimpleIndex(const G4DNAMolecularReactionData*) const override = 0;
0065
0066     void RecordReaction(std::unique_ptr<G4VMoleculeReactionCounterIndex>, G4double, G4int = 1) override;
0067
0068     std::set<const G4DNAMolecularReactionData*> GetRecordedReactions() const override;
0069     std::set<G4double> GetRecordedTimes() const override;
0070
0071   protected:
0072     std::map<TIndex, InnerCounterMapType> fCounterMap{};
0073
0074   public:
0075     const std::map<TIndex, InnerCounterMapType>& GetCounterMap() const { return fCounterMap; }
0076     std::vector<TIndex> GetMapIndices() const;
0077
0078     virtual G4int GetNbReactionsAtTime(const TIndex&, G4double) const;
0079     virtual G4int GetNbReactionsAtTime(Search&, const TIndex&, G4double) const;
0080     virtual std::vector<G4int> GetNbReactionsAtTimes(const TIndex&,
0081                                                      const std::vector<G4double>&) const;
0082
0083     //-SEARCH-----------------------------------------------------------------------
0084   protected:
0085     struct Search
0086     {
0087         Search() : fLowerBoundSet(false) {}
0088         typename std::map<TIndex, InnerCounterMapType>::const_iterator fLastIndexSearched;
0089         InnerCounterMapType::const_iterator fLowerBoundTime;
0090         G4bool fLowerBoundSet;
0091     };
0092     G4bool SearchIndexUpdated(Search&, const TIndex&) const;
0093     G4int SearchUpperBoundTime(Search&, G4double, G4bool) const;
0094 };
0095
0096 //------------------------------------------------------------------------------
0097
0098 // #include "G4VUserMoleculeReactionCounter.icc"
0099
0100 //------------------------------------------------------------------------------
0101
0102 template<typename T>
0103 G4VUserMoleculeReactionCounter<T>::G4VUserMoleculeReactionCounter() : G4VMoleculeReactionCounter()
0104 {}
0105
0106 //------------------------------------------------------------------------------
0107
0108 template<typename T>
0109 G4VUserMoleculeReactionCounter<T>::G4VUserMoleculeReactionCounter(const G4String& name,
0110                                                                   MoleculeReactionCounterType type)
0111   : G4VMoleculeReactionCounter(name, type)
0112 {}
0113
0114 //------------------------------------------------------------------------------
0115
0116 template<typename TIndex>
0117 void G4VUserMoleculeReactionCounter<TIndex>::Initialize()
0118 {
0119   InitializeUser();
0120   fIsInitialized = true;
0121 }
0122
0123 //------------------------------------------------------------------------------
0124
0125 template<typename TIndex>
0126 G4int G4VUserMoleculeReactionCounter<TIndex>::GetNbReactionsAtTime(const TIndex& index, G4double time) const
0127 {
0128   Search search = {};
0129   return GetNbReactionsAtTime(search, index, time);
0130 }
0131
0132 //------------------------------------------------------------------------------
0133
0134 template<typename TIndex>
0135 G4int G4VUserMoleculeReactionCounter<TIndex>::GetNbReactionsAtTime(Search& search,
0136                                                                    const TIndex& index,
0137                                                                    G4double time) const
0138 {
0139   G4bool sameIndex = !SearchIndexUpdated(search, index);
0140   return SearchUpperBoundTime(search, time, sameIndex);
0141 }
0142
0143 //------------------------------------------------------------------------------
0144
0145 template<typename TIndex>
0146 std::vector<G4int> G4VUserMoleculeReactionCounter<TIndex>::GetNbReactionsAtTimes(
0147   const TIndex& index, const std::vector<G4double>& times) const
0148 {
0149   Search search = {};
0150   std::vector<G4int> counts = {};
0151   for (auto time : times)
0152     counts.push_back(GetNbReactionsAtTime(search, index, time));
0153   return counts;
0154 }
0155
0156 //------------------------------------------------------------------------------
0157
0158 template<typename TIndex>
0159 void G4VUserMoleculeReactionCounter<TIndex>::RecordReaction(
0160   std::unique_ptr<G4VMoleculeReactionCounter::G4VMoleculeReactionCounterIndex> pIndex,
0161   G4double time, G4int number)
0162 {
0163   const TIndex* mapIndex = dynamic_cast<TIndex*>(pIndex.get());
0164
0165   if (IsTimeAboveUpperBound(time)) {
0166     if (fVerbose > 3) {
0167       G4cout << "G4VUserMoleculeReactionCounter<"
0168              << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">(" << GetName()
0169              << ")::RecordReaction : " << mapIndex->GetReactionData()->GetReactionID()
0170              << " at time : " << G4BestUnit(time, "Time") << G4endl;
0171       G4cout << ":: [IsTimeAboveUpperBound] Skipping since IsTimeAboveUpperBound == true" << G4endl;
0172     }
0173     return;
0174   }
0175   else if (IsTimeBelowLowerBound(time)) {
0176     if (fVerbose > 3) {
0177       G4cout << "G4VUserMoleculeReactionCounter<"
0178              << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">(" << GetName()
0179              << ")::RecordReaction : " << mapIndex->GetReactionData()->GetReactionID()
0180              << " at time : " << G4BestUnit(time, "Time") << G4endl;
0181       G4cout << ":: [IsTimeBelowLowerBound] Skipping since IsTimeBelowLowerBound == true" << G4endl;
0182     }
0183     return;
0184   }
0185
0186   if (fVerbose > 2) {
0187     G4cout << "G4VUserMoleculeReactionCounter<"
0188            << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">(" << GetName()
0189            << ")::RecordReaction : " << mapIndex->GetReactionData()->GetReactionID()
0190            << " at time : " << G4BestUnit(time, "Time") << G4endl;
0191   }
0192
0193   auto [it, indexIsNew] = fCounterMap.emplace(*mapIndex, InnerCounterMapType{fTimeComparer});
0194
0195   if (indexIsNew || it->second.empty()) {
0196     it->second.emplace(time, number);
0197     // it->second[time] = number;
0198   }
0199   else {
0200     if (G4MoleculeCounterManager::Instance()->GetResetCountersBeforeEvent())
0201     // can only do consistency check if the counters are cleared before each event
0202     {
0203       auto end = it->second.rbegin();
0204
0205       if ((end->first <= time || std::fabs(end->first - time) <= fTimeComparer.GetPrecisionAtTime(time)))
0206       // Case 1 = new time comes after last recorded data
0207       // Case 2 = new time is about the same as the last recorded one
0208       {
0209         // it->second[time] = end->second + number;
0210         auto [it_time, _] = it->second.emplace(time, end->second);
0211         it_time->second += number;
0212       }
0213       else {
0214         G4ExceptionDescription errMsg;
0215         errMsg << "Time of reaction " << mapIndex->GetReactionData()->GetReactionID() << " is "
0216                << G4BestUnit(time, "Time") << "while the global time is "
0217                << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
0218                << "(last counter time: " << G4BestUnit(end->first, "Time") << ")" << G4endl;
0219         G4Exception(G4String("G4VUserMoleculeReactionCounter<"
0220                              + G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0221                              + ">::RecordReaction"),
0222                     "TIME_DONT_MATCH", FatalException, errMsg);
0223       }
0224     }
0225     else {
0226       // since counters are not cleared after chemical run (i.e., after event)
0227       // there will already be numbers in the map, so...
0228       // (1) find the closest time
0229       // (2) emplace entry using closest value as init + number
0230       // (3) add number to all "future" entries as well
0231       auto it_closest = G4::MoleculeCounter::FindClosestEntryForKey(it->second, time);
0232       auto [it_new, _] = it->second.emplace(time, it_closest->second);
0233       do {
0234         it_new->second += number;
0235       } while (++it_new != it->second.end());
0236     }
0237   }
0238 }
0239
0240 //------------------------------------------------------------------------------
0241
0242 template<typename TIndex>
0243 std::vector<TIndex> G4VUserMoleculeReactionCounter<TIndex>::GetMapIndices() const
0244 {
0245   if (fVerbose > 2) {
0246     G4cout << "Entering in G4VUserMoleculeReactionCounter::GetMapIndices" << G4endl;
0247   }
0248   return G4::MoleculeCounter::GetMapIndices(fCounterMap);
0249 }
0250
0251 //------------------------------------------------------------------------------
0252
0253 template<typename T>
0254 std::set<const G4DNAMolecularReactionData*> G4VUserMoleculeReactionCounter<T>::GetRecordedReactions() const
0255 {
0256   if (fVerbose > 2) {
0257     G4cout << "Entering in G4VUserMoleculeReactionCounter::GetRecordedReactions" << G4endl;
0258   }
0259   std::set<const G4DNAMolecularReactionData*> output{};
0260   for (const auto& it : fCounterMap) {
0261     output.insert(it.first.GetReactionData());
0262   }
0263   return output;
0264 }
0265
0266 //------------------------------------------------------------------------------
0267
0268 template<typename T>
0269 std::set<G4double> G4VUserMoleculeReactionCounter<T>::GetRecordedTimes() const
0270 {
0271   return G4::MoleculeCounter::GetRecordedTimes<T>(fCounterMap);
0272 }
0273
0274 //------------------------------------------------------------------------------
0275
0276 template<typename T>
0277 void G4VUserMoleculeReactionCounter<T>::Dump() const
0278 {
0279   DumpCounterMapIndices();
0280   G4::MoleculeCounter::DumpCounterMapContents<T>(fCounterMap);
0281 }
0282
0283 template<typename T>
0284 void G4VUserMoleculeReactionCounter<T>::DumpCounterMapIndices() const
0285 {
0286   G4::MoleculeCounter::DumpCounterMapIndices<T>(fCounterMap);
0287 }
0288
0289 //------------------------------------------------------------------------------
0290
0291 template<typename T>
0292 void G4VUserMoleculeReactionCounter<T>::ResetCounter()
0293 {
0294   if (fVerbose > 1) {
0295     G4cout << "G4VUserMoleculeReactionCounter<" << G4::MoleculeCounter::GetTemplateTypeName<T>()
0296            << ">(" << GetName() << ")::ResetCounter" << G4endl;
0297   }
0298   fCounterMap.clear();
0299 }
0300
0301 //------------------------------------------------------------------------------
0302
0303 template<typename TIndex>
0304 G4bool G4VUserMoleculeReactionCounter<TIndex>::SearchIndexUpdated(Search& search, const TIndex& index) const
0305 {
0306   if (search.fLowerBoundSet && !(search.fLastIndexSearched->first < index)
0307       && !(index < search.fLastIndexSearched->first))
0308   {
0309     return true;
0310   }
0311
0312   auto mol_it = fCounterMap.find(index);
0313   search.fLastIndexSearched = mol_it;
0314
0315   if (mol_it != fCounterMap.end()) {
0316     search.fLowerBoundTime = search.fLastIndexSearched->second.end();
0317     search.fLowerBoundSet = true;
0318   }
0319   else {
0320     search.fLowerBoundSet = false;
0321   }
0322
0323   return false;
0324 }
0325
0326 //------------------------------------------------------------------------------
0327
0328 template<typename T>
0329 G4int G4VUserMoleculeReactionCounter<T>::SearchUpperBoundTime(Search& search, G4double time, G4bool sameIndex) const
0330 {
0331   auto mol_it = search.fLastIndexSearched;
0332   if (mol_it == fCounterMap.end()) {
0333     return 0;
0334   }
0335
0336   InnerCounterMapType const& timeMap = mol_it->second;
0337   if (timeMap.empty()) {
0338     return 0;
0339   }
0340
0341   if (sameIndex) {
0342     if (search.fLowerBoundSet && search.fLowerBoundTime != timeMap.end()) {
0343       if (search.fLowerBoundTime->first < time) {
0344         auto upperToLast = search.fLowerBoundTime;
0345         upperToLast++;
0346
0347         if (upperToLast == timeMap.end()) {
0348           return search.fLowerBoundTime->second;
0349         }
0350
0351         if (upperToLast->first > time) {
0352           return search.fLowerBoundTime->second;
0353         }
0354       }
0355     }
0356   }
0357
0358   auto up_time_it = timeMap.upper_bound(time);
0359
0360   if (up_time_it == timeMap.end()) {
0361     auto last_time = timeMap.rbegin();
0362     return last_time->second;
0363   }
0364   if (up_time_it == timeMap.begin()) {
0365     return 0;
0366   }
0367
0368   up_time_it--;
0369
0370   search.fLowerBoundTime = up_time_it;
0371   search.fLowerBoundSet = true;
0372
0373   return search.fLowerBoundTime->second;
0374 }
0375
0376 //------------------------------------------------------------------------------
0377
0378 template<typename TIndex>
0379 void G4VUserMoleculeReactionCounter<TIndex>::AbsorbCounter(
0380   const G4VMoleculeCounterInternalBase* pCounterBase)
0381 {
0382   if (pCounterBase == nullptr) {
0383     G4ExceptionDescription errMsg;
0384     errMsg << "Could not cast the pointer to type G4VUserMoleculeReactionCounter<"
0385            << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">!\n"
0386            << "Because the pointer is nullptr!" << G4endl;
0387     G4Exception(G4String("G4VUserMoleculeReactionCounter<"
0388                          + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AbsorbCounter"),
0389                 "BAD_REFERENCE", FatalException, errMsg);
0390   }
0391
0392   auto pCounter = dynamic_cast<G4VUserMoleculeReactionCounter<TIndex> const*>(pCounterBase);
0393
0394   if (pCounter == nullptr) {
0395     G4ExceptionDescription errMsg;
0396     errMsg << "Could not cast the pointer to type G4VUserMoleculeReactionCounter<"
0397            << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">!\n"
0398            << "Because the objects aren't of the same type!" << G4endl;
0399     G4Exception(G4String("G4VUserMoleculeReactionCounter<"
0400                          + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AbsorbCounter"),
0401                 "BAD_REFERENCE", FatalException, errMsg);
0402   }
0403
0404   if (pCounter->GetType() != GetType()) {
0405     G4ExceptionDescription errMsg;
0406     errMsg << "You are trying to absorb a counter with different type!" << G4endl;
0407     G4Exception(G4String("G4VUserMoleculeReactionCounter<"
0408                          + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AbsorbCounter"),
0409                 "TYPE_DIFF", JustWarning, errMsg);
0410   }
0411
0412   for (auto const& worker_it : pCounter->GetCounterMap()) {
0413     auto [master_it, indexIsNew] =
0414       fCounterMap.emplace(worker_it.first, InnerCounterMapType{fTimeComparer});
0415
0416     G4int currentNumber = 0, previousNumber = 0;
0417     for (auto const& [time, number] : worker_it.second) {
0418       currentNumber = number - previousNumber;
0419       previousNumber = number;
0420
0421       if (master_it->second.empty()) {
0422         master_it->second.emplace(time, currentNumber);
0423       }
0424       else {  // at least one element exists, so we can try to find the closest key
0425         auto it_closest = G4::MoleculeCounter::FindClosestEntryForKey(master_it->second, time);
0426         auto [it, _] = master_it->second.emplace(time, it_closest->second);
0427         do {
0428           it->second += currentNumber;
0429         } while (++it != master_it->second.end());
0430       }
0431     }
0432   }
0433 }
0434
0435 //------------------------------------------------------------------------------
0436
0437 #endif