Geant4_10
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 64374 2012-10-31 16:37:23Z 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 "G4UnitsTable.hh"
44 
45 using namespace std;
46 
49  fMolecularReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)),
50  fReactionModel(0)
51 {
52  fVerbose = 0;
53  fHasAlreadyReachedNullTime = false;
54 }
55 
56 G4DNAMoleculeEncounterStepper& G4DNAMoleculeEncounterStepper::operator=(const G4DNAMoleculeEncounterStepper& rhs)
57 {
58  if(this == &rhs) return *this;
59  fReactionModel = 0;
60  fVerbose = rhs.fVerbose;
61  fMolecularReactionTable = rhs.fMolecularReactionTable;
62  fHasAlreadyReachedNullTime = false;
63  return *this;
64 }
65 
67 {}
68 
70  G4VITTimeStepper(right),
71  fMolecularReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
72 {
73  fVerbose = right.fVerbose ;
74  fMolecularReactionTable = right.fMolecularReactionTable;
75  fReactionModel = 0;
76  fHasAlreadyReachedNullTime = false;
77 }
78 
80 {
81  // DEBUG
82  // G4cout << "G4DNAMoleculeEncounterStepper::PrepareForAllProcessors" << G4endl;
84  G4ITManager<G4Molecule>::Instance()->UpdatePositionMap();
85 }
86 
88 {
89  // DEBUG
90  // G4cout << "G4MoleculeEncounterStepper::CalculateStep, time :" << G4ITTrackHolder::Instance()->GetGlobalTime() << G4endl;
91 
92  G4Molecule* moleculeA = GetMolecule(trackA);
93 
94 
95 #ifdef G4VERBOSE
96  if(fVerbose)
97  {
98  G4cout << "_______________________________________________________________________" << G4endl;
99  G4cout << "G4DNAMoleculeEncounterStepper::CalculateStep" << G4endl;
100  G4cout << "Incident molecule : " << moleculeA->GetName()
101  << " (" << trackA.GetTrackID() << ") "
102  << G4endl;
103  }
104 #endif
105 
106  //__________________________________________________________________
107  // Retrieve general informations for making reactions
108  const vector<const G4Molecule*>* reactivesVector =
109  fMolecularReactionTable -> CanReactWith(moleculeA);
110 
111  if(!reactivesVector)
112  {
113 #ifdef G4VERBOSE
114  // DEBUG
115  if(fVerbose > 1)
116  {
117  G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
118  G4cout << "!!! WARNING" << G4endl;
119  G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity for the reaction because the molecule "
120  << moleculeA->GetName()
121  << " does not have any reactants given in the reaction table."
122  << G4endl;
123  G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
124  }
125 #endif
126  return DBL_MAX;
127  }
128 
129  G4int nbReactives = reactivesVector->size();
130 
131  if(nbReactives == 0)
132  {
133 #ifdef G4VERBOSE
134  // DEBUG
135  if(fVerbose)
136  {
137  G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
138  G4cout << "!!! WARNING" << G4endl;
139  G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity for the reaction because the molecule "
140  << moleculeA->GetName()
141  << " does not have any reactants given in the reaction table."
142  << "This message can also result from a wrong implementation of the reaction table."
143  << G4endl;
144  G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
145  }
146 #endif
147  return DBL_MAX;
148  }
149  // DEBUG
150  // else
151  // {
152  // G4cout << "nb reactants : " << nbReactives << " pour mol "<< moleculeA -> GetName () << G4endl;
153  // for(int k=0 ; k < nbReactives ; k++)
154  // {
155  // G4cout << (*reactivesVector)[k]->GetName() << G4endl;
156  // }
157  // }
158 
159  fUserMinTimeStep = userMinTimeStep ;
160  if(fReactants) fReactants = 0 ;
161  fReactants = new vector<G4Track*>();
162 
164  fHasAlreadyReachedNullTime = false;
165 
166  fReactionModel -> Initialise(moleculeA, trackA) ;
167 
168  //__________________________________________________________________
169  // Start looping on possible reactants
170  for (G4int i=0 ; i<nbReactives ; i++)
171  {
172  const G4Molecule* moleculeB = (*reactivesVector)[i];
173 
174  //______________________________________________________________
175  // Retrieve reaction range
176  G4double R = -1 ; // reaction Range
177  R = fReactionModel -> GetReactionRadius(i);
178 
179  //______________________________________________________________
180  // Use KdTree algorithm to find closest reactants
182  -> FindNearest(moleculeA, moleculeB));
183 
184  RetrieveResults(trackA,moleculeA,moleculeB,R,results, true);
185  }
186 
187 #ifdef G4VERBOSE
188  // DEBUG
189  if(fVerbose)
190  {
191  G4cout << "G4MoleculeEncounterStepper::CalculateStep will finally return :"
192  << G4BestUnit(fSampledMinTimeStep, "Time")<< G4endl;
193 
194  if(fVerbose > 1)
195  {
196  G4cout << "TrackA: " << moleculeA->GetName() << " (" << trackA.GetTrackID() << ") can react with: ";
197 
198  vector<G4Track*>::iterator it;
199  for(it = fReactants->begin() ; it != fReactants->end() ; it++)
200  {
201  G4Track* trackB = *it;
202  G4cout << GetMolecule(trackB)->GetName() << " ("
203  << trackB->GetTrackID() << ") \t ";
204  }
205  G4cout << G4endl;
206  }
207  }
208 #endif
209  return fSampledMinTimeStep ;
210 }
211 
212 
213 void G4DNAMoleculeEncounterStepper::RetrieveResults(const G4Track& trackA, const G4Molecule* moleculeA,
214  const G4Molecule* moleculeB, const G4double R,
215  G4KDTreeResultHandle& results, G4bool iterate)
216 {
217  if(!results)
218  {
219 #ifdef G4VERBOSE
220  // DEBUG
221  if(fVerbose > 1)
222  {
223  G4cout << "No molecule " << moleculeB->GetName()
224  << " found to react with "
225  << moleculeA->GetName()
226  << G4endl;
227  }
228 #endif
229  return ;
230  }
231 
232  G4double DA = moleculeA->GetDiffusionCoefficient() ;
233  G4double DB = moleculeB->GetDiffusionCoefficient() ;
234 
235  for(results->Rewind();
236  !results->End();
237  results->Next())
238  {
239 
240  G4IT* reactiveB = (G4IT*) results->GetItemData() ;
241 
242  if (reactiveB==0)
243  {
244  // DEBUG
245  // G4cout<<"Continue 1"<<G4endl;
246  continue ;
247  }
248 
249  G4Track *trackB = reactiveB->GetTrack();
250 
251  if(trackB == 0)
252  {
253  G4ExceptionDescription exceptionDescription ;
254  exceptionDescription << "The reactant B found using the ITManager does not have a valid track "
255  << " attached to it. If this is done on purpose, please do not record this "
256  << " molecule in the ITManager."
257  << G4endl;
258  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults","MoleculeEncounterStepper001",
259  FatalErrorInArgument,exceptionDescription);
260  continue ;
261  }
262 
263  if (trackB->GetTrackStatus() != fAlive)
264  {
265  G4ExceptionDescription exceptionDescription ;
266  exceptionDescription << "The track status of one of the nearby reactants is not fAlive" << G4endl;
267  exceptionDescription << "The incomming trackID "
268  << "(trackA entering in G4DNAMoleculeEncounterStepper and "
269  << "for which you are looking reactant for) is : "
270  << trackA.GetTrackID() <<"("<< GetMolecule(trackA)->GetName()<<")" << G4endl;
271  exceptionDescription << "And the trackID of the reactant (trackB) is: "
272  << trackB->GetTrackID() <<"("<< GetMolecule(trackB)->GetName()<<")" << G4endl;
273  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults","MoleculeEncounterStepper002",
274  FatalErrorInArgument,exceptionDescription);
275  continue ;
276  }
277 
278  if(trackB == &trackA)
279  {
280  // DEBUG
281  G4ExceptionDescription exceptionDescription ;
282  exceptionDescription << "A track is reacting with itself (which is impossible) ie trackA == trackB"
283  << G4endl ;
284  exceptionDescription << "Molecule A (and B) is of type : "
285  << moleculeA->GetName()
286  << " with trackID : "
287  << trackA.GetTrackID() << G4endl;
288 
289  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults","MoleculeEncounterStepper003",
290  FatalErrorInArgument,exceptionDescription);
291 
292  }
293 
294  if(fabs(trackB->GetGlobalTime() - trackA.GetGlobalTime()) > trackA.GetGlobalTime()*(1-1/100) )
295  {
296  // DEBUG
297  G4ExceptionDescription exceptionDescription;
298  exceptionDescription << "The interacting tracks are not synchronized in time"<< G4endl;
299  exceptionDescription<< "trackB->GetGlobalTime() != trackA.GetGlobalTime()"
300  << G4endl;
301 
302  exceptionDescription << "trackA : trackID : " << trackA.GetTrackID()
303  << "\t Name :" << moleculeA->GetName()
304  <<"\t trackA->GetGlobalTime() = "
305  << G4BestUnit(trackA.GetGlobalTime(), "Time") << G4endl;
306 
307  exceptionDescription << "trackB : trackID : " << trackB->GetTrackID()
308  << "\t Name :" << moleculeB->GetName()
309  << "\t trackB->GetGlobalTime() = "
310  << G4BestUnit(trackB->GetGlobalTime(), "Time")<< G4endl;
311 
312  G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults","MoleculeEncounterStepper004",
313  FatalErrorInArgument,exceptionDescription);
314  }
315 
316  G4double r2 = results->GetDistanceSqr() ;
317 #ifdef G4VERBOSE
318  if(fVerbose > 1)
319  {
320  G4cout << "\t ************************************************** " << G4endl;
321  G4cout <<"\t Reaction between "
322  << moleculeA->GetName() << " (" << trackA.GetTrackID() << ") "
323  << " & " << moleculeB->GetName() << " (" << trackB->GetTrackID() << "), "
324  << "Interaction Range = "
325  << G4BestUnit(R, "Length")<<G4endl;
326  G4cout <<"\t Real distance between reactants = "
327  << G4BestUnit((trackA.GetPosition() - trackB->GetPosition()).mag(), "Length")<<G4endl;
328  G4cout <<"\t Distance between reactants calculated by nearest neighbor algorithm = "
329  << G4BestUnit(sqrt(r2), "Length")<<G4endl;
330 // G4cout << " ***** " << G4endl;
331  }
332 #endif
333 
334  if(r2 <= R*R)
335  {
336  // Entering in this condition may due to the fact that molecules are very close
337  // to each other
338  // Therefore, if we consider only the nearby reactant, this one may have already
339  // react. Instead, we will take all possible reactants that satisfy the condition r<R
340 
341  if(fHasAlreadyReachedNullTime == false)
342  {
343  fReactants->clear();
344  fHasAlreadyReachedNullTime = true;
345  }
346 
347  if(iterate) // First call (call from outside this method)
348  {
349  fSampledMinTimeStep = 0.;
351  -> FindNearestInRange(moleculeA, moleculeB,R));
352  RetrieveResults(trackA, moleculeA, moleculeB, R, results2, false);
353  }
354  else // Other calls (call from this method)
355  {
356  fReactants->push_back(trackB);
357  }
358 
359  continue;
360  }
361  else
362  {
363  G4double r = sqrt(r2);
364  G4double tempMinET = pow(r - R,2)
365  /(16 * (DA + DB + 2*sqrt(DA*DB)));
366 
367  if(tempMinET <= fSampledMinTimeStep)
368  {
369  if(tempMinET <= fUserMinTimeStep)
370  {
372  fReactants->clear();
374  fReactants->push_back(trackB);
375 
376  }
377  else
378  {
379  fSampledMinTimeStep = tempMinET;
380  if(tempMinET < fSampledMinTimeStep)
381  fReactants->clear();
382  fReactants->push_back(trackB);
383  }
384  }
385  }
386  }
387 }
Definition: G4IT.hh:82
static G4ThreadLocal G4double fUserMinTimeStep
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
const G4ThreeVector & GetPosition() const
G4TrackStatus GetTrackStatus() const
ReturnType & reference_cast(OriginalType &source)
G4double GetDiffusionCoefficient() const
Definition: G4Molecule.cc:389
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
const G4ThreeVector const G4double const
int G4int
Definition: G4Types.hh:78
G4TrackVectorHandle fReactants
virtual G4double CalculateStep(const G4Track &, const G4double &)
const G4String & GetName() const
Definition: G4Molecule.cc:243
G4GLOB_DLL std::ostream G4cout
static G4ITManager< T > * Instance()
bool G4bool
Definition: G4Types.hh:79
G4int GetTrackID() const
G4double GetGlobalTime() const
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:67
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
jump r
Definition: plot.C:36
G4Track * GetTrack()
Definition: G4IT.hh:208
#define G4endl
Definition: G4ios.hh:61
G4double fSampledMinTimeStep
double G4double
Definition: G4Types.hh:76
virtual void Prepare()
#define DBL_MAX
Definition: templates.hh:83