G4DNAMoleculeEncounterStepper Class Reference

#include <G4DNAMoleculeEncounterStepper.hh>

Inheritance diagram for G4DNAMoleculeEncounterStepper:

Collaboration diagram for G4DNAMoleculeEncounterStepper:

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 ()

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::G4DNAMoleculeEncounterStepper

(

)

Definition at line 68 of file G4DNAMoleculeEncounterStepper.cc.

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

G4DNAMoleculeEncounterStepper::~G4DNAMoleculeEncounterStepper ( )

virtual

Definition at line 88 of file G4DNAMoleculeEncounterStepper.cc.

{
}

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

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:

G4VDNAReactionModel * G4DNAMoleculeEncounterStepper::GetReactionModel ( )

inline

Definition at line 132 of file G4DNAMoleculeEncounterStepper.hh.

{
  return fReactionModel;
}

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:

void G4DNAMoleculeEncounterStepper::SetReactionModel ( G4VDNAReactionModel *  reactionModel )

inline

Definition at line 127 of file G4DNAMoleculeEncounterStepper.hh.

{
  fReactionModel = reactionModel;
}

void G4DNAMoleculeEncounterStepper::SetVerbose ( int  flag )

inline

Definition at line 137 of file G4DNAMoleculeEncounterStepper.hh.

{
  fVerbose = flag;
}

---

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

• source/geant4.10.03.p03/source/processes/electromagnetic/dna/models/include/G4DNAMoleculeEncounterStepper.hh
• source/geant4.10.03.p03/source/processes/electromagnetic/dna/models/src/G4DNAMoleculeEncounterStepper.cc