Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAMoleculeEncounterStepper.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4DNAMoleculeEncounterStepper.cc 93616 2015-10-27 08:59:17Z gcosmo $
27 //
28 // Author: Mathieu Karamitros (kara@cenbg.in2p3.fr)
29 //
30 // WARNING : This class is released as a prototype.
31 // It might strongly evolve or even disapear in the next releases.
32 //
33 // History:
34 // -----------
35 // 10 Oct 2011 M.Karamitros created
36 //
37 // -------------------------------------------------------------------
38 
40 #include "G4VDNAReactionModel.hh"
42 #include "G4H2O.hh"
43 #include "G4memory.hh"
44 #include "G4UnitsTable.hh"
45 #include "G4MoleculeFinder.hh"
47 
48 using namespace std;
49 using namespace CLHEP;
50 
51 //#define DEBUG_MEM
52 
53 #ifdef DEBUG_MEM
54 #include "G4MemStat.hh"
55 using namespace G4MemStat;
56 #endif
57 
58 G4DNAMoleculeEncounterStepper::Utils::Utils(const G4Track& tA,
60  trackA(tA), moleculeB(mB)
61 {
62  moleculeA = GetMolecule(tA);
63  DA = moleculeA->GetDiffusionCoefficient();
64  DB = moleculeB->GetDiffusionCoefficient();
65  Constant = 8 * (DA + DB + 2 * sqrt(DA * DB));
66 }
67 
70  fMolecularReactionTable(
72  fReactionModel(0)
73 {
74  fVerbose = 0;
75  fHasAlreadyReachedNullTime = false;
76 }
77 
78 G4DNAMoleculeEncounterStepper& G4DNAMoleculeEncounterStepper::operator=(const G4DNAMoleculeEncounterStepper& rhs)
79 {
80  if (this == &rhs) return *this;
81  fReactionModel = 0;
82  fVerbose = rhs.fVerbose;
83  fMolecularReactionTable = rhs.fMolecularReactionTable;
84  fHasAlreadyReachedNullTime = false;
85  return *this;
86 }
87 
89 {
90 }
91 
93  G4VITTimeStepComputer(right),
94  fMolecularReactionTable(
96 {
97  fVerbose = right.fVerbose;
98  fMolecularReactionTable = right.fMolecularReactionTable;
99  fReactionModel = 0;
100  fHasAlreadyReachedNullTime = false;
101 }
102 
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 }
126 
127 void G4DNAMoleculeEncounterStepper::InitializeForNewTrack()
128 {
129 // if(fReactants) fReactants = 0 ;
130  if (fReactants) fReactants.reset();
132  fHasAlreadyReachedNullTime = false;
133 }
134 
135 template<typename T>
136  inline bool IsInf(T value)
137  {
138  return std::numeric_limits<T>::has_infinity
139  && value == std::numeric_limits<T>::infinity();
140  }
141 
142 G4double
144  const G4double& userMinTimeStep)
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 }
373 
374 void G4DNAMoleculeEncounterStepper::CheckAndRecordResults(const Utils& utils,
375 #ifdef G4VERBOSE
376  const G4double R,
377 #endif
378  G4KDTreeResultHandle& results)
379 {
380  if (results == 0)
381  {
382 #ifdef G4VERBOSE
383  // DEBUG
384  if (fVerbose > 1)
385  {
386  G4cout << "No molecule " << utils.moleculeB->GetName()
387  << " found to react with " << utils.moleculeA->GetName()
388  << G4endl;
389  }
390 #endif
391  return;
392  }
393 
394  for (results->Rewind(); !results->End(); results->Next())
395  {
396 
397  G4IT* reactiveB = results->GetItem<G4IT>();
398 
399  if (reactiveB == 0)
400  {
401  // DEBUG
402  // G4cout<<"Continue 1"<<G4endl;
403  continue;
404  }
405 
406  G4Track *trackB = reactiveB->GetTrack();
407 
408  if (trackB == 0)
409  {
410  G4ExceptionDescription exceptionDescription;
411  exceptionDescription
412  <<"The reactant B found using the MoleculeFinder does not have a valid "
413  "track attached to it. If this is done on purpose, please do "
414  "not record this molecule in the MoleculeFinder."
415  << G4endl;
416  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
417  "MoleculeEncounterStepper001", FatalErrorInArgument,
418  exceptionDescription);
419  continue;
420  }
421 
422  if (trackB->GetTrackStatus() != fAlive)
423  {
424 // G4ExceptionDescription exceptionDescription;
425 // exceptionDescription
426 // << "The track status of one of the nearby reactants is not fAlive"
427 // << G4endl;
428 // exceptionDescription << "The incomming trackID "
429 // << "(trackA entering in G4DNAMoleculeEncounterStepper and "
430 // << "for which you are looking reactant for) is : "
431 // << utils.trackA.GetTrackID() << "("
432 // << GetMolecule(utils.trackA)->GetName() << ")" << G4endl;
433 // exceptionDescription << "And the trackID of the reactant (trackB) is: "
434 // << trackB->GetTrackID() << "(" << GetMolecule(trackB)->GetName()
435 // << ")" << G4endl;
436 // G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
437 // "MoleculeEncounterStepper002", FatalErrorInArgument,
438 // exceptionDescription);
439  continue;
440  }
441 
442  if (trackB == &utils.trackA)
443  {
444  // DEBUG
445  G4ExceptionDescription exceptionDescription;
446  exceptionDescription
447  << "A track is reacting with itself (which is impossible) ie trackA == trackB"
448  << G4endl;
449  exceptionDescription << "Molecule A (and B) is of type : "
450  << utils.moleculeA->GetName() << " with trackID : "
451  << utils.trackA.GetTrackID() << G4endl;
452 
453  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
454  "MoleculeEncounterStepper003", FatalErrorInArgument,
455  exceptionDescription);
456 
457  }
458 
459  if (fabs(trackB->GetGlobalTime() - utils.trackA.GetGlobalTime())
460  > utils.trackA.GetGlobalTime() * (1 - 1 / 100))
461  {
462  // DEBUG
463  G4ExceptionDescription exceptionDescription;
464  exceptionDescription
465  << "The interacting tracks are not synchronized in time" << G4endl;
466  exceptionDescription
467  << "trackB->GetGlobalTime() != trackA.GetGlobalTime()" << G4endl;
468 
469  exceptionDescription << "trackA : trackID : " << utils.trackA.GetTrackID()
470  << "\t Name :" << utils.moleculeA->GetName()
471  << "\t trackA->GetGlobalTime() = "
472  << G4BestUnit(utils.trackA.GetGlobalTime(), "Time") << G4endl;
473 
474  exceptionDescription << "trackB : trackID : " << trackB->GetTrackID()
475  << "\t Name :" << utils.moleculeB->GetName()
476  << "\t trackB->GetGlobalTime() = "
477  << G4BestUnit(trackB->GetGlobalTime(), "Time") << G4endl;
478 
479  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
480  "MoleculeEncounterStepper004", FatalErrorInArgument,
481  exceptionDescription);
482  }
483 
484 #ifdef G4VERBOSE
485  if(fVerbose > 1)
486  {
487  G4double r2 = results->GetDistanceSqr();
488  G4cout << "\t ************************************************** " << G4endl;
489  G4cout <<"\t Reaction between "
490  << utils.moleculeA->GetName() << " (" << utils.trackA.GetTrackID() << ") "
491  << " & " << utils.moleculeB->GetName() << " (" << trackB->GetTrackID() << "), "
492  << "Interaction Range = "
493  << G4BestUnit(R, "Length")<<G4endl;
494  G4cout <<"\t Real distance between reactants = "
495  << G4BestUnit((utils.trackA.GetPosition() - trackB->GetPosition()).mag(), "Length")<<G4endl;
496  G4cout <<"\t Distance between reactants calculated by nearest neighbor algorithm = "
497  << G4BestUnit(sqrt(r2), "Length")<<G4endl;
498  // G4cout << " ***** " << G4endl;
499  }
500 #endif
501 
502  fReactants->push_back(trackB);
503  }
504 }
Definition: G4IT.hh:88
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4TrackVectorHandle fReactants
const G4ThreeVector & GetPosition() const
G4TrackStatus GetTrackStatus() const
ReturnType & reference_cast(OriginalType &source)
static G4ITFinder * Instance()
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
G4MolecularConfiguration * GetMolecularConfiguration() const
Definition: G4Molecule.cc:576
const G4ThreeVector const G4double const
int G4int
Definition: G4Types.hh:78
MemStat MemoryUsage()
Definition: G4MemStat.cc:55
virtual G4double CalculateStep(const G4Track &, const G4double &)
const G4String & GetName() const
Definition: G4Molecule.cc:356
G4GLOB_DLL std::ostream G4cout
static G4ThreadLocal G4double fUserMinTimeStep
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4ParticleDefinition const G4Material *G4double range
virtual G4double GetReactionRadius(G4MolecularConfiguration *, G4MolecularConfiguration *)=0
virtual void UpdatePositionMap()
G4int GetTrackID() const
G4double GetGlobalTime() const
virtual void Initialise(G4MolecularConfiguration *, const G4Track &)
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:76
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4Track * GetTrack()
Definition: G4IT.hh:219
#define G4endl
Definition: G4ios.hh:61
bool IsInf(T value)
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
const std::vector< G4MolecularConfiguration * > * CanReactWith(G4MolecularConfiguration *) const