109{
110
112 InitializeForNewTrack();
114
115#ifdef G4VERBOSE
116 if (fVerbose != 0)
117 {
119 << "_______________________________________________________________________"
121 G4cout <<
"G4DNAMoleculeEncounterStepper::CalculateStep" <<
G4endl;
122 G4cout <<
"Check done for molecule : " << pMoleculeA->GetName()
125 }
126#endif
127
128
129
130 auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
131
132 const auto pReactantList = fMolecularReactionTable->
CanReactWith(pMolConfA);
133
134 if (pReactantList == nullptr)
135 {
136#ifdef G4VERBOSE
137
138 if (fVerbose > 1)
139 {
142 G4cout <<
"G4MoleculeEncounterStepper::CalculateStep will return infinity "
143 "for the reaction because the molecule "
144 << pMoleculeA->GetName()
145 << " does not have any reactants given in the reaction table."
148 }
149#endif
151 }
152
153 auto nbReactives = (
G4int)pReactantList->size();
154
155 if (nbReactives == 0)
156 {
157#ifdef G4VERBOSE
158
159 if (fVerbose != 0)
160 {
161
164 G4cout <<
"G4MoleculeEncounterStepper::CalculateStep will return infinity "
165 "for the reaction because the molecule "
166 << pMoleculeA->GetName()
167 << " does not have any reactants given in the reaction table."
168 << "This message can also result from a wrong implementation of the reaction table."
171 }
172#endif
174 }
175
176 fReactants = std::make_shared<vector<G4Track*>>();
177 fReactionModel->
Initialise(pMolConfA, trackA);
178
179
180
181 for (
G4int i = 0; i < nbReactives; ++i)
182 {
183 auto pMoleculeB = (*pReactantList)[i];
184
185
186
188
189
190
193 pMoleculeB->GetMoleculeID()));
194
195 if (static_cast<int>(resultsNearest) == 0) continue;
196
197 G4double r2 = resultsNearest->GetDistanceSqr();
198 Utils utils(trackA, pMoleculeB);
199
200 if (r2 <= R * R)
201 {
202
203
204
205
206
207 if (!fHasAlreadyReachedNullTime)
208 {
210 fHasAlreadyReachedNullTime = true;
211 }
212
216 pMoleculeB->GetMoleculeID(),
217 R));
218 CheckAndRecordResults(utils,
219#ifdef G4VERBOSE
220 R,
221#endif
222 resultsInRange);
223 }
224 else
225 {
227 G4double tempMinET = pow(r - R, 2) / utils.fConstant;
228
229
231 {
234 {
236 {
238 }
239
241
243
246 FindNearestInRange(pMoleculeA,
247 pMoleculeB->GetMoleculeID(),
248 range));
249
250 CheckAndRecordResults(utils,
251#ifdef G4VERBOSE
252 range,
253#endif
254 resultsInRange);
255 }
256 else
257 {
259
260 {
263 }
264
265 CheckAndRecordResults(utils,
266#ifdef G4VERBOSE
267 R,
268#endif
269 resultsNearest);
270 }
271 }
272 }
273 }
274
275#ifdef G4VERBOSE
276 if (fVerbose != 0)
277 {
278 G4cout <<
"G4MoleculeEncounterStepper::CalculateStep will finally return :"
280
281 if (fVerbose > 1)
282 {
283 G4cout <<
"Selected reactants for trackA: " << pMoleculeA->GetName()
285
286 vector<G4Track*>::iterator it;
288 {
292 }
294 }
295 }
296#endif
298}
G4Molecule * GetMolecule(const G4Track &track)
G4GLOB_DLL std::ostream G4cout
const ReactantList * CanReactWith(Reactant *) const
static G4ITFinder * Instance()
const G4String & GetName() const override
virtual void Initialise(const G4MolecularConfiguration *, const G4Track &)
virtual G4double GetReactionRadius(const G4MolecularConfiguration *, const G4MolecularConfiguration *)=0
G4TrackVectorHandle fReactants