ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAMoleculeEncounterStepper.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAMoleculeEncounterStepper.cc
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 //
27 // Author: Mathieu Karamitros (kara@cenbg.in2p3.fr)
28 //
29 // WARNING : This class is released as a prototype.
30 // It might strongly evolve or even disapear in the next releases.
31 //
32 // History:
33 // -----------
34 // 10 Oct 2011 M.Karamitros created
35 //
36 // -------------------------------------------------------------------
37 
39 #include "G4VDNAReactionModel.hh"
41 #include "G4H2O.hh"
42 #include "G4memory.hh"
43 #include "G4UnitsTable.hh"
44 #include "G4MoleculeFinder.hh"
46 
47 using namespace std;
48 using namespace CLHEP;
49 
50 //#define DEBUG_MEM
51 
52 #ifdef DEBUG_MEM
53 #include "G4MemStat.hh"
54 using namespace G4MemStat;
55 #endif
56 
58  const G4MolecularConfiguration* pMoleculeB)
59  : fpTrackA(tA)
60  , fpMoleculeB(pMoleculeB)
61 {
65  fConstant = 8 * (fDA + fDB + 2 * sqrt(fDA * fDB));
66 }
67 
72  , fReactionModel(nullptr)
73  , fVerbose(0)
74 {
75 }
76 
78 
80 {
82 
83 #if defined (DEBUG_MEM)
84  MemStat mem_first, mem_second, mem_diff;
85 #endif
86 
87 #if defined (DEBUG_MEM)
88  mem_first = MemoryUsage();
89 #endif
91 
92 #if defined (DEBUG_MEM)
93  mem_second = MemoryUsage();
94  mem_diff = mem_second - mem_first;
95  G4cout << "\t || MEM || G4DNAMoleculeEncounterStepper::Prepare || "
96  "After computing G4ITManager<G4Molecule>::Instance()->"
97  "UpdatePositionMap, diff is : " << mem_diff << G4endl;
98 #endif
99 }
100 
102 {
103  if (fReactants)
104  {
105  fReactants.reset();
106  }
109 }
110 
111 template<typename T>
112 inline bool IsInf(T value)
113 {
114  return std::numeric_limits<T>::has_infinity
115  && value == std::numeric_limits<T>::infinity();
116 }
117 
118 G4double
120  const G4double& userMinTimeStep)
121 {
122  auto pMoleculeA = GetMolecule(trackA);
124  fUserMinTimeStep = userMinTimeStep;
125 
126 #ifdef G4VERBOSE
127  if (fVerbose)
128  {
129  G4cout
130  << "_______________________________________________________________________"
131  << G4endl;
132  G4cout << "G4DNAMoleculeEncounterStepper::CalculateStep" << G4endl;
133  G4cout << "Check done for molecule : " << pMoleculeA->GetName()
134  << " (" << trackA.GetTrackID() << ") "
135  << G4endl;
136  }
137 #endif
138 
139  //__________________________________________________________________
140  // Retrieve general informations for making reactions
141  auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
142 
143  const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
144 
145  if (!pReactantList)
146  {
147 #ifdef G4VERBOSE
148  // DEBUG
149  if (fVerbose > 1)
150  {
151  G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
152  G4cout << "!!! WARNING" << G4endl;
153  G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
154  "for the reaction because the molecule "
155  << pMoleculeA->GetName()
156  << " does not have any reactants given in the reaction table."
157  << G4endl;
158  G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
159  }
160 #endif
161  return DBL_MAX;
162  }
163 
164  G4int nbReactives = pReactantList->size();
165 
166  if (nbReactives == 0)
167  {
168 #ifdef G4VERBOSE
169  // DEBUG
170  if (fVerbose)
171  {
172  // TODO replace with the warning mode of G4Exception
173  G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
174  G4cout << "!!! WARNING" << G4endl;
175  G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity "
176  "for the reaction because the molecule "
177  << pMoleculeA->GetName()
178  << " does not have any reactants given in the reaction table."
179  << "This message can also result from a wrong implementation of the reaction table."
180  << G4endl;
181  G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
182  }
183 #endif
184  return DBL_MAX;
185  }
186 
187  fReactants.reset(new vector<G4Track*>());
188  fReactionModel->Initialise(pMolConfA, trackA);
189 
190  //__________________________________________________________________
191  // Start looping on possible reactants
192  for (G4int i = 0; i < nbReactives; i++)
193  {
194  auto pMoleculeB = (*pReactantList)[i];
195 
196  //______________________________________________________________
197  // Retrieve reaction range
199 
200  //______________________________________________________________
201  // Use KdTree algorithm to find closest reactants
202  G4KDTreeResultHandle resultsNearest(
203  G4MoleculeFinder::Instance()->FindNearest(pMoleculeA,
204  pMoleculeB->GetMoleculeID()));
205 
206  if (resultsNearest == 0) continue;
207 
208  G4double r2 = resultsNearest->GetDistanceSqr();
209  Utils utils(trackA, pMoleculeB);
210 
211  if (r2 <= R * R) // ==> Record in range
212  {
213  // Entering in this condition may due to the fact that molecules are very close
214  // to each other
215  // Therefore, if we only take the nearby reactant into account, it might have already
216  // reacted. Instead, we will take all possible reactants that satisfy the condition r<R
217 
218  if (fHasAlreadyReachedNullTime == false)
219  {
220  fReactants->clear();
222  }
223 
224  fSampledMinTimeStep = 0.;
225  G4KDTreeResultHandle resultsInRange(
226  G4MoleculeFinder::Instance()->FindNearestInRange(pMoleculeA,
227  pMoleculeB->GetMoleculeID(),
228  R));
229  CheckAndRecordResults(utils,
230 #ifdef G4VERBOSE
231  R,
232 #endif
233  resultsInRange);
234  }
235  else
236  {
237  G4double r = sqrt(r2);
238  G4double tempMinET = pow(r - R, 2) / utils.fConstant;
239  // constant = 16 * (fDA + fDB + 2*sqrt(fDA*fDB))
240 
241  if (tempMinET <= fSampledMinTimeStep)
242  {
243  if (fUserMinTimeStep < DBL_MAX/*IsInf(fUserMinTimeStep) == false*/
244  && tempMinET <= fUserMinTimeStep) // ==> Record in range
245  {
247  {
248  fReactants->clear();
249  }
250 
252 
253  G4double range = R + sqrt(fUserMinTimeStep*utils.fConstant);
254 
255  G4KDTreeResultHandle resultsInRange(
257  FindNearestInRange(pMoleculeA,
258  pMoleculeB->GetMoleculeID(),
259  range));
260 
261  CheckAndRecordResults(utils,
262 #ifdef G4VERBOSE
263  range,
264 #endif
265  resultsInRange);
266  }
267  else // ==> Record nearest
268  {
269  if (tempMinET < fSampledMinTimeStep)
270  // to avoid cases where fSampledMinTimeStep == tempMinET
271  {
272  fSampledMinTimeStep = tempMinET;
273  fReactants->clear();
274  }
275 
276  CheckAndRecordResults(utils,
277 #ifdef G4VERBOSE
278  R,
279 #endif
280  resultsNearest);
281  }
282  }
283  }
284  }
285 
286 #ifdef G4VERBOSE
287  if (fVerbose)
288  {
289  G4cout << "G4MoleculeEncounterStepper::CalculateStep will finally return :"
290  << G4BestUnit(fSampledMinTimeStep, "Time") << G4endl;
291 
292  if (fVerbose > 1)
293  {
294  G4cout << "Selected reactants for trackA: " << pMoleculeA->GetName()
295  << " (" << trackA.GetTrackID() << ") are: ";
296 
297  vector<G4Track*>::iterator it;
298  for (it = fReactants->begin(); it != fReactants->end(); it++)
299  {
300  G4Track* trackB = *it;
301  G4cout << GetMolecule(trackB)->GetName() << " ("
302  << trackB->GetTrackID() << ") \t ";
303  }
304  G4cout << G4endl;
305  }
306  }
307 #endif
308  return fSampledMinTimeStep;
309 }
310 
312 #ifdef G4VERBOSE
313  const G4double R,
314 #endif
315  G4KDTreeResultHandle& results)
316 {
317  if (results == 0)
318  {
319 #ifdef G4VERBOSE
320  if (fVerbose > 1)
321  {
322  G4cout << "No molecule " << utils.fpMoleculeB->GetName()
323  << " found to react with " << utils.fpMoleculeA->GetName()
324  << G4endl;
325  }
326 #endif
327  return;
328  }
329 
330  for (results->Rewind(); !results->End(); results->Next())
331  {
332  G4IT* reactiveB = results->GetItem<G4IT>();
333 
334  if (reactiveB == 0)
335  {
336  continue;
337  }
338 
339  G4Track *trackB = reactiveB->GetTrack();
340 
341  if (trackB == 0)
342  {
343  G4ExceptionDescription exceptionDescription;
344  exceptionDescription
345  << "The reactant B found using the MoleculeFinder does not have a valid "
346  "track attached to it. If this is done on purpose, please do "
347  "not record this molecule in the MoleculeFinder."
348  << G4endl;
349  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
350  "MoleculeEncounterStepper001", FatalErrorInArgument,
351  exceptionDescription);
352  continue;
353  }
354 
355  if (trackB->GetTrackStatus() != fAlive)
356  {
357  continue;
358  }
359 
360  if (trackB == &utils.fpTrackA)
361  {
362  G4ExceptionDescription exceptionDescription;
363  exceptionDescription
364  << "A track is reacting with itself (which is impossible) ie fpTrackA == trackB"
365  << G4endl;
366  exceptionDescription << "Molecule A (and B) is of type : "
367  << utils.fpMoleculeA->GetName() << " with trackID : "
368  << utils.fpTrackA.GetTrackID() << G4endl;
369 
370  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
371  "MoleculeEncounterStepper003", FatalErrorInArgument,
372  exceptionDescription);
373 
374  }
375 
376  if (fabs(trackB->GetGlobalTime() - utils.fpTrackA.GetGlobalTime())
377  > utils.fpTrackA.GetGlobalTime() * (1 - 1 / 100))
378  {
379  // DEBUG
380  G4ExceptionDescription exceptionDescription;
381  exceptionDescription
382  << "The interacting tracks are not synchronized in time" << G4endl;
383  exceptionDescription
384  << "trackB->GetGlobalTime() != fpTrackA.GetGlobalTime()" << G4endl;
385 
386  exceptionDescription << "fpTrackA : trackID : " << utils.fpTrackA.GetTrackID()
387  << "\t Name :" << utils.fpMoleculeA->GetName()
388  << "\t fpTrackA->GetGlobalTime() = "
389  << G4BestUnit(utils.fpTrackA.GetGlobalTime(), "Time") << G4endl;
390 
391  exceptionDescription << "trackB : trackID : " << trackB->GetTrackID()
392  << "\t Name :" << utils.fpMoleculeB->GetName()
393  << "\t trackB->GetGlobalTime() = "
394  << G4BestUnit(trackB->GetGlobalTime(), "Time") << G4endl;
395 
396  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults",
397  "MoleculeEncounterStepper004", FatalErrorInArgument,
398  exceptionDescription);
399  }
400 
401 #ifdef G4VERBOSE
402  if (fVerbose > 1)
403  {
404  G4double r2 = results->GetDistanceSqr();
405  G4cout << "\t ************************************************** " << G4endl;
406  G4cout << "\t Reaction between "
407  << utils.fpMoleculeA->GetName() << " (" << utils.fpTrackA.GetTrackID() << ") "
408  << " & " << utils.fpMoleculeB->GetName() << " (" << trackB->GetTrackID() << "), "
409  << "Interaction Range = "
410  << G4BestUnit(R, "Length") << G4endl;
411  G4cout << "\t Real distance between reactants = "
412  << G4BestUnit((utils.fpTrackA.GetPosition() - trackB->GetPosition()).mag(), "Length") << G4endl;
413  G4cout << "\t Distance between reactants calculated by nearest neighbor algorithm = "
414  << G4BestUnit(sqrt(r2), "Length") << G4endl;
415  }
416 #endif
417 
418  fReactants->push_back(trackB);
419  }
420 }
421 
423 {
424  fReactionModel = pReactionModel;
425 }
426 
428 {
429  return fReactionModel;
430 }
431 
433 {
434  fVerbose = flag;
435 }