G4DNAMoleculeEncounterStepper Class Reference

#include <G4DNAMoleculeEncounterStepper.hh>

Inheritance diagram for G4DNAMoleculeEncounterStepper:

Public Member Functions
	G4DNAMoleculeEncounterStepper ()
virtual	~G4DNAMoleculeEncounterStepper ()
	G4DNAMoleculeEncounterStepper (const G4DNAMoleculeEncounterStepper &)=delete
G4DNAMoleculeEncounterStepper &	operator= (const G4DNAMoleculeEncounterStepper &)=delete
virtual void	Prepare ()
virtual G4double	CalculateStep (const G4Track &, const G4double &)
virtual G4double	CalculateMinTimeStep (G4double, G4double)
void	SetReactionModel (G4VDNAReactionModel *)
G4VDNAReactionModel *	GetReactionModel ()
void	SetVerbose (int)
Public Member Functions inherited from G4VITTimeStepComputer
	G4VITTimeStepComputer ()
virtual	~G4VITTimeStepComputer ()
	G4VITTimeStepComputer (const G4VITTimeStepComputer &)
G4VITTimeStepComputer &	operator= (const G4VITTimeStepComputer &other)
virtual void	Initialize ()
virtual void	Prepare ()
virtual G4double	CalculateStep (const G4Track &, const G4double &)=0
virtual G4double	CalculateMinTimeStep (G4double, G4double)=0
G4TrackVectorHandle	GetReactants ()
virtual void	ResetReactants ()
G4double	GetSampledMinTimeStep ()
void	SetReactionTable (const G4ITReactionTable *)
const G4ITReactionTable *	GetReactionTable ()

Additional Inherited Members
Static Public Member Functions inherited from G4VITTimeStepComputer
static void	SetTimes (const G4double &, const G4double &)
Protected Attributes inherited from G4VITTimeStepComputer
G4double	fSampledMinTimeStep
G4TrackVectorHandle	fReactants
const G4ITReactionTable *	fpReactionTable
Static Protected Attributes inherited from G4VITTimeStepComputer
static G4ThreadLocal G4double	fCurrentGlobalTime = -1
static G4ThreadLocal G4double	fUserMinTimeStep = -1

Detailed Description

Given a molecule G4DNAMoleculeEncounterStepper will calculate for its possible reactants what will be the minimum encounter time and the associated molecules.*

This model includes dynamical time steps as explained in "Computer-Aided Stochastic Modeling of the Radiolysis of Liquid Water", V. Michalik, M. Begusová, E. A. Bigildeev, Radiation Research, Vol. 149, No. 3 (Mar., 1998), pp. 224-236

Definition at line 70 of file G4DNAMoleculeEncounterStepper.hh.

Constructor & Destructor Documentation

G4DNAMoleculeEncounterStepper() [1/2]

G4DNAMoleculeEncounterStepper::G4DNAMoleculeEncounterStepper

(

)

Definition at line 68 of file G4DNAMoleculeEncounterStepper.cc.

    : G4VITTimeStepComputer()
    , fHasAlreadyReachedNullTime(false)
    , fMolecularReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
    , fReactionModel(nullptr)
    , fVerbose(0)
{
    fpTrackContainer = G4ITTrackHolder::Instance();
    fReactionSet = G4ITReactionSet::Instance();
}

~G4DNAMoleculeEncounterStepper()

G4DNAMoleculeEncounterStepper::~G4DNAMoleculeEncounterStepper ( )

virtualdefault

G4DNAMoleculeEncounterStepper() [2/2]

G4DNAMoleculeEncounterStepper::G4DNAMoleculeEncounterStepper ( const G4DNAMoleculeEncounterStepper &  )

delete

Member Function Documentation

CalculateMinTimeStep()

G4double G4DNAMoleculeEncounterStepper::CalculateMinTimeStep ( G4double  , G4double  definedMinTimeStep  )

virtual

Implements G4VITTimeStepComputer.

Definition at line 441 of file G4DNAMoleculeEncounterStepper.cc.

                                                                                                   {
 
    G4double fTSTimeStep = DBL_MAX;
 
    for (auto pTrack : *fpTrackContainer->GetMainList())
    {
        if (pTrack == nullptr)
        {
            G4ExceptionDescription exceptionDescription;
            exceptionDescription << "No track found.";
            G4Exception("G4Scheduler::CalculateMinStep", "ITScheduler006",
                        FatalErrorInArgument, exceptionDescription);
            continue;
        }
 
        G4TrackStatus trackStatus = pTrack->GetTrackStatus();
        if (trackStatus == fStopAndKill || trackStatus == fStopButAlive)
        {
            continue;
        }
 
        G4double sampledMinTimeStep = CalculateStep(*pTrack, definedMinTimeStep);
        G4TrackVectorHandle reactants = GetReactants();
 
        if (sampledMinTimeStep < fTSTimeStep)
        {
            fTSTimeStep = sampledMinTimeStep;
            fReactionSet->CleanAllReaction();
            if (reactants)
            {
                fReactionSet->AddReactions(fTSTimeStep,
                                           const_cast<G4Track*>(pTrack),
                                           reactants);
                ResetReactants();
            }
         }
         else if (fTSTimeStep == sampledMinTimeStep && bool(reactants))
         {
             fReactionSet->AddReactions(fTSTimeStep,
                                        const_cast<G4Track*>(pTrack),
                                        reactants);
             ResetReactants();
         }
         else if (reactants)
         {
             ResetReactants();
         }
    }
 
    return fTSTimeStep;
}

CalculateStep()

G4double G4DNAMoleculeEncounterStepper::CalculateStep ( const G4Track &  trackA, const G4double &  userMinTimeStep  )

virtual

Implements G4VITTimeStepComputer.

Definition at line 121 of file G4DNAMoleculeEncounterStepper.cc.

{
    auto pMoleculeA = GetMolecule(trackA);
    InitializeForNewTrack();
    fUserMinTimeStep = userMinTimeStep;
 
#ifdef G4VERBOSE
    if (fVerbose)
    {
        G4cout
            << "_______________________________________________________________________"
            << G4endl;
        G4cout << "G4DNAMoleculeEncounterStepper::CalculateStep" << G4endl;
        G4cout << "Check done for molecule : " << pMoleculeA->GetName()
            << " (" << trackA.GetTrackID() << ") "
            << G4endl;
    }
#endif
 
    //__________________________________________________________________
    // Retrieve general informations for making reactions
    auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
 
    const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
 
    if (!pReactantList)
    {
#ifdef G4VERBOSE
        //    DEBUG
        if (fVerbose > 1)
        {
            G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
            G4cout << "!!! WARNING" << G4endl;
            G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
                "for the reaction because the molecule "
                << pMoleculeA->GetName()
                << " does not have any reactants given in the reaction table."
                << G4endl;
            G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
        }
#endif
        return DBL_MAX;
    }
 
    G4int nbReactives = pReactantList->size();
 
    if (nbReactives == 0)
    {
#ifdef G4VERBOSE
        //    DEBUG
        if (fVerbose)
        {
            // TODO replace with the warning mode of G4Exception
            G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
            G4cout << "!!! WARNING" << G4endl;
            G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
                "for the reaction because the molecule "
                << pMoleculeA->GetName()
                << " does not have any reactants given in the reaction table."
                << "This message can also result from a wrong implementation of the reaction table."
                << G4endl;
            G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
        }
#endif
        return DBL_MAX;
    }
 
    fReactants.reset(new vector<G4Track*>());
    fReactionModel->Initialise(pMolConfA, trackA);
 
    //__________________________________________________________________
    // Start looping on possible reactants
    for (G4int i = 0; i < nbReactives; i++)
    {
        auto pMoleculeB = (*pReactantList)[i];
 
        //______________________________________________________________
        // Retrieve reaction range
        const G4double R = fReactionModel->GetReactionRadius(i);
 
        //______________________________________________________________
        // Use KdTree algorithm to find closest reactants
        G4KDTreeResultHandle resultsNearest(
            G4MoleculeFinder::Instance()->FindNearest(pMoleculeA,
                                                      pMoleculeB->GetMoleculeID()));
 
        if (resultsNearest == 0) continue;
 
        G4double r2 = resultsNearest->GetDistanceSqr();
        Utils utils(trackA, pMoleculeB);
 
        if (r2 <= R * R) // ==> Record in range
        {
            // Entering in this condition may due to the fact that molecules are very close
            // to each other
            // Therefore, if we only take the nearby reactant into account, it might have already
            // reacted. Instead, we will take all possible reactants that satisfy the condition r<R
 
            if (fHasAlreadyReachedNullTime == false)
            {
                fReactants->clear();
                fHasAlreadyReachedNullTime = true;
            }
 
            fSampledMinTimeStep = 0.;
            G4KDTreeResultHandle resultsInRange(
                G4MoleculeFinder::Instance()->FindNearestInRange(pMoleculeA,
                                                                 pMoleculeB->GetMoleculeID(),
                                                                 R));
            CheckAndRecordResults(utils,
#ifdef G4VERBOSE
                                  R,
#endif
                                  resultsInRange);
        }
        else
        {
            G4double r = sqrt(r2);
            G4double tempMinET = pow(r - R, 2) / utils.fConstant;
            // constant = 16 * (fDA + fDB + 2*sqrt(fDA*fDB))
 
            if (tempMinET <= fSampledMinTimeStep)
            {
                if (fUserMinTimeStep < DBL_MAX/*IsInf(fUserMinTimeStep) == false*/
                    && tempMinET <= fUserMinTimeStep) // ==> Record in range
                {
                    if (fSampledMinTimeStep > fUserMinTimeStep)
                    {
                        fReactants->clear();
                    }
 
                    fSampledMinTimeStep = fUserMinTimeStep;
 
                    G4double range = R + sqrt(fUserMinTimeStep*utils.fConstant);
 
                    G4KDTreeResultHandle resultsInRange(
                        G4MoleculeFinder::Instance()->
                        FindNearestInRange(pMoleculeA,
                                           pMoleculeB->GetMoleculeID(),
                                           range));
 
                    CheckAndRecordResults(utils,
#ifdef G4VERBOSE
                                          range,
#endif
                                          resultsInRange);
                }
                else // ==> Record nearest
                {
                    if (tempMinET < fSampledMinTimeStep)
                        // to avoid cases where fSampledMinTimeStep == tempMinET
                    {
                        fSampledMinTimeStep = tempMinET;
                        fReactants->clear();
                    }
 
                    CheckAndRecordResults(utils,
#ifdef G4VERBOSE
                                          R,
#endif
                                          resultsNearest);
                }
            }
        }
    }
 
#ifdef G4VERBOSE
    if (fVerbose)
    {
        G4cout << "G4MoleculeEncounterStepper::CalculateStep will finally return :"
            << G4BestUnit(fSampledMinTimeStep, "Time") << G4endl;
 
        if (fVerbose > 1)
        {
            G4cout << "Selected reactants for trackA: " << pMoleculeA->GetName()
                << " (" << trackA.GetTrackID() << ") are: ";
 
            vector<G4Track*>::iterator it;
            for (it = fReactants->begin(); it != fReactants->end(); it++)
            {
                G4Track* trackB = *it;
                G4cout << GetMolecule(trackB)->GetName() << " ("
                    << trackB->GetTrackID() << ") \t ";
            }
            G4cout << G4endl;
        }
    }
#endif
    return fSampledMinTimeStep;
}

Referenced by CalculateMinTimeStep().

GetReactionModel()

G4VDNAReactionModel * G4DNAMoleculeEncounterStepper::GetReactionModel

(

)

Definition at line 431 of file G4DNAMoleculeEncounterStepper.cc.

{
    return fReactionModel;
}

operator=()

G4DNAMoleculeEncounterStepper & G4DNAMoleculeEncounterStepper::operator= ( const G4DNAMoleculeEncounterStepper &  )

delete

Prepare()

void G4DNAMoleculeEncounterStepper::Prepare ( )

virtual

Reimplemented from G4VITTimeStepComputer.

Definition at line 81 of file G4DNAMoleculeEncounterStepper.cc.

{
    G4VITTimeStepComputer::Prepare();
 
#if defined (DEBUG_MEM)
    MemStat mem_first, mem_second, mem_diff;
#endif
 
#if defined (DEBUG_MEM)
    mem_first = MemoryUsage();
#endif
    G4MoleculeFinder::Instance()->UpdatePositionMap();
 
#if defined (DEBUG_MEM)
    mem_second = MemoryUsage();
    mem_diff = mem_second - mem_first;
    G4cout << "\t || MEM || G4DNAMoleculeEncounterStepper::Prepare || "
        "After computing G4ITManager<G4Molecule>::Instance()->"
        "UpdatePositionMap, diff is : " << mem_diff << G4endl;
#endif
}

SetReactionModel()

void G4DNAMoleculeEncounterStepper::SetReactionModel

(

G4VDNAReactionModel *

pReactionModel

)

Definition at line 426 of file G4DNAMoleculeEncounterStepper.cc.

{
    fReactionModel = pReactionModel;
}

SetVerbose()

void G4DNAMoleculeEncounterStepper::SetVerbose

(

int

flag

)

Definition at line 436 of file G4DNAMoleculeEncounterStepper.cc.

{
    fVerbose = flag;
}

---

The documentation for this class was generated from the following files:

• G4DNAMoleculeEncounterStepper.hh
• G4DNAMoleculeEncounterStepper.cc