G4DNAIndependentReactionTimeStepper Class Reference

#include <G4DNAIndependentReactionTimeStepper.hh>

Inheritance diagram for G4DNAIndependentReactionTimeStepper:

Public Member Functions
	G4DNAIndependentReactionTimeStepper ()
	~G4DNAIndependentReactionTimeStepper () override=default
	G4DNAIndependentReactionTimeStepper (const G4DNAIndependentReactionTimeStepper &)=delete
G4DNAIndependentReactionTimeStepper &	operator= (const G4DNAIndependentReactionTimeStepper &)=delete
void	Prepare () override
G4double	CalculateStep (const G4Track &, const G4double &) override
G4double	CalculateMinTimeStep (G4double, G4double) override
void	SetReactionModel (G4VDNAReactionModel *)
G4VDNAReactionModel *	GetReactionModel ()
std::unique_ptr< G4ITReactionChange >	FindReaction (G4ITReactionSet *pReactionSet, const G4double &currentStepTime=0, const G4double &previousStepTime=0, const G4bool &reachedUserStepTimeLimit=false)
void	SetReactionProcess (G4VITReactionProcess *pReactionProcess)
void	SetVerbose (G4int)
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

Definition at line 51 of file G4DNAIndependentReactionTimeStepper.hh.

Constructor & Destructor Documentation

G4DNAIndependentReactionTimeStepper() [1/2]

G4DNAIndependentReactionTimeStepper::G4DNAIndependentReactionTimeStepper

(

)

Definition at line 54 of file G4DNAIndependentReactionTimeStepper.cc.

  : G4VITTimeStepComputer()
{
  fReactionSet->SortByTime();
}

~G4DNAIndependentReactionTimeStepper()

G4DNAIndependentReactionTimeStepper::~G4DNAIndependentReactionTimeStepper ( )

overridedefault

G4DNAIndependentReactionTimeStepper() [2/2]

G4DNAIndependentReactionTimeStepper::G4DNAIndependentReactionTimeStepper ( const G4DNAIndependentReactionTimeStepper &  )

delete

Member Function Documentation

CalculateMinTimeStep()

G4double G4DNAIndependentReactionTimeStepper::CalculateMinTimeStep ( G4double  , G4double  definedMinTimeStep  )

overridevirtual

Implements G4VITTimeStepComputer.

Definition at line 413 of file G4DNAIndependentReactionTimeStepper.cc.

{
  G4double fTSTimeStep = DBL_MAX;
  fCheckedTracks.clear();
 
  for(auto pTrack : *fpTrackContainer->GetMainList())
  {
    if(pTrack == nullptr)
    {
      G4ExceptionDescription exceptionDescription;
      exceptionDescription << "No track found.";
      G4Exception("G4DNAIndependentReactionTimeStepper::CalculateMinTimeStep",
                  "G4DNAIndependentReactionTimeStepper006",
                  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);
 
        fSampledPositions[pTrack->GetTrackID()] = pTrack->GetPosition();
        for(const auto& it : *fReactants)
        {
          auto pTrackB = it;
          // G4cout<<"position : "<<pTrackB->GetTrackID()<<G4endl;
          fSampledPositions[pTrackB->GetTrackID()] = pTrackB->GetPosition();
        }
        ResetReactants();
      }
    }
    else if(fTSTimeStep == sampledMinTimeStep && G4bool(reactants))
    {
      fReactionSet->AddReactions(fTSTimeStep, const_cast<G4Track*>(pTrack),
                                 reactants);
 
      fSampledPositions[pTrack->GetTrackID()] = pTrack->GetPosition();
      for(const auto& it : *fReactants)
      {
        auto pTrackB = it;
        // G4cout<<"position : "<<pTrackB->GetTrackID()<<G4endl;
        fSampledPositions[pTrackB->GetTrackID()] = pTrackB->GetPosition();
      }
      ResetReactants();
    }
    else if(reactants)
    {
      ResetReactants();
    }
  }
 
  return fTSTimeStep;
}

CalculateStep()

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

overridevirtual

Implements G4VITTimeStepComputer.

Definition at line 77 of file G4DNAIndependentReactionTimeStepper.cc.

{
  auto pMoleculeA = GetMolecule(trackA);
  InitializeForNewTrack();
  fUserMinTimeStep = userMinTimeStep;
  fCheckedTracks.insert(trackA.GetTrackID());
 
#ifdef G4VERBOSE
  if(fVerbose != 0)
  {
    G4cout << "________________________________________________________________"
              "_______"
           << G4endl;
    G4cout << "G4DNAIndependentReactionTimeStepper::CalculateStep" << G4endl;
    G4cout << "Check done for molecule : " << pMoleculeA->GetName() << " ("
           << trackA.GetTrackID() << ") " << G4endl;
  }
#endif
 
  auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
 
  const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
 
  if(pReactantList == nullptr)
  {
#ifdef G4VERBOSE
    if(fVerbose > 1)
    {
      G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
      G4cout << "!!! WARNING" << G4endl;
      G4cout << "G4DNAIndependentReactionTimeStepper::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 = (G4int)pReactantList->size();
 
  if(nbReactives == 0)
  {
#ifdef G4VERBOSE
    //    DEBUG
    if(fVerbose != 0)
    {
      G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
      G4cout << "!!! WARNING" << G4endl;
      G4cout << "G4DNAIndependentReactionTimeStepper::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 = std::make_shared<vector<G4Track*>>();
  fReactionModel->Initialise(pMolConfA, trackA);
  for(G4int i = 0; i < nbReactives; ++i)
  {
    auto pMoleculeB = (*pReactantList)[i];
    G4int key       = pMoleculeB->GetMoleculeID();
 
    // fRCutOff = G4IRTUtils::GetRCutOff(1 * ps);
    fRCutOff = G4IRTUtils::GetRCutOff();
    //______________________________________________________________
    // Retrieve reaction range
    const G4double Reff = fReactionModel->GetReactionRadius(i);
    std::vector<std::pair<G4TrackList::iterator, G4double>> resultIndices;
    resultIndices.clear();
    G4ChemicalMoleculeFinder::Instance()->FindNearestInRange(
      trackA, key, fRCutOff, resultIndices);
 
    if(resultIndices.empty())
    {
      continue;
    }
    for(auto& it : resultIndices)
    {
      G4Track* pTrackB = *(std::get<0>(it));
 
      if(pTrackB == &trackA)
      {
        continue;
      }
      if(pTrackB == nullptr)
      {
        G4ExceptionDescription exceptionDescription;
        exceptionDescription << "No trackB no valid";
        G4Exception("G4DNAIndependentReactionTimeStepper"
                    "::CalculateStep()",
                    "G4DNAIndependentReactionTimeStepper007", FatalException,
                    exceptionDescription);
      }else
      {
        if(fCheckedTracks.find(pTrackB->GetTrackID()) != fCheckedTracks.end())
        {
          continue;
        }
 
        Utils utils(trackA, *pTrackB);
 
        auto pMolB        = GetMolecule(pTrackB);
        auto pMolConfB    = pMolB->GetMolecularConfiguration();
        G4double distance = (trackA.GetPosition() - pTrackB->GetPosition()).mag();
        if(distance * distance < Reff * Reff)
        {
          auto reactionData =
            fMolecularReactionTable->GetReactionData(pMolConfA, pMolConfB);
          if(G4Scheduler::Instance()->GetGlobalTime() == G4Scheduler::Instance()->GetStartTime())
          {
            if(reactionData->GetProbability() > G4UniformRand())
            {
              if(!fHasAlreadyReachedNullTime)
              {
                fReactants->clear();
                fHasAlreadyReachedNullTime = true;
              }
              fSampledMinTimeStep = 0.;
              CheckAndRecordResults(utils);
            }
          }
        }
        else
        {
          G4double tempMinET = GetTimeToEncounter(trackA, *pTrackB);
          if(tempMinET < 0 || tempMinET > G4Scheduler::Instance()->GetEndTime())
          {
            continue;
          }
          if(tempMinET >= fSampledMinTimeStep)
          {
            continue;
          }
          fSampledMinTimeStep = tempMinET;
          fReactants->clear();
          CheckAndRecordResults(utils);
        }
      }
    }
  }
 
#ifdef G4VERBOSE
  if(fVerbose != 0)
  {
    G4cout << "G4DNAIndependentReactionTimeStepper::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().

FindReaction()

std::unique_ptr< G4ITReactionChange > G4DNAIndependentReactionTimeStepper::FindReaction	(	G4ITReactionSet *	pReactionSet,
		const G4double &	currentStepTime = 0,
		const G4double &	previousStepTime = 0,
		const G4bool &	reachedUserStepTimeLimit = false
	)

Definition at line 311 of file G4DNAIndependentReactionTimeStepper.cc.

{
  if(pReactionSet == nullptr)
  {
    return nullptr;
  }
 
  G4ITReactionPerTime& reactionPerTime = pReactionSet->GetReactionsPerTime();
  if(reactionPerTime.empty())
  {
    return nullptr;
  }
 
  for(auto reaction_i                                 = reactionPerTime.begin();
      reaction_i != reactionPerTime.end(); reaction_i = reactionPerTime.begin())
  {
    G4Track* pTrackA = (*reaction_i)->GetReactants().first;
    if(pTrackA->GetTrackStatus() == fStopAndKill)
    {
      continue;
    }
    G4Track* pTrackB = (*reaction_i)->GetReactant(pTrackA);
    if(pTrackB->GetTrackStatus() == fStopAndKill)
    {
      continue;
    }
 
    if(pTrackB == pTrackA)
    {
      G4ExceptionDescription exceptionDescription;
      exceptionDescription << "The IT reaction process sent back a reaction "
                              "between trackA and trackB. ";
      exceptionDescription << "The problem is trackA == trackB";
      G4Exception("G4DNAIndependentReactionTimeStepper::FindReaction",
                  "G4DNAIndependentReactionTimeStepper02", FatalErrorInArgument,
                  exceptionDescription);
    }
    pReactionSet->SelectThisReaction(*reaction_i);
    if(fpReactionProcess != nullptr &&
       fpReactionProcess->TestReactibility(*pTrackA, *pTrackB, currentStepTime,
                                           false))
    {
      if((fSampledPositions.find(pTrackA->GetTrackID()) ==
            fSampledPositions.end() &&
          (fSampledPositions.find(pTrackB->GetTrackID()) ==
           fSampledPositions.end())))
      {
        G4ExceptionDescription exceptionDescription;
        exceptionDescription
          << "The positions of trackA and trackB have no counted ";
        G4Exception("G4DNAIndependentReactionTimeStepper::FindReaction",
                    "G4DNAIndependentReactionTimeStepper0001",
                    FatalErrorInArgument, exceptionDescription);
      }
 
      pTrackA->SetPosition(fSampledPositions[pTrackA->GetTrackID()]);
      pTrackB->SetPosition(fSampledPositions[pTrackB->GetTrackID()]);
      auto pReactionChange =
        fpReactionProcess->MakeReaction(*pTrackA, *pTrackB);
      if(pReactionChange == nullptr)
      {
        return nullptr;
      }
      return pReactionChange;
    }
  }
  return nullptr;
}

GetReactionModel()

G4VDNAReactionModel * G4DNAIndependentReactionTimeStepper::GetReactionModel

(

)

Definition at line 389 of file G4DNAIndependentReactionTimeStepper.cc.

{
  return fReactionModel;
}

operator=()

G4DNAIndependentReactionTimeStepper & G4DNAIndependentReactionTimeStepper::operator= ( const G4DNAIndependentReactionTimeStepper &  )

delete

Prepare()

void G4DNAIndependentReactionTimeStepper::Prepare ( )

overridevirtual

Reimplemented from G4VITTimeStepComputer.

Definition at line 60 of file G4DNAIndependentReactionTimeStepper.cc.

{
  G4VITTimeStepComputer::Prepare();
  fSampledPositions.clear();
  BuildChemicalMoleculeFinder()
}

SetReactionModel()

void G4DNAIndependentReactionTimeStepper::SetReactionModel

(

G4VDNAReactionModel *

pReactionModel

)

Definition at line 383 of file G4DNAIndependentReactionTimeStepper.cc.

{
  fReactionModel = pReactionModel;
}

SetReactionProcess()

void G4DNAIndependentReactionTimeStepper::SetReactionProcess

(

G4VITReactionProcess *

pReactionProcess

)

Definition at line 408 of file G4DNAIndependentReactionTimeStepper.cc.

{
  fpReactionProcess = pReactionProcess;
}

SetVerbose()

void G4DNAIndependentReactionTimeStepper::SetVerbose

(

G4int

flag

)

Definition at line 394 of file G4DNAIndependentReactionTimeStepper.cc.

{
  fVerbose = flag;
}

---

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

• G4DNAIndependentReactionTimeStepper.hh
• G4DNAIndependentReactionTimeStepper.cc