G4MoleculeCounter.cc

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.          *
// ********************************************************************
//
//
 
#include "G4MoleculeCounter.hh"
 
#include "G4MolecularConfiguration.hh"
#include "G4MoleculeDefinition.hh"
#include "G4MoleculeTable.hh"
#include "G4Scheduler.hh" // TODO: remove this dependency
#include "G4SystemOfUnits.hh"
#include "G4UIcommand.hh"
#include "G4UnitsTable.hh"
 
#include <iomanip>
#include <memory>
 
using namespace std;
 
namespace G4{
namespace MoleculeCounter {
 
bool TimePrecision::operator()(const double& a, const double& b) const
{
    if (std::fabs(a - b) < fPrecision)
    {
        return false;
    }
    
    return a < b;
}
 
G4ThreadLocal double TimePrecision::fPrecision = 0.5 * picosecond;
}
}
 
//------------------------------------------------------------------------------
G4MoleculeCounter* G4MoleculeCounter::Instance()
{
    if (fpInstance == nullptr)
    {
        fpInstance = new G4MoleculeCounter();
    }
    return dynamic_cast<G4MoleculeCounter*>(fpInstance);
}
 
//------------------------------------------------------------------------------
 
G4MoleculeCounter::G4MoleculeCounter()
{
    fVerbose = 0;
    fCheckTimeIsConsistentWithScheduler = true;
}
 
//------------------------------------------------------------------------------
 
G4MoleculeCounter::~G4MoleculeCounter() = default;
 
//------------------------------------------------------------------------------
 
void G4MoleculeCounter::Initialize()
{
    auto mol_iterator = G4MoleculeTable::Instance()->GetConfigurationIterator();
    while ((mol_iterator)())
    {
        if (!IsRegistered(mol_iterator.value()->GetDefinition()))
        {
            continue;
        }
 
        fCounterMap[mol_iterator.value()]; // initialize the second map
    }
}
 
//------------------------------------------------------------------------------
 
void G4MoleculeCounter::SetTimeSlice(double timeSlice)
{
    G4::MoleculeCounter::TimePrecision::fPrecision = timeSlice;
}
 
//------------------------------------------------------------------------------
 
G4bool G4MoleculeCounter::SearchTimeMap(Reactant* molecule)
{
    if (fpLastSearch == nullptr)
    {
        fpLastSearch = std::make_unique<Search>();
    }
    else
    {
        if (fpLastSearch->fLowerBoundSet &&
            fpLastSearch->fLastMoleculeSearched->first == molecule)
        {
            return true;
        }
    }
 
    auto mol_it = fCounterMap.find(molecule);
    fpLastSearch->fLastMoleculeSearched = mol_it;
 
    if (mol_it != fCounterMap.end())
    {
        fpLastSearch->fLowerBoundTime = fpLastSearch->fLastMoleculeSearched->second
                .end();
        fpLastSearch->fLowerBoundSet = true;
    }
    else
    {
        fpLastSearch->fLowerBoundSet = false;
    }
 
    return false;
}
 
//------------------------------------------------------------------------------
 
int G4MoleculeCounter::SearchUpperBoundTime(double time,
                                            bool sameTypeOfMolecule)
{
    auto mol_it = fpLastSearch->fLastMoleculeSearched;
    if (mol_it == fCounterMap.end())
    {
        return 0;
    }
 
    NbMoleculeAgainstTime& timeMap = mol_it->second;
    if (timeMap.empty())
    {
        return 0;
    }
 
    if (sameTypeOfMolecule)
    {
        if (fpLastSearch->fLowerBoundSet && fpLastSearch->fLowerBoundTime != timeMap.end())
        {
            if (fpLastSearch->fLowerBoundTime->first < time)
            {
                auto upperToLast = fpLastSearch->fLowerBoundTime;
                upperToLast++;
 
                if (upperToLast == timeMap.end())
                {
                    return fpLastSearch->fLowerBoundTime->second;
                }
 
                if (upperToLast->first > time)
                {
                    return fpLastSearch->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--;
 
    fpLastSearch->fLowerBoundTime = up_time_it;
    fpLastSearch->fLowerBoundSet = true;
 
    return fpLastSearch->fLowerBoundTime->second;
}
 
//------------------------------------------------------------------------------
 
int G4MoleculeCounter::GetNMoleculesAtTime(Reactant* molecule,
                                           double time)
{
    G4bool sameTypeOfMolecule = SearchTimeMap(molecule);
    return SearchUpperBoundTime(time, sameTypeOfMolecule);
}
 
//------------------------------------------------------------------------------
 
void G4MoleculeCounter::AddAMoleculeAtTime(Reactant* molecule,
                                           G4double time,
                                           const G4ThreeVector* /*position*/,
                                           int number)
{
    if (fDontRegister[molecule->GetDefinition()])
    {
        return;
    }
 
    if (fVerbose != 0)
    {
        G4cout << "G4MoleculeCounter::AddAMoleculeAtTime : " << molecule->GetName()
               << " at time : " << G4BestUnit(time, "Time") << G4endl;
    }
 
    auto counterMap_i = fCounterMap.find(molecule);
 
    if (counterMap_i == fCounterMap.end())
    {
        fCounterMap[molecule][time] = number;
    }
    else if (counterMap_i->second.empty())
    {
        counterMap_i->second[time] = number;
    }
    else
    {
        auto end = counterMap_i->second.rbegin();
 
        if (end->first <= time ||
            fabs(end->first - time) <= G4::MoleculeCounter::TimePrecision::fPrecision)
            // Case 1 = new time comes after last recorded data
            // Case 2 = new time is about the same as the last recorded one
        {
            double newValue = end->second + number;
            counterMap_i->second[time] = newValue;
        }
        else
        {
            //      if(fabs(time - G4Scheduler::Instance()->GetGlobalTime()) >
            //         G4Scheduler::Instance()->GetTimeTolerance())
            {
                G4ExceptionDescription errMsg;
                errMsg << "Time of species "
                       << molecule->GetName() << " is "
                       << G4BestUnit(time, "Time") << " while "
                       << " global time is "
                       << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
                       << G4endl;
                G4Exception("G4MoleculeCounter::AddAMoleculeAtTime",
                            "TIME_DONT_MATCH",
                            FatalException, errMsg);
            }
        }
    }
}
 
//------------------------------------------------------------------------------
 
void G4MoleculeCounter::RemoveAMoleculeAtTime(const G4MolecularConfiguration* pMolecule,
                                              G4double time,
                                              const G4ThreeVector* /*position*/,
                                              int number)
{
    if (fDontRegister[pMolecule->GetDefinition()])
    {
        return;
    }
 
    if (fVerbose != 0)
    {
        G4cout << "G4MoleculeCounter::RemoveAMoleculeAtTime : "
               << pMolecule->GetName() << " at time : " << G4BestUnit(time, "Time")
               << G4endl;
    }
 
    if (fCheckTimeIsConsistentWithScheduler)
    {
        if (fabs(time - G4Scheduler::Instance()->GetGlobalTime()) >
            G4Scheduler::Instance()->GetTimeTolerance())
        {
            G4ExceptionDescription errMsg;
            errMsg << "Time of species "
                   << pMolecule->GetName() << " is "
                   << G4BestUnit(time, "Time") << " while "
                   << " global time is "
                   << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
                   << G4endl;
            G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
                        "TIME_DONT_MATCH",
                        FatalException, errMsg);
        }
    }
 
    NbMoleculeAgainstTime& nbMolPerTime = fCounterMap[pMolecule];
 
    if (nbMolPerTime.empty())
    {
        pMolecule->PrintState();
        Dump();
        G4String errMsg =
                "You are trying to remove molecule " + pMolecule->GetName() +
                " from the counter while this kind of molecules has not been registered yet";
        G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime", "",
                    FatalErrorInArgument, errMsg);
 
        return;
    }
    
    auto it = nbMolPerTime.rbegin();
 
    if (it == nbMolPerTime.rend())
    {
        it--;
 
        G4String errMsg =
                "There was no " + pMolecule->GetName() + " recorded at the time or even before the time asked";
        G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime", "",
                    FatalErrorInArgument, errMsg);
    }
 
    if (time - it->first < -G4::MoleculeCounter::TimePrecision::fPrecision)
    {
        Dump();
        G4ExceptionDescription errMsg;
        errMsg << "Is time going back?? " << pMolecule->GetName()
               << " is being removed at time " << G4BestUnit(time, "Time")
               << " while last recorded time was "
               << G4BestUnit(it->first, "Time") << ".";
        G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
                    "RETURN_TO_THE_FUTUR",
                    FatalErrorInArgument,
                    errMsg);
    }
 
    double finalN = it->second - number;
 
    if (finalN < 0)
    {
        Dump();
        G4ExceptionDescription errMsg;
        errMsg << "After removal of " << number << " species of "
               << pMolecule->GetName() << " the final number at time "
               << G4BestUnit(time, "Time") << " is less than zero and so not valid."
               << " Global time is "
               << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
               << ". Previous selected time is "
               << G4BestUnit(it->first, "Time")
               << G4endl;
        G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
                    "N_INF_0",
                    FatalException, errMsg);
    }
 
    nbMolPerTime[time] = finalN;
}
 
//------------------------------------------------------------------------------
 
G4MoleculeCounter::RecordedMolecules G4MoleculeCounter::GetRecordedMolecules()
{
    if (fVerbose > 1)
    {
        G4cout << "Entering in G4MoleculeCounter::RecordMolecules" << G4endl;
    }
 
    RecordedMolecules output(new ReactantList());
 
    for (const auto & it : fCounterMap)
    {
        output->push_back(it.first);
    }
    return output;
}
 
//------------------------------------------------------------------------------
 
RecordedTimes G4MoleculeCounter::GetRecordedTimes()
{
    RecordedTimes output(new std::set<G4double>);
 
    for(const auto& it : fCounterMap)
    {
        for(const auto& it2 : it.second)
        {
            //time = it2->first;
            output->insert(it2.first);
        }
    }
 
    return output;
}
 
//------------------------------------------------------------------------------
 
// >>DEV<<
//void G4MoleculeCounter::SignalReceiver(G4SpeciesInCM* /*speciesInCM*/,
//                                       size_t moleculeID,
//                                       int /*number*/,
//                                       G4SpeciesInCM::SpeciesChange speciesChange,
//                                       int diff)
//{
//  switch(speciesChange)
//  {
//    case G4SpeciesInCM::eAdd:
//      AddAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
//                         G4Scheduler::Instance()->GetGlobalTime(),
//                         diff);
//      break;
//    case G4SpeciesInCM::eRemove:
//      RemoveAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
//                         G4Scheduler::Instance()->GetGlobalTime(),
//                         diff);
//      break;
//  }
//}
 
//------------------------------------------------------------------------------
 
void G4MoleculeCounter::Dump()
{
    for (const auto& it : fCounterMap)
    {
        auto pReactant = it.first;
 
        G4cout << " --- > For " << pReactant->GetName() << G4endl;
 
        for (const auto& it2 : it.second)
        {
            G4cout << " " << G4BestUnit(it2.first, "Time")
                   << "    " << it2.second << G4endl;
        }
    }
}
 
//------------------------------------------------------------------------------
 
void G4MoleculeCounter::ResetCounter()
{
    if (fVerbose != 0)
    {
        G4cout << " ---> G4MoleculeCounter::ResetCounter" << G4endl;
    }
    fCounterMap.clear();
    fpLastSearch.reset(nullptr);
}
 
//------------------------------------------------------------------------------
 
const NbMoleculeAgainstTime& G4MoleculeCounter::GetNbMoleculeAgainstTime(Reactant* molecule)
{
    return fCounterMap[molecule];
}
 
//------------------------------------------------------------------------------
 
void G4MoleculeCounter::SetVerbose(G4int level)
{
    fVerbose = level;
}
 
//------------------------------------------------------------------------------
 
G4int G4MoleculeCounter::GetVerbose()
{
    return fVerbose;
}
 
//------------------------------------------------------------------------------
 
void G4MoleculeCounter::DontRegister(const G4MoleculeDefinition* molDef)
{
    fDontRegister[molDef] = true;
}
 
//------------------------------------------------------------------------------
 
bool G4MoleculeCounter::IsRegistered(const G4MoleculeDefinition* molDef)
{
    return fDontRegister.find(molDef) == fDontRegister.end();
}
 
//------------------------------------------------------------------------------
 
void G4MoleculeCounter::RegisterAll()
{
    fDontRegister.clear();
}
 
G4bool G4MoleculeCounter::IsTimeCheckedForConsistency() const
{
    return fCheckTimeIsConsistentWithScheduler;
}
 
void G4MoleculeCounter::CheckTimeForConsistency(G4bool flag)
{
    fCheckTimeIsConsistentWithScheduler = flag;
}