G4DNAMoleculeEncounterStepper Class Reference

#include <G4DNAMoleculeEncounterStepper.hh>

Inheritance diagram for G4DNAMoleculeEncounterStepper:

Collaboration diagram for G4DNAMoleculeEncounterStepper:

Classes
class	Utils

Public Member Functions
	G4DNAMoleculeEncounterStepper ()
virtual	~G4DNAMoleculeEncounterStepper ()
	G4DNAMoleculeEncounterStepper (const G4DNAMoleculeEncounterStepper &)
virtual void	Prepare ()
virtual G4double	CalculateStep (const G4Track &, const 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 ()
G4TrackVectorHandle	GetReactants ()
virtual void	ResetReactants ()
G4double	GetSampledMinTimeStep ()
void	SetReactionTable (const G4ITReactionTable *)
const G4ITReactionTable *	GetReactionTable ()

Private Member Functions
void	InitializeForNewTrack ()
void	CheckAndRecordResults (const Utils &, G4KDTreeResultHandle &)
G4DNAMoleculeEncounterStepper &	operator= (const G4DNAMoleculeEncounterStepper &)

Private Attributes
G4bool	fHasAlreadyReachedNullTime
const G4DNAMolecularReactionTable *&	fMolecularReactionTable
G4VDNAReactionModel *	fReactionModel
G4int	fVerbose

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(),
    fMolecularReactionTable(
        reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)),
    fReactionModel(0)
{
  fVerbose = 0;
  fHasAlreadyReachedNullTime = false;
}

~G4DNAMoleculeEncounterStepper()

G4DNAMoleculeEncounterStepper::~G4DNAMoleculeEncounterStepper ( )

virtual

Definition at line 88 of file G4DNAMoleculeEncounterStepper.cc.

{
}

G4DNAMoleculeEncounterStepper() [2/2]

G4DNAMoleculeEncounterStepper::G4DNAMoleculeEncounterStepper

(

const G4DNAMoleculeEncounterStepper &

right

)

Definition at line 92 of file G4DNAMoleculeEncounterStepper.cc.

                                                                                                       :
    G4VITTimeStepComputer(right),
    fMolecularReactionTable(
        reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
{
  fVerbose = right.fVerbose;
  fMolecularReactionTable = right.fMolecularReactionTable;
  fReactionModel = 0;
  fHasAlreadyReachedNullTime = false;
}

Member Function Documentation

CalculateStep()

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

virtual

Implements G4VITTimeStepComputer.

Definition at line 143 of file G4DNAMoleculeEncounterStepper.cc.

{
  // DEBUG
//      G4cout << "G4MoleculeEncounterStepper::CalculateStep, time :"
//  << G4ITTrackHolder::Instance()->GetGlobalTime()  << G4endl;

  G4Molecule* moleculeA = GetMolecule(trackA);
  InitializeForNewTrack();
  fUserMinTimeStep = userMinTimeStep;

#ifdef G4VERBOSE
  if (fVerbose)
  {
    G4cout
    << "_______________________________________________________________________"
    << G4endl;
    G4cout << "G4DNAMoleculeEncounterStepper::CalculateStep" << G4endl;
    G4cout << "Check done for molecule : " << moleculeA->GetName()
           << " (" << trackA.GetTrackID() << ") "
           << G4endl;
  }
#endif

  //__________________________________________________________________
  // Retrieve general informations for making reactions
  G4MolecularConfiguration* molConfA = moleculeA->GetMolecularConfiguration();

  const vector<G4MolecularConfiguration*>* reactivesVector = fMolecularReactionTable
      ->CanReactWith(molConfA);

  if (!reactivesVector)
  {
#ifdef G4VERBOSE
    //    DEBUG
    if (fVerbose > 1)
    {
      G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
      G4cout << "!!! WARNING" << G4endl;
      G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
          "for the reaction because the molecule "
          << moleculeA->GetName()
          << " does not have any reactants given in the reaction table."
          << G4endl;
      G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
    }
#endif
    return DBL_MAX;
  }

  G4int nbReactives = reactivesVector->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 "
      << moleculeA->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;
  }
  // DEBUG
  //    else
  //    {
  //        G4cout << "nb reactants : " << nbReactives << " pour mol "<< moleculeA -> GetName () << G4endl;
  //        for(int k=0 ; k < nbReactives ; k++)
  //        {
  //            G4cout << (*reactivesVector)[k]->GetName() << G4endl;
  //        }
  //    }

// fReactants = new vector<G4Track*>();
  fReactants.reset(new vector<G4Track*>());
  fReactionModel->Initialise(molConfA, trackA);

  //__________________________________________________________________
  // Start looping on possible reactants
  for (G4int i = 0; i < nbReactives; i++)
  {
    G4MolecularConfiguration* moleculeB = (*reactivesVector)[i];

    //______________________________________________________________
    // Retrieve reaction range
    const G4double R = fReactionModel->GetReactionRadius(i);

    //G4cout << "Reaction range = " << G4BestUnit(R, "Length") << G4endl;

    //______________________________________________________________
    // Use KdTree algorithm to find closest reactants
    G4KDTreeResultHandle resultsNearest(
        G4MoleculeFinder::Instance()->FindNearest(moleculeA,
                                                  moleculeB->GetMoleculeID()));

    if (resultsNearest == 0) continue;

    G4double r2 = resultsNearest->GetDistanceSqr();
    Utils utils(trackA, moleculeB);

    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(moleculeA,
                                                           moleculeB->GetMoleculeID(),
                                                           R));
      CheckAndRecordResults(utils,
#ifdef G4VERBOSE
                            R,
#endif
                            resultsInRange);
    }
    else
    {
      G4double r = sqrt(r2);
      G4double tempMinET = pow(r - R, 2) / utils.Constant;
      // constant = 16 * (DA + DB + 2*sqrt(DA*DB))

      //G4cout << tempMinET << G4endl;
//   G4cout << "fSampledMinTimeStep =" << fSampledMinTimeStep << G4endl;
//   G4cout << "fUserMinTimeStep =" << fUserMinTimeStep
//      << " isInf = "<< IsInf(fUserMinTimeStep) << G4endl;

      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.Constant);

          G4KDTreeResultHandle resultsInRange (
                G4MoleculeFinder::Instance()->
                  FindNearestInRange(moleculeA,
                                     moleculeB->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);
        }
      }
    }

    // DEBUG
//  if(bool(fReactants))
//  {
//     G4cout << "Potential reactions :" << G4endl;
//        G4cout << GetMolecule(trackA)->GetName()
//    << " ("<< trackA.GetTrackID()<< ") " << " + ..." << G4endl;
//    //<< " | " << trackB->GetTrackID() << G4endl;
//
//        for(int j = 0 ; j < (int) fReactants->size() ; j++)
//        {
//         G4cout << GetMolecule(fReactants->at(j) )->GetName()
//    <<" ("<< fReactants->at(j)->GetTrackID() << ")" << G4endl;
//        }
//  }

 }

#ifdef G4VERBOSE
  //    DEBUG
  if (fVerbose)
  {
    G4cout << "G4MoleculeEncounterStepper::CalculateStep will finally return :"
           << G4BestUnit(fSampledMinTimeStep, "Time") << G4endl;

    if(fVerbose > 1)
    {
      G4cout << "Selected reactants for trackA: " << moleculeA->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;
}

Here is the call graph for this function:

CheckAndRecordResults()

void G4DNAMoleculeEncounterStepper::CheckAndRecordResults ( const Utils &  utils, G4KDTreeResultHandle &  results  )

private

Definition at line 374 of file G4DNAMoleculeEncounterStepper.cc.

{
  if (results == 0)
  {
#ifdef G4VERBOSE
    // DEBUG
    if (fVerbose > 1)
    {
      G4cout << "No molecule " << utils.moleculeB->GetName()
             << " found to react with " << utils.moleculeA->GetName()
             << G4endl;
    }
#endif
    return;
  }

  for (results->Rewind(); !results->End(); results->Next())
  {

    G4IT* reactiveB = results->GetItem<G4IT>();

    if (reactiveB == 0)
    {
     //  DEBUG
     //  G4cout<<"Continue 1"<<G4endl;
      continue;
    }

    G4Track *trackB = reactiveB->GetTrack();

    if (trackB == 0)
    {
      G4ExceptionDescription exceptionDescription;
      exceptionDescription
      <<"The reactant B found using the MoleculeFinder does not have a valid "
      "track attached to it. If this is done on purpose, please do "
      "not record this molecule in the MoleculeFinder."
      << G4endl;
      G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
          "MoleculeEncounterStepper001", FatalErrorInArgument,
          exceptionDescription);
      continue;
    }

    if (trackB->GetTrackStatus() != fAlive)
    {
//      G4ExceptionDescription exceptionDescription;
//      exceptionDescription
//          << "The track status of one of the nearby reactants is not fAlive"
//          << G4endl;
//      exceptionDescription << "The incomming trackID "
//          << "(trackA entering in G4DNAMoleculeEncounterStepper and "
//          << "for which you are looking reactant for) is : "
//          << utils.trackA.GetTrackID() << "("
//          << GetMolecule(utils.trackA)->GetName() << ")" << G4endl;
//      exceptionDescription << "And the trackID of the reactant (trackB) is: "
//          << trackB->GetTrackID() << "(" << GetMolecule(trackB)->GetName()
//          << ")" << G4endl;
//      G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
//          "MoleculeEncounterStepper002", FatalErrorInArgument,
//          exceptionDescription);
      continue;
    }

    if (trackB == &utils.trackA)
    {
      // DEBUG
      G4ExceptionDescription exceptionDescription;
      exceptionDescription
          << "A track is reacting with itself (which is impossible) ie trackA == trackB"
          << G4endl;
      exceptionDescription << "Molecule A (and B) is of type : "
          << utils.moleculeA->GetName() << " with trackID : "
          << utils.trackA.GetTrackID() << G4endl;

      G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
          "MoleculeEncounterStepper003", FatalErrorInArgument,
          exceptionDescription);

    }

    if (fabs(trackB->GetGlobalTime() - utils.trackA.GetGlobalTime())
        > utils.trackA.GetGlobalTime() * (1 - 1 / 100))
    {
      // DEBUG
      G4ExceptionDescription exceptionDescription;
      exceptionDescription
          << "The interacting tracks are not synchronized in time" << G4endl;
      exceptionDescription
          << "trackB->GetGlobalTime() != trackA.GetGlobalTime()" << G4endl;

      exceptionDescription << "trackA : trackID : " << utils.trackA.GetTrackID()
          << "\t Name :" << utils.moleculeA->GetName()
          << "\t trackA->GetGlobalTime() = "
          << G4BestUnit(utils.trackA.GetGlobalTime(), "Time") << G4endl;

      exceptionDescription << "trackB : trackID : " << trackB->GetTrackID()
          << "\t Name :" << utils.moleculeB->GetName()
          << "\t trackB->GetGlobalTime() = "
          << G4BestUnit(trackB->GetGlobalTime(), "Time") << G4endl;

      G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
          "MoleculeEncounterStepper004", FatalErrorInArgument,
          exceptionDescription);
    }

#ifdef G4VERBOSE
    if(fVerbose > 1)
    {
      G4double r2 = results->GetDistanceSqr();
      G4cout << "\t ************************************************** " << G4endl;
      G4cout <<"\t Reaction between "
      << utils.moleculeA->GetName() << " (" << utils.trackA.GetTrackID() << ") "
      << " & " << utils.moleculeB->GetName() << " (" << trackB->GetTrackID() << "), "
      << "Interaction Range = "
      << G4BestUnit(R, "Length")<<G4endl;
      G4cout <<"\t Real distance between reactants  = "
      << G4BestUnit((utils.trackA.GetPosition() - trackB->GetPosition()).mag(), "Length")<<G4endl;
      G4cout <<"\t Distance between reactants calculated by nearest neighbor algorithm = "
      << G4BestUnit(sqrt(r2), "Length")<<G4endl;
      //            G4cout << " ***** " << G4endl;
    }
#endif

    fReactants->push_back(trackB);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetReactionModel()

G4VDNAReactionModel * G4DNAMoleculeEncounterStepper::GetReactionModel ( )

inline

Definition at line 132 of file G4DNAMoleculeEncounterStepper.hh.

{
  return fReactionModel;
}

InitializeForNewTrack()

void G4DNAMoleculeEncounterStepper::InitializeForNewTrack ( )

private

Definition at line 127 of file G4DNAMoleculeEncounterStepper.cc.

{
// if(fReactants) fReactants = 0 ;
  if (fReactants) fReactants.reset();
  fSampledMinTimeStep = DBL_MAX;
  fHasAlreadyReachedNullTime = false;
}

Here is the caller graph for this function:

operator=()

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

private

Definition at line 78 of file G4DNAMoleculeEncounterStepper.cc.

{
  if (this == &rhs) return *this;
  fReactionModel = 0;
  fVerbose = rhs.fVerbose;
  fMolecularReactionTable = rhs.fMolecularReactionTable;
  fHasAlreadyReachedNullTime = false;
  return *this;
}

Prepare()

void G4DNAMoleculeEncounterStepper::Prepare ( )

virtual

Reimplemented from G4VITTimeStepComputer.

Definition at line 103 of file G4DNAMoleculeEncounterStepper.cc.

{
  // DEBUG
  //    G4cout << "G4DNAMoleculeEncounterStepper::PrepareForAllProcessors" << G4endl;
  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
}

Here is the call graph for this function:

SetReactionModel()

void G4DNAMoleculeEncounterStepper::SetReactionModel ( G4VDNAReactionModel *  reactionModel )

inline

Definition at line 127 of file G4DNAMoleculeEncounterStepper.hh.

{
  fReactionModel = reactionModel;
}

SetVerbose()

void G4DNAMoleculeEncounterStepper::SetVerbose ( int  flag )

inline

Definition at line 137 of file G4DNAMoleculeEncounterStepper.hh.

{
  fVerbose = flag;
}

Member Data Documentation

fHasAlreadyReachedNullTime

G4bool G4DNAMoleculeEncounterStepper::fHasAlreadyReachedNullTime

private

Definition at line 99 of file G4DNAMoleculeEncounterStepper.hh.

fMolecularReactionTable

const G4DNAMolecularReactionTable*& G4DNAMoleculeEncounterStepper::fMolecularReactionTable

private

Definition at line 102 of file G4DNAMoleculeEncounterStepper.hh.

fReactionModel

G4VDNAReactionModel* G4DNAMoleculeEncounterStepper::fReactionModel

private

Definition at line 103 of file G4DNAMoleculeEncounterStepper.hh.

fVerbose

G4int G4DNAMoleculeEncounterStepper::fVerbose

private

Definition at line 104 of file G4DNAMoleculeEncounterStepper.hh.

---

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

• G4DNAMoleculeEncounterStepper.hh
• G4DNAMoleculeEncounterStepper.cc