Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/Geant4/G4VUserMoleculeCounter.hh was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

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 G4VUSERMOLECULECOUNTER_HH
0029 #define G4VUSERMOLECULECOUNTER_HH 1
0030 
0031 #include "G4MoleculeCounterManager.hh"
0032 #include "G4MolecularConfiguration.hh"
0033 #include "G4MoleculeCounterTemplates.hh"
0034 #include "G4Scheduler.hh"
0035 #include "G4UnitsTable.hh"
0036 #include "G4VMoleculeCounter.hh"
0037 
0038 //------------------------------------------------------------------------------
0039 
0040 template<class TIndex>
0041 class G4VUserMoleculeCounter : public G4VMoleculeCounter
0042 {
0043     static_assert(std::is_base_of<G4VMoleculeCounter::G4VMoleculeCounterIndex, TIndex>::value,
0044                   "TIndex must be derived from G4VMoleculeCounter::G4VMoleculeCounterIndex! "
0045                   "No forward declaration is allowed.");
0046 
0047   protected:
0048     struct Search;
0049 
0050   public:
0051     G4VUserMoleculeCounter();
0052     G4VUserMoleculeCounter(G4String, MoleculeCounterType = MoleculeCounterType::Other);
0053     ~G4VUserMoleculeCounter() 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<G4VMoleculeCounter::G4VMoleculeCounterIndex> BuildIndex(const G4Track*) const override = 0;
0065     std::unique_ptr<G4VMoleculeCounter::G4VMoleculeCounterIndex> BuildIndex(const G4Track*, const G4StepPoint*) const override = 0;
0066     std::unique_ptr<G4VMoleculeCounter::G4VMoleculeCounterIndex> BuildSimpleIndex(const G4MolecularConfiguration*) const override = 0;
0067 
0068     void AddMolecule(std::unique_ptr<G4VMoleculeCounter::G4VMoleculeCounterIndex>, G4double, G4int = 1) override;
0069     void RemoveMolecule(std::unique_ptr<G4VMoleculeCounter::G4VMoleculeCounterIndex>, G4double, G4int = 1) override;
0070 
0071     std::set<const G4MolecularConfiguration*> GetRecordedMolecules() const override;
0072     std::set<G4double> GetRecordedTimes() const override;
0073 
0074     void SchedulerFinalizedTracking() override;
0075 
0076   protected:
0077     std::map<TIndex, InnerCounterMapType> fCounterMap{};
0078     std::map<TIndex, G4int> fShadowCounterMap{};
0079 
0080   public:
0081     const std::map<TIndex, InnerCounterMapType>& GetCounterMap() const { return fCounterMap; }
0082     std::vector<TIndex> GetMapIndices() const;
0083 
0084     virtual G4int GetNbMoleculesAtTime(const TIndex&, G4double) const;
0085     virtual G4int GetNbMoleculesAtTime(Search&, const TIndex&, G4double) const;
0086     virtual std::vector<G4int> GetNbMoleculesAtTimes(const TIndex&, const std::vector<G4double>&) const;
0087 
0088     //-SEARCH-----------------------------------------------------------------------
0089   protected:
0090     struct Search
0091     {
0092         Search() : fLowerBoundSet(false) {}
0093         typename std::map<TIndex, InnerCounterMapType>::const_iterator fLastIndexSearched;
0094         InnerCounterMapType::const_iterator fLowerBoundTime;
0095         G4bool fLowerBoundSet;
0096     };
0097     G4bool SearchIndexUpdated(Search&, const TIndex&) const;
0098     G4int SearchUpperBoundTime(Search&, G4double, G4bool) const;
0099 };
0100 
0101 //------------------------------------------------------------------------------
0102 
0103 // #include "G4VUserMoleculeCounter.icc"
0104 
0105 //------------------------------------------------------------------------------
0106 
0107 template<typename T>
0108 G4VUserMoleculeCounter<T>::G4VUserMoleculeCounter() : G4VMoleculeCounter()
0109 {}
0110 
0111 //------------------------------------------------------------------------------
0112 
0113 template<typename T>
0114 G4VUserMoleculeCounter<T>::G4VUserMoleculeCounter(G4String name, MoleculeCounterType type)
0115   : G4VMoleculeCounter(name, type)
0116 {}
0117 
0118 //------------------------------------------------------------------------------
0119 
0120 template<typename TIndex>
0121 void G4VUserMoleculeCounter<TIndex>::Initialize()
0122 {
0123   InitializeUser();
0124   fIsInitialized = true;
0125 }
0126 
0127 //------------------------------------------------------------------------------
0128 
0129 template<typename TIndex>
0130 G4int G4VUserMoleculeCounter<TIndex>::GetNbMoleculesAtTime(const TIndex& index, G4double time) const
0131 {
0132   Search search = {};
0133   return GetNbMoleculesAtTime(search, index, time);
0134 }
0135 
0136 //------------------------------------------------------------------------------
0137 
0138 template<typename TIndex>
0139 G4int G4VUserMoleculeCounter<TIndex>::GetNbMoleculesAtTime(Search& search, const TIndex& index,
0140                                                            G4double time) const
0141 {
0142   G4bool sameIndex = !SearchIndexUpdated(search, index);
0143   return SearchUpperBoundTime(search, time, sameIndex);
0144 }
0145 
0146 //------------------------------------------------------------------------------
0147 
0148 template<typename TIndex>
0149 std::vector<G4int>
0150 G4VUserMoleculeCounter<TIndex>::GetNbMoleculesAtTimes(const TIndex& index,
0151                                                       const std::vector<G4double>& times) const
0152 {
0153   Search search = {};
0154   std::vector<G4int> counts = {};
0155   for (auto time : times)
0156     counts.push_back(GetNbMoleculesAtTime(search, index, time));
0157   return counts;
0158 }
0159 
0160 //------------------------------------------------------------------------------
0161 
0162 template<typename TIndex>
0163 void G4VUserMoleculeCounter<TIndex>::AddMolecule(
0164   std::unique_ptr<G4VMoleculeCounter::G4VMoleculeCounterIndex> pIndex, G4double time, G4int number)
0165 {
0166   const TIndex* mapIndex = dynamic_cast<TIndex*>(pIndex.get());
0167 
0168   if(mapIndex == nullptr)
0169   {
0170     G4ExceptionDescription errMsg;
0171     errMsg << "mapIndex is not found "<< G4endl;
0172     G4Exception(G4String("G4VUserMoleculeCounter<"
0173                          + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AddMolecule"),
0174                 "mapIndex == nullptr", FatalException, errMsg);
0175   }else{
0176     if (G4::MoleculeCounter::Contains(fIgnoredMolecules, mapIndex->GetMolecule()->GetDefinition())
0177         || G4::MoleculeCounter::Contains(fIgnoredReactants, mapIndex->GetMolecule()))
0178     {
0179       return;
0180     }
0181 
0182     if (fCheckTimeIsConsistentWithScheduler && G4Scheduler::Instance()->IsRunning()
0183         && std::fabs(time - G4Scheduler::Instance()->GetGlobalTime())
0184              > G4Scheduler::Instance()->GetTimeTolerance())
0185     {
0186       G4ExceptionDescription errMsg;
0187       errMsg << "Time of species " << mapIndex->GetMolecule()->GetName() << " is "
0188              << G4BestUnit(time, "Time") << "while the global time is "
0189              << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time") << G4endl;
0190       G4Exception(G4String("G4VUserMoleculeCounter<"
0191                            + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AddMolecule"),
0192                   "TIME_DONT_MATCH", FatalException, errMsg);
0193     }
0194 
0195     if (IsTimeAboveUpperBound(time)) {
0196       if (fVerbose > 3) {
0197         G4cout << "G4VUserMoleculeCounter<" << G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0198                << ">(" << GetName() << ")::AddMolecule : " << mapIndex->GetMolecule()->GetName()
0199                << " at time : " << G4BestUnit(time, "Time") << G4endl;
0200         G4cout << ":: [IsTimeAboveUpperBound] Skipping since IsTimeAboveUpperBound == true"
0201                << G4endl;
0202       }
0203       return;
0204     }
0205     else if (IsTimeBelowLowerBound(time)) {
0206       // put into shadow counter
0207       auto [it, indexIsNew] = fShadowCounterMap.emplace(*mapIndex, number);
0208       if (!indexIsNew) it->second += number;
0209       if (fVerbose > 3) {
0210         G4cout << ":: [IsTimeBelowLowerBound] Adding " << mapIndex->GetInfo()
0211                << " shadow count: " << it->second - number << " + " << number << G4endl;
0212       }
0213       return;
0214     }
0215 
0216     if (fVerbose > 1) {
0217       G4cout << "G4VUserMoleculeCounter<" << G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0218              << ">(" << GetName() << ")::AddMolecule : " << mapIndex->GetMolecule()->GetName()
0219              << " at time : " << G4BestUnit(time, "Time") << G4endl;
0220     }
0221 
0222     // within time bounds && not-ignored molecule
0223     // -> continue
0224 
0225     auto it_shadow = fShadowCounterMap.find(*mapIndex);
0226     auto [it, indexIsNew] = fCounterMap.emplace(*mapIndex, InnerCounterMapType{fTimeComparer});
0227 
0228     if (it_shadow != fShadowCounterMap.end()) {
0229       // entry found in shadow counter
0230       if (indexIsNew) {
0231         // mapIndex is new, initialize with shadow counter
0232         it->second[fActiveLowerBound] = it_shadow->second;
0233       }
0234       else {
0235         // mapIndex existed, we need to add the shadow count to it
0236         // this happens if we are in a subsequent event and have just crossed over the lower
0237         // activity bound the counter has then already entries from the previous event
0238         InnerCounterMapType::iterator it_time;
0239         G4bool timeIsNew;
0240         std::tie(it_time, timeIsNew) = it->second.emplace(fActiveLowerBound, 0);
0241         do {
0242           it_time->second += it_shadow->second;
0243         } while (++it_time != it->second.end());
0244       }
0245       // either way, remove the shadow count
0246       fShadowCounterMap.erase(it_shadow);
0247     }
0248 
0249     // map for index existed (and was not empty)
0250     if (G4MoleculeCounterManager::Instance()->GetResetCountersBeforeEvent())
0251     // can only do consistency check if the counters are cleared for each event (= chem run)
0252     {
0253       auto end = it->second.rbegin();
0254       auto init_n = end == it->second.rend() ? 0 : end->second;
0255 
0256       auto [it_time, timeIsNew] = it->second.emplace(time, init_n);
0257       it_time->second += number;
0258 
0259       if (fCheckRecordedTimesAreConsistent
0260           && !(end->first <= time
0261                || std::fabs(end->first - time) <= fTimeComparer.GetPrecisionAtTime(time)))
0262       // Case 1 = new time comes after last recorded data
0263       // Case 2 = new time is about the same as the last recorded one
0264       {
0265         G4ExceptionDescription errMsg;
0266         errMsg << "Time of species " << mapIndex->GetMolecule()->GetName() << " is "
0267                << G4BestUnit(time, "Time") << "while the global time is "
0268                << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
0269                << "(last counter time: " << G4BestUnit(end->first, "Time") << ")" << G4endl;
0270         G4Exception(G4String("G4VUserMoleculeCounter<"
0271                              + G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0272                              + ">::AddMolecule"),
0273                     "TIME_DONT_MATCH", FatalException, errMsg);
0274       }
0275     }
0276     else  // counters are not (automatically) reset by manager
0277     {
0278       // since counters are not cleared after chemical run (i.e., after event)
0279       // there will already be numbers in the map, so...
0280       // (1) find the closest time
0281       // (2) emplace entry using closest value as init + number
0282       // (3) add number to all "future" entries as well
0283       if (it->second.empty()) {
0284         it->second.emplace(time, number);
0285       }
0286       else {  // at least one element exists, so we can try to find the closest key
0287         auto it_closest = G4::MoleculeCounter::FindClosestEntryForKey(it->second, time);
0288         auto [it_time, _] = it->second.emplace(time, it_closest->second);
0289         do {
0290           it_time->second += number;
0291         } while (++it_time != it->second.end());
0292       }
0293     }
0294   }
0295 }
0296 
0297 //------------------------------------------------------------------------------
0298 
0299 template<typename TIndex>
0300 void G4VUserMoleculeCounter<TIndex>::RemoveMolecule(
0301   std::unique_ptr<G4VMoleculeCounter::G4VMoleculeCounterIndex> pIndex, G4double time, G4int number)
0302 {
0303   const TIndex* mapIndex = dynamic_cast<TIndex*>(pIndex.get());
0304 
0305   if(mapIndex == nullptr)
0306   {
0307     G4ExceptionDescription errMsg;
0308     errMsg << "mapIndex is not found "<< G4endl;
0309     G4Exception(G4String("G4VUserMoleculeCounter<"
0310                          + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AddMolecule"),
0311                 "mapIndex == nullptr", FatalException, errMsg);
0312   }else{
0313     if (G4::MoleculeCounter::Contains(fIgnoredMolecules, mapIndex->GetMolecule()->GetDefinition())
0314         || G4::MoleculeCounter::Contains(fIgnoredReactants, mapIndex->GetMolecule()))
0315     {
0316       return;
0317     }
0318 
0319     if (fCheckTimeIsConsistentWithScheduler && G4Scheduler::Instance()->IsRunning()
0320         && std::fabs(time - G4Scheduler::Instance()->GetGlobalTime())
0321              > G4Scheduler::Instance()->GetTimeTolerance())
0322     {
0323       G4ExceptionDescription errMsg;
0324       errMsg << "Time of species " << mapIndex->GetMolecule()->GetName() << " is "
0325              << G4BestUnit(time, "Time") << "while the global time is "
0326              << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time") << G4endl;
0327       G4Exception(G4String("G4VUserMoleculeCounter<"
0328                            + G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0329                            + ">::RemoveMolecule"),
0330                   "TIME_DONT_MATCH", FatalException, errMsg);
0331     }
0332 
0333     if (IsTimeBelowLowerBound(time)) {
0334       auto it = fShadowCounterMap.find(*mapIndex);
0335       if (it == fShadowCounterMap.end()) {
0336         G4ExceptionDescription errMsg;
0337         errMsg << "There was no " << mapIndex->GetMolecule()->GetName()
0338                << " recorded at the time or even before the time asked" << G4endl;
0339         G4Exception(G4String("G4VUserMoleculeCounter<"
0340                              + G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0341                              + ">::RemoveMolecule"),
0342                     "", FatalErrorInArgument, errMsg);
0343       }
0344       else {
0345         if (fVerbose > 3) {
0346           G4cout << ":: [IsTimeBelowLowerBound] Removing " << mapIndex->GetInfo()
0347                  << " shadow count: " << it->second << " - " << number << G4endl;
0348         }
0349         it->second -= number;
0350         return;
0351       }
0352     }
0353     else if (IsTimeAboveUpperBound(time)) {
0354       // if the "active" counter was not filled, add the shadow counter to it at the lower bound
0355       // only needed for remove since remove will always be called at the end when the molecule is
0356       // destroyed
0357       auto [it, indexIsNew] = fCounterMap.emplace(*mapIndex, InnerCounterMapType{fTimeComparer});
0358       if (indexIsNew || it->second.empty()) {
0359         auto it_shadow = fShadowCounterMap.find(*mapIndex);
0360         if (it_shadow != fShadowCounterMap.end()) {
0361           it->second[fActiveLowerBound] = it_shadow->second;
0362           if (fVerbose > 3) {
0363             G4cout << ":: [IsTimeAboveUpperBound] Set " << mapIndex->GetInfo()
0364                    << " Map[ActiveLowerBound] with shadow count:" << it_shadow->second << G4endl;
0365           }
0366           fShadowCounterMap.erase(it_shadow);
0367         }
0368         else if (fVerbose > 3) {
0369           G4cout << ":: [IsTimeAboveUpperBound] Not updating with shadow count since"
0370                     " no shadow count was found!"
0371                  << G4endl;
0372         }
0373       }
0374       else if (fVerbose > 3) {
0375         G4cout << ":: [IsTimeAboveUpperBound] Not updating with shadow count"
0376                   " since it already exists "
0377                << G4endl;
0378       }
0379       return;
0380     }
0381 
0382     if (fVerbose > 2) {
0383       G4cout << "G4VUserMoleculeCounter<" << G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0384              << ">(" << GetName() << ")::RemoveMolecule : " << mapIndex->GetMolecule()->GetName()
0385              << " at time : " << G4BestUnit(time, "Time") << G4endl;
0386     }
0387 
0388     // within time bounds && not-ignored molecule
0389     // -> continue
0390 
0391     auto it_shadow = fShadowCounterMap.find(*mapIndex);
0392     auto it = fCounterMap.find(*mapIndex);
0393 
0394     if (it_shadow != fShadowCounterMap.end()) {
0395       // entry found in shadow counter
0396       if (it == fCounterMap.end()) {
0397         // no mapIndex found, initialize with shadow counter
0398         G4bool indexIsNew = false;
0399         std::tie(it, indexIsNew) =
0400           fCounterMap.emplace(*mapIndex, InnerCounterMapType{fTimeComparer});
0401         if (!indexIsNew) {
0402           G4ExceptionDescription errMsg;
0403           errMsg << "We tried to emplace the index after it was found to not exist, but now it "
0404                     "says it existed!?"
0405                  << G4endl;
0406           G4Exception(G4String("G4VUserMoleculeCounter<"
0407                                + G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0408                                + ">::RemoveMolecule"),
0409                       "NONSENSICAL", FatalErrorInArgument, errMsg);
0410         }
0411         it->second[fActiveLowerBound] = it_shadow->second;
0412       }
0413       else {  // it != fCounterMap.end()
0414         // mapIndex exists, we need to add the shadow count to it
0415         // this happens if we are in a subsequent event and have just crossed over the lower
0416         // activity bound the counter has then already entries from the previous event
0417         auto [it_time, _] = it->second.emplace(fActiveLowerBound, 0);
0418         do {
0419           it_time->second += it_shadow->second;
0420         } while (++it_time != it->second.end());
0421       }
0422 
0423       // either way, remove the shadow count
0424       fShadowCounterMap.erase(it_shadow);
0425 
0426       //      if (it_shadow != fShadowCounterMap.end()) {
0427       //      G4bool indexIsNew = false;
0428       //      // auto [it_, indexIsNew] = fCounterMap.emplace(*mapIndex,
0429       //      // InnerCounterMapType{fTimeComparer});
0430       //      std::tie(it, indexIsNew) = fCounterMap.emplace(*mapIndex,
0431       //      InnerCounterMapType{fTimeComparer}); if (!indexIsNew) {
0432       //        G4ExceptionDescription errMsg;
0433       //        errMsg << "We tried to emplace the index after it was found to not exist, but now it
0434       //        "
0435       //                  "says it existed!?"
0436       //               << G4endl;
0437       //        G4Exception(G4String("G4VUserMoleculeCounter<"
0438       //                             + G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0439       //                             + ">::RemoveMolecule"),
0440       //                    "NONSENSICAL", FatalErrorInArgument, errMsg);
0441       //      }
0442       //      it->second[fActiveLowerBound] = it_shadow->second;
0443       //      // it = it_;
0444       //      if (fVerbose > 3) {
0445       //        G4cout << ":: Initialize " << mapIndex->GetInfo()
0446       //               << " Map[ActiveLowerBound] with shadow count:" << it_shadow->second <<
0447       //               G4endl;
0448       //      }
0449       //      fShadowCounterMap.erase(it_shadow);
0450       //    }
0451     }
0452 
0453     InnerCounterMapType& nbMolPerTime = it->second;
0454     InnerCounterMapType::iterator it_time;
0455     G4bool isNewTime = false;
0456     G4double oldTime = 0;
0457 
0458     if (G4MoleculeCounterManager::Instance()->GetResetCountersBeforeEvent()) {
0459       auto end = nbMolPerTime.rbegin();  // get last entry
0460       oldTime = end->first;
0461 
0462       // CHECK: no molecules have been recorded for this index
0463       if (end == nbMolPerTime.rend()) {
0464         if (fVerbose > 2) {
0465           mapIndex->GetMolecule()->PrintState();
0466           Dump();
0467         }
0468         G4ExceptionDescription errMsg;
0469         errMsg << "There was no " << mapIndex->GetMolecule()->GetName()
0470                << " recorded at the time or even before the time asked" << G4endl;
0471         G4Exception("G4VUserMoleculeCounter::RemoveMolecule", "", FatalErrorInArgument, errMsg);
0472       }
0473       // CHECK: current time is less (by more than counter precision) than the most recently
0474       // recorded index
0475       if (fCheckRecordedTimesAreConsistent
0476           && time - end->first < -fTimeComparer.GetPrecisionAtTime(time))
0477       {
0478         if (fVerbose > 2) {
0479           mapIndex->GetMolecule()->PrintState();
0480           Dump();
0481         }
0482         G4ExceptionDescription errMsg;
0483         errMsg << "Is time going back?? " << mapIndex->GetMolecule()->GetName()
0484                << " is being removed at time " << G4BestUnit(time, "Time")
0485                << "while last recorded time was " << G4BestUnit(end->first, "Time") << ".";
0486         G4Exception("G4VUserMoleculeCounter::RemoveMolecule", "RETURN_TO_THE_FUTUR",
0487                     FatalErrorInArgument, errMsg);
0488       }
0489       std::tie(it_time, isNewTime) = nbMolPerTime.emplace(time, end->second);
0490       // auto oldNumber = it_time->second;
0491       it_time->second -= number;
0492       // if (time > 0.001)
0493       //   G4cout << "(t=" << time << ")=" << number << " | old_it(t=" << it_time->first
0494       //          << ") = " << oldNumber << " | timeIsNew=" << isNewTime
0495       //          << " | new_it(t=" << it_time->first << ") = " << it_time->second << G4endl;
0496     }
0497     else {
0498       // since counters are not cleared after chemical run (i.e., after event)
0499       // there will already be numbers in the map, so...
0500       // (1) find the closest time
0501       // (2) emplace entry using closest value as init - number
0502       // (3) remove number from all "future" entries as well
0503       auto it_closest = G4::MoleculeCounter::FindClosestEntryForKey(nbMolPerTime, time);
0504       std::tie(it_time, isNewTime) = nbMolPerTime.emplace(time, it_closest->second);
0505       auto _it = it_time;
0506       do {
0507         _it->second -= number;
0508       } while (++_it != it->second.end());
0509     }
0510 
0511     // Check that count at new time is >= 0
0512     // This currently throws tons of errors for non-basic counters.
0513     //    auto it_time = nbMolPerTime.find(time);
0514     if (it_time == nbMolPerTime.end() || it_time->second < 0) {
0515       if (fVerbose > 2) Dump();
0516       G4ExceptionDescription errMsg;
0517       errMsg << "After removal of " << number << " species of "
0518              << mapIndex->GetMolecule()->GetName() << " the final number at time "
0519              << G4BestUnit(time, "Time") << " is less than zero and so not valid."
0520              << "\nIndex was :" << mapIndex->GetInfo() << "\nGlobal time is "
0521              << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
0522              << "\nPrevious selected time is " << G4BestUnit(oldTime, "Time") << G4endl;
0523       if (fNegativeCountsAreFatal) {
0524         G4Exception(G4String("G4VUserMoleculeCounter<"
0525                              + G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0526                              + ">::RemoveMolecule"),
0527                     "N_INF_0", FatalException, errMsg);
0528       }
0529       else if (fVerbose > 0) {
0530         G4Exception(G4String("G4VUserMoleculeCounter<"
0531                              + G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0532                              + ">::RemoveMolecule"),
0533                     "N_INF_0", JustWarning, errMsg);
0534       }
0535     }
0536   }
0537 }
0538 
0539 //------------------------------------------------------------------------------
0540 
0541 template<typename TIndex>
0542 void G4VUserMoleculeCounter<TIndex>::SchedulerFinalizedTracking()
0543 {
0544   // Add record to fCounterMap for each fShadowCounterMap index unless they exist already
0545   for (auto& it_shadow : fShadowCounterMap) {
0546     auto [it, indexIsNew] =
0547       fCounterMap.emplace(it_shadow.first, InnerCounterMapType{fTimeComparer});
0548     if (indexIsNew || it->second.empty()) {
0549       it->second[fActiveLowerBound] = it_shadow.second;
0550       if (fVerbose > 3) {
0551         G4cout << "G4VUserMoleculeCounter<" << G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0552                << ">(" << GetName() << ")::SchedulerEndedTracking : " << "setting map index '"
0553                << it_shadow.first.GetInfo() << "' from shadow counter to n = " << it_shadow.second
0554                << G4endl;
0555       }
0556     }
0557     else if (fVerbose > 2) {
0558       G4cout << "G4VUserMoleculeCounter<" << G4::MoleculeCounter::GetTemplateTypeName<TIndex>()
0559              << ">(" << GetName() << ")::SchedulerEndedTracking : "
0560              << "encountered dangling shadow counter iterator for index '"
0561              << it_shadow.first.GetInfo() << "'" << G4endl;
0562     }
0563   }
0564   fShadowCounterMap.clear();
0565 }
0566 
0567 //------------------------------------------------------------------------------
0568 
0569 template<typename TIndex>
0570 std::vector<TIndex> G4VUserMoleculeCounter<TIndex>::GetMapIndices() const
0571 {
0572   if (fVerbose > 2) {
0573     G4cout << "Entering in G4VUserMoleculeCounter::GetMapIndices" << G4endl;
0574   }
0575   return G4::MoleculeCounter::GetMapIndices(fCounterMap);
0576 }
0577 
0578 //------------------------------------------------------------------------------
0579 
0580 template<typename T>
0581 std::set<const G4MolecularConfiguration*> G4VUserMoleculeCounter<T>::GetRecordedMolecules() const
0582 {
0583   if (fVerbose > 2) {
0584     G4cout << "Entering in G4MoleculeCounter::RecordMolecules" << G4endl;
0585   }
0586   std::set<const G4MolecularConfiguration*> output{};
0587   for (const auto& it : fCounterMap) {
0588     output.insert(it.first.GetMolecule());
0589   }
0590   return output;
0591 }
0592 
0593 //------------------------------------------------------------------------------
0594 
0595 template<typename T>
0596 std::set<G4double> G4VUserMoleculeCounter<T>::GetRecordedTimes() const
0597 {
0598   return G4::MoleculeCounter::GetRecordedTimes<T>(fCounterMap);
0599 }
0600 
0601 //------------------------------------------------------------------------------
0602 
0603 template<typename T>
0604 void G4VUserMoleculeCounter<T>::Dump() const
0605 {
0606   DumpCounterMapIndices();
0607   G4::MoleculeCounter::DumpCounterMapContents<T>(fCounterMap);
0608 }
0609 
0610 template<typename T>
0611 void G4VUserMoleculeCounter<T>::DumpCounterMapIndices() const
0612 {
0613   G4::MoleculeCounter::DumpCounterMapIndices<T>(fCounterMap);
0614 }
0615 
0616 //------------------------------------------------------------------------------
0617 
0618 template<typename T>
0619 void G4VUserMoleculeCounter<T>::ResetCounter()
0620 {
0621   if (fVerbose > 1) {
0622     G4cout << "G4VUserMoleculeCounter<" << G4::MoleculeCounter::GetTemplateTypeName<T>() << ">("
0623            << GetName() << ")::ResetCounter" << G4endl;
0624   }
0625   fCounterMap.clear();
0626 }
0627 
0628 //------------------------------------------------------------------------------
0629 
0630 template<typename TIndex>
0631 G4bool G4VUserMoleculeCounter<TIndex>::SearchIndexUpdated(Search& search, const TIndex& index) const
0632 {
0633   if (search.fLowerBoundSet && !(search.fLastIndexSearched->first < index)
0634       && !(index < search.fLastIndexSearched->first))
0635   {
0636     return true;
0637   }
0638 
0639   auto mol_it = fCounterMap.find(index);
0640   search.fLastIndexSearched = mol_it;
0641 
0642   if (mol_it != fCounterMap.end()) {
0643     search.fLowerBoundTime = search.fLastIndexSearched->second.end();
0644     search.fLowerBoundSet = true;
0645   }
0646   else {
0647     search.fLowerBoundSet = false;
0648   }
0649 
0650   return false;
0651 }
0652 
0653 //------------------------------------------------------------------------------
0654 
0655 template<typename T>
0656 G4int G4VUserMoleculeCounter<T>::SearchUpperBoundTime(Search& search, G4double time,
0657                                                       G4bool sameIndex) const
0658 {
0659   auto mol_it = search.fLastIndexSearched;
0660   if (mol_it == fCounterMap.end()) {
0661     return 0;
0662   }
0663 
0664   InnerCounterMapType const& timeMap = mol_it->second;
0665   if (timeMap.empty()) {
0666     return 0;
0667   }
0668 
0669   if (sameIndex) {
0670     if (search.fLowerBoundSet && search.fLowerBoundTime != timeMap.end()) {
0671       if (search.fLowerBoundTime->first < time) {
0672         auto upperToLast = search.fLowerBoundTime;
0673         upperToLast++;
0674 
0675         if (upperToLast == timeMap.end()) {
0676           return search.fLowerBoundTime->second;
0677         }
0678 
0679         if (upperToLast->first > time) {
0680           return search.fLowerBoundTime->second;
0681         }
0682       }
0683     }
0684   }
0685 
0686   auto up_time_it = timeMap.upper_bound(time);
0687 
0688   if (up_time_it == timeMap.end()) {
0689     auto last_time = timeMap.rbegin();
0690     return last_time->second;
0691   }
0692   if (up_time_it == timeMap.begin()) {
0693     return 0;
0694   }
0695 
0696   up_time_it--;
0697 
0698   search.fLowerBoundTime = up_time_it;
0699   search.fLowerBoundSet = true;
0700 
0701   return search.fLowerBoundTime->second;
0702 }
0703 
0704 //------------------------------------------------------------------------------
0705 
0706 template<typename TIndex>
0707 void G4VUserMoleculeCounter<TIndex>::AbsorbCounter(const G4VMoleculeCounterInternalBase* pCounterBase)
0708 {
0709   if (pCounterBase == nullptr) {
0710     G4ExceptionDescription errMsg;
0711     errMsg << "Could not cast the pointer to type G4VUserMoleculeCounter<"
0712            << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">!\n"
0713            << "Because the pointer is nullptr!" << G4endl;
0714     G4Exception(G4String("G4VUserMoleculeCounter<"
0715                          + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AbsorbCounter"),
0716                 "BAD_REFERENCE", FatalException, errMsg);
0717   }
0718 
0719   auto pCounter = dynamic_cast<G4VUserMoleculeCounter<TIndex> const*>(pCounterBase);
0720 
0721   if (pCounter == nullptr) {
0722     G4ExceptionDescription errMsg;
0723     errMsg << "Could not cast the pointer to type G4VUserMoleculeCounter<"
0724            << G4::MoleculeCounter::GetTemplateTypeName<TIndex>() << ">!\n"
0725            << "Because the objects aren't of the same type!" << G4endl;
0726     G4Exception(G4String("G4VUserMoleculeCounter<"
0727                          + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AbsorbCounter"),
0728                 "BAD_REFERENCE", FatalException, errMsg);
0729   }
0730 
0731   if (pCounter->GetType() != GetType()) {
0732     G4ExceptionDescription errMsg;
0733     errMsg << "You are trying to absorb a counter with different type!" << G4endl;
0734     G4Exception(G4String("G4VUserMoleculeCounter<"
0735                          + G4::MoleculeCounter::GetTemplateTypeName<TIndex>() + ">::AbsorbCounter"),
0736                 "TYPE_DIFF", JustWarning, errMsg);
0737   }
0738 
0739   for (auto const& worker_it : pCounter->GetCounterMap()) {
0740     auto [master_it, indexIsNew] =
0741       fCounterMap.emplace(worker_it.first, InnerCounterMapType{fTimeComparer});
0742 
0743     G4int currentNumber = 0, previousNumber = 0;
0744     for (auto const& [time, number] : worker_it.second) {
0745       currentNumber = number - previousNumber;
0746       previousNumber = number;
0747 
0748       if (master_it->second.empty()) {
0749         master_it->second.emplace(time, currentNumber);
0750       }
0751       else {  // at least one element exists, so we can try to find the closest key
0752         auto it_closest = G4::MoleculeCounter::FindClosestEntryForKey(master_it->second, time);
0753         auto [it, _] = master_it->second.emplace(time, it_closest->second);
0754         do {
0755           it->second += currentNumber;
0756         } while (++it != master_it->second.end());
0757       }
0758     }
0759   }
0760 }
0761 
0762 //------------------------------------------------------------------------------
0763 
0764 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!