Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 *)
 
G4VDNAReactionModelGetReactionModel ()
 
void SetVerbose (int)
 
- Public Member Functions inherited from G4VITTimeStepComputer
 G4VITTimeStepComputer ()
 
virtual ~G4VITTimeStepComputer ()
 
 G4VITTimeStepComputer (const G4VITTimeStepComputer &)
 
G4VITTimeStepComputeroperator= (const G4VITTimeStepComputer &other)
 
virtual void Initialize ()
 
G4TrackVectorHandle GetReactants ()
 
virtual void ResetReactants ()
 
G4double GetSampledMinTimeStep ()
 
void SetReactionTable (const G4ITReactionTable *)
 
const G4ITReactionTableGetReactionTable ()
 

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 G4ITReactionTablefpReactionTable
 
- 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.

68  :
70  fMolecularReactionTable(
71  reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)),
72  fReactionModel(0)
73 {
74  fVerbose = 0;
75  fHasAlreadyReachedNullTime = false;
76 }
const G4ITReactionTable * fpReactionTable
G4DNAMoleculeEncounterStepper::~G4DNAMoleculeEncounterStepper ( )
virtual

Definition at line 88 of file G4DNAMoleculeEncounterStepper.cc.

89 {
90 }
G4DNAMoleculeEncounterStepper::G4DNAMoleculeEncounterStepper ( const G4DNAMoleculeEncounterStepper right)

Definition at line 92 of file G4DNAMoleculeEncounterStepper.cc.

92  :
93  G4VITTimeStepComputer(right),
94  fMolecularReactionTable(
95  reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
96 {
97  fVerbose = right.fVerbose;
98  fMolecularReactionTable = right.fMolecularReactionTable;
99  fReactionModel = 0;
100  fHasAlreadyReachedNullTime = false;
101 }
const G4ITReactionTable * fpReactionTable

Member Function Documentation

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

Implements G4VITTimeStepComputer.

Definition at line 143 of file G4DNAMoleculeEncounterStepper.cc.

145 {
146  // DEBUG
147 // G4cout << "G4MoleculeEncounterStepper::CalculateStep, time :"
148 // << G4ITTrackHolder::Instance()->GetGlobalTime() << G4endl;
149 
150  G4Molecule* moleculeA = GetMolecule(trackA);
151  InitializeForNewTrack();
152  fUserMinTimeStep = userMinTimeStep;
153 
154 #ifdef G4VERBOSE
155  if (fVerbose)
156  {
157  G4cout
158  << "_______________________________________________________________________"
159  << G4endl;
160  G4cout << "G4DNAMoleculeEncounterStepper::CalculateStep" << G4endl;
161  G4cout << "Check done for molecule : " << moleculeA->GetName()
162  << " (" << trackA.GetTrackID() << ") "
163  << G4endl;
164  }
165 #endif
166 
167  //__________________________________________________________________
168  // Retrieve general informations for making reactions
169  G4MolecularConfiguration* molConfA = moleculeA->GetMolecularConfiguration();
170 
171  const vector<G4MolecularConfiguration*>* reactivesVector = fMolecularReactionTable
172  ->CanReactWith(molConfA);
173 
174  if (!reactivesVector)
175  {
176 #ifdef G4VERBOSE
177  // DEBUG
178  if (fVerbose > 1)
179  {
180  G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
181  G4cout << "!!! WARNING" << G4endl;
182  G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
183  "for the reaction because the molecule "
184  << moleculeA->GetName()
185  << " does not have any reactants given in the reaction table."
186  << G4endl;
187  G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
188  }
189 #endif
190  return DBL_MAX;
191  }
192 
193  G4int nbReactives = reactivesVector->size();
194 
195  if (nbReactives == 0)
196  {
197 #ifdef G4VERBOSE
198  // DEBUG
199  if (fVerbose)
200  {
201  // TODO replace with the warning mode of G4Exception
202  G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
203  G4cout << "!!! WARNING" << G4endl;
204  G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
205  "for the reaction because the molecule "
206  << moleculeA->GetName()
207  << " does not have any reactants given in the reaction table."
208  << "This message can also result from a wrong implementation of the reaction table."
209  << G4endl;
210  G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
211  }
212 #endif
213  return DBL_MAX;
214  }
215  // DEBUG
216  // else
217  // {
218  // G4cout << "nb reactants : " << nbReactives << " pour mol "<< moleculeA -> GetName () << G4endl;
219  // for(int k=0 ; k < nbReactives ; k++)
220  // {
221  // G4cout << (*reactivesVector)[k]->GetName() << G4endl;
222  // }
223  // }
224 
225 // fReactants = new vector<G4Track*>();
226  fReactants.reset(new vector<G4Track*>());
227  fReactionModel->Initialise(molConfA, trackA);
228 
229  //__________________________________________________________________
230  // Start looping on possible reactants
231  for (G4int i = 0; i < nbReactives; i++)
232  {
233  G4MolecularConfiguration* moleculeB = (*reactivesVector)[i];
234 
235  //______________________________________________________________
236  // Retrieve reaction range
237  const G4double R = fReactionModel->GetReactionRadius(i);
238 
239  //G4cout << "Reaction range = " << G4BestUnit(R, "Length") << G4endl;
240 
241  //______________________________________________________________
242  // Use KdTree algorithm to find closest reactants
243  G4KDTreeResultHandle resultsNearest(
244  G4MoleculeFinder::Instance()->FindNearest(moleculeA,
245  moleculeB->GetMoleculeID()));
246 
247  if (resultsNearest == 0) continue;
248 
249  G4double r2 = resultsNearest->GetDistanceSqr();
250  Utils utils(trackA, moleculeB);
251 
252  if (r2 <= R * R) // ==> Record in range
253  {
254  // Entering in this condition may due to the fact that molecules are very close
255  // to each other
256  // Therefore, if we only take the nearby reactant into account, it might have already
257  // reacted. Instead, we will take all possible reactants that satisfy the condition r<R
258 
259  if (fHasAlreadyReachedNullTime == false)
260  {
261  fReactants->clear();
262  fHasAlreadyReachedNullTime = true;
263  }
264 
265  fSampledMinTimeStep = 0.;
266  G4KDTreeResultHandle resultsInRange(
267  G4MoleculeFinder::Instance()->FindNearestInRange(moleculeA,
268  moleculeB->GetMoleculeID(),
269  R));
270  CheckAndRecordResults(utils,
271 #ifdef G4VERBOSE
272  R,
273 #endif
274  resultsInRange);
275  }
276  else
277  {
278  G4double r = sqrt(r2);
279  G4double tempMinET = pow(r - R, 2) / utils.Constant;
280  // constant = 16 * (DA + DB + 2*sqrt(DA*DB))
281 
282  //G4cout << tempMinET << G4endl;
283 // G4cout << "fSampledMinTimeStep =" << fSampledMinTimeStep << G4endl;
284 // G4cout << "fUserMinTimeStep =" << fUserMinTimeStep
285 // << " isInf = "<< IsInf(fUserMinTimeStep) << G4endl;
286 
287  if (tempMinET <= fSampledMinTimeStep)
288  {
289  if (fUserMinTimeStep < DBL_MAX/*IsInf(fUserMinTimeStep) == false*/
290  && tempMinET <= fUserMinTimeStep) // ==> Record in range
291  {
293  {
294  fReactants->clear();
295  }
296 
298 
299  G4double range = R + sqrt(fUserMinTimeStep*utils.Constant);
300 
301  G4KDTreeResultHandle resultsInRange (
303  FindNearestInRange(moleculeA,
304  moleculeB->GetMoleculeID(),
305  range));
306 
307  CheckAndRecordResults(utils,
308 #ifdef G4VERBOSE
309  range,
310 #endif
311  resultsInRange);
312  }
313  else // ==> Record nearest
314  {
315  if(tempMinET < fSampledMinTimeStep)
316  // to avoid cases where fSampledMinTimeStep == tempMinET
317  {
318  fSampledMinTimeStep = tempMinET;
319  fReactants->clear();
320  }
321 
322  CheckAndRecordResults(utils,
323 #ifdef G4VERBOSE
324  R,
325 #endif
326  resultsNearest);
327  }
328  }
329  }
330 
331  // DEBUG
332 // if(bool(fReactants))
333 // {
334 // G4cout << "Potential reactions :" << G4endl;
335 // G4cout << GetMolecule(trackA)->GetName()
336 // << " ("<< trackA.GetTrackID()<< ") " << " + ..." << G4endl;
337 // //<< " | " << trackB->GetTrackID() << G4endl;
338 //
339 // for(int j = 0 ; j < (int) fReactants->size() ; j++)
340 // {
341 // G4cout << GetMolecule(fReactants->at(j) )->GetName()
342 // <<" ("<< fReactants->at(j)->GetTrackID() << ")" << G4endl;
343 // }
344 // }
345 
346  }
347 
348 #ifdef G4VERBOSE
349  // DEBUG
350  if (fVerbose)
351  {
352  G4cout << "G4MoleculeEncounterStepper::CalculateStep will finally return :"
353  << G4BestUnit(fSampledMinTimeStep, "Time") << G4endl;
354 
355  if(fVerbose > 1)
356  {
357  G4cout << "Selected reactants for trackA: " << moleculeA->GetName()
358  << " (" << trackA.GetTrackID() << ") are: ";
359 
360  vector<G4Track*>::iterator it;
361  for(it = fReactants->begin(); it != fReactants->end(); it++)
362  {
363  G4Track* trackB = *it;
364  G4cout << GetMolecule(trackB)->GetName() << " ("
365  << trackB->GetTrackID() << ") \t ";
366  }
367  G4cout << G4endl;
368  }
369  }
370 #endif
371  return fSampledMinTimeStep;
372 }
G4TrackVectorHandle fReactants
static G4ITFinder * Instance()
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
G4MolecularConfiguration * GetMolecularConfiguration() const
Definition: G4Molecule.cc:576
int G4int
Definition: G4Types.hh:78
const G4String & GetName() const
Definition: G4Molecule.cc:356
G4GLOB_DLL std::ostream G4cout
static G4ThreadLocal G4double fUserMinTimeStep
const G4ParticleDefinition const G4Material *G4double range
virtual G4double GetReactionRadius(G4MolecularConfiguration *, G4MolecularConfiguration *)=0
G4int GetTrackID() const
virtual void Initialise(G4MolecularConfiguration *, const G4Track &)
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:76
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
const std::vector< G4MolecularConfiguration * > * CanReactWith(G4MolecularConfiguration *) const

Here is the call graph for this function:

G4VDNAReactionModel * G4DNAMoleculeEncounterStepper::GetReactionModel ( )
inline

Definition at line 132 of file G4DNAMoleculeEncounterStepper.hh.

133 {
134  return fReactionModel;
135 }
void G4DNAMoleculeEncounterStepper::Prepare ( )
virtual

Reimplemented from G4VITTimeStepComputer.

Definition at line 103 of file G4DNAMoleculeEncounterStepper.cc.

104 {
105  // DEBUG
106  // G4cout << "G4DNAMoleculeEncounterStepper::PrepareForAllProcessors" << G4endl;
108 
109 #if defined (DEBUG_MEM)
110  MemStat mem_first, mem_second, mem_diff;
111 #endif
112 
113 #if defined (DEBUG_MEM)
114  mem_first = MemoryUsage();
115 #endif
117 
118 #if defined (DEBUG_MEM)
119  mem_second = MemoryUsage();
120  mem_diff = mem_second-mem_first;
121  G4cout << "\t || MEM || G4DNAMoleculeEncounterStepper::Prepare || "
122  "After computing G4ITManager<G4Molecule>::Instance()->"
123  "UpdatePositionMap, diff is : " << mem_diff << G4endl;
124 #endif
125 }
static G4ITFinder * Instance()
MemStat MemoryUsage()
Definition: G4MemStat.cc:55
G4GLOB_DLL std::ostream G4cout
virtual void UpdatePositionMap()
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

void G4DNAMoleculeEncounterStepper::SetReactionModel ( G4VDNAReactionModel reactionModel)
inline

Definition at line 127 of file G4DNAMoleculeEncounterStepper.hh.

128 {
129  fReactionModel = reactionModel;
130 }
void G4DNAMoleculeEncounterStepper::SetVerbose ( int  flag)
inline

Definition at line 137 of file G4DNAMoleculeEncounterStepper.hh.

138 {
139  fVerbose = flag;
140 }

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