| Geant4 Cross Reference |
1 // ******************************************* 1 // ********************************************************************
2 // * License and Disclaimer 2 // * License and Disclaimer *
3 // * 3 // * *
4 // * The Geant4 software is copyright of th 4 // * The Geant4 software is copyright of the Copyright Holders of *
5 // * the Geant4 Collaboration. It is provided 5 // * the Geant4 Collaboration. It is provided under the terms and *
6 // * conditions of the Geant4 Software License 6 // * conditions of the Geant4 Software License, included in the file *
7 // * LICENSE and available at http://cern.ch/ 7 // * LICENSE and available at http://cern.ch/geant4/license . These *
8 // * include a list of copyright holders. 8 // * include a list of copyright holders. *
9 // * 9 // * *
10 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
11 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
12 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
13 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
14 // * use. Please see the license in the file 14 // * use. Please see the license in the file LICENSE and URL above *
15 // * for the full disclaimer and the limitatio 15 // * for the full disclaimer and the limitation of liability. *
16 // * 16 // * *
17 // * This code implementation is the result 17 // * This code implementation is the result of the scientific and *
18 // * technical work of the GEANT4 collaboratio 18 // * technical work of the GEANT4 collaboration. *
19 // * By using, copying, modifying or distri 19 // * By using, copying, modifying or distributing the software (or *
20 // * any work based on the software) you ag 20 // * any work based on the software) you agree to acknowledge its *
21 // * use in resulting scientific publicati 21 // * use in resulting scientific publications, and indicate your *
22 // * acceptance of all terms of the Geant4 Sof 22 // * acceptance of all terms of the Geant4 Software license. *
23 // ******************************************* 23 // ********************************************************************
24 // 24 //
25 25
>> 26 #include "G4DNAMolecularReactionTable.hh"
26 #include "G4DNAGillespieDirectMethod.hh" 27 #include "G4DNAGillespieDirectMethod.hh"
27 #include "Randomize.hh" 28 #include "Randomize.hh"
28 #include "G4PhysicalConstants.hh" 29 #include "G4PhysicalConstants.hh"
29 #include <memory> 30 #include <memory>
30 #include <tuple> 31 #include <tuple>
31 #include "G4DNAEventSet.hh" 32 #include "G4DNAEventSet.hh"
32 #include "G4UnitsTable.hh" 33 #include "G4UnitsTable.hh"
33 #include "G4DNAScavengerMaterial.hh" 34 #include "G4DNAScavengerMaterial.hh"
34 #include "G4Scheduler.hh" 35 #include "G4Scheduler.hh"
35 #include "G4DNAMolecularReactionTable.hh" << 36 #include <cassert>
36 #include "G4ChemEquilibrium.hh" <<
37 37
38 G4DNAGillespieDirectMethod::G4DNAGillespieDire 38 G4DNAGillespieDirectMethod::G4DNAGillespieDirectMethod()
39 : fMolecularReactions(G4DNAMolecularReaction 39 : fMolecularReactions(G4DNAMolecularReactionTable::Instance())
>> 40 , fpMesh(nullptr)
>> 41 , fTimeStep(0)
>> 42 , fpEventSet(nullptr)
>> 43 , fVerbose(0)
>> 44 , fpScavengerMaterial(nullptr)
40 {} 45 {}
41 46
42 G4DNAGillespieDirectMethod::~G4DNAGillespieDir 47 G4DNAGillespieDirectMethod::~G4DNAGillespieDirectMethod() = default;
43 48
44 void G4DNAGillespieDirectMethod::SetEventSet(G 49 void G4DNAGillespieDirectMethod::SetEventSet(G4DNAEventSet* pEventSet)
45 { 50 {
46 fpEventSet = pEventSet; 51 fpEventSet = pEventSet;
47 } 52 }
48 53
49 //#define DEBUG 1 54 //#define DEBUG 1
50 55
51 G4double G4DNAGillespieDirectMethod::VolumeOfN << 56 G4double G4DNAGillespieDirectMethod::VolumeOfNode(const Index& index)
52 { 57 {
53 auto box = std::get<1>(voxel); << 58 auto LengthY = fpMesh->GetBoundingBox(index).Getyhi() -
54 auto LengthY = box.Getyhi() - box.Getylo(); << 59 fpMesh->GetBoundingBox(index).Getylo();
55 auto LengthX = box.Getxhi() - box.Getxlo(); << 60 auto LengthX = fpMesh->GetBoundingBox(index).Getxhi() -
56 auto LengthZ = box.Getzhi() - box.Getzlo(); << 61 fpMesh->GetBoundingBox(index).Getxlo();
57 G4double V = LengthY * LengthX * LengthZ; << 62 auto LengthZ = fpMesh->GetBoundingBox(index).Getzhi() -
58 if(V <= 0) << 63 fpMesh->GetBoundingBox(index).Getzlo();
59 { << 64 G4double V = LengthY * LengthX * LengthZ;
60 G4ExceptionDescription exceptionDescriptio << 65 assert(V > 0);
61 exceptionDescription << "V > 0 !! "; <<
62 G4Exception("G4DNAGillespieDirectMethod::V <<
63 "G4DNAGillespieDirectMethod03" <<
64 exceptionDescription); <<
65 } <<
66 return V; 66 return V;
67 } 67 }
68 G4double G4DNAGillespieDirectMethod::Propensit << 68 G4double G4DNAGillespieDirectMethod::PropensityFunction(const Index& index,
69 69 MolType moleType)
70 { 70 {
71 if(moleType->GetDiffusionCoefficient() == 0) 71 if(moleType->GetDiffusionCoefficient() == 0)
72 { 72 {
73 return 0.; 73 return 0.;
74 } 74 }
75 const auto& node = std::get<2>(voxel); << 75 const auto& node = fpMesh->GetVoxelMapList(index);
76 const auto& box = std::get<1>(voxel); << 76 G4double alpha = 0;
77 << 77 auto it = node.find(moleType);
78 G4double alpha = 0; <<
79 auto it = node.find(moleType); <<
80 if(it != node.end()) 78 if(it != node.end())
81 { 79 {
82 auto LengthY = box.Getyhi() - box.Getylo() << 80 auto LengthY = fpMesh->GetBoundingBox(index).Getyhi() -
83 G4double d = it->first->GetDiffusionCoef << 81 fpMesh->GetBoundingBox(index).Getylo();
84 alpha = d * it->second; << 82 G4double d = it->first->GetDiffusionCoefficient() / std::pow(LengthY, 2);
>> 83 alpha = d * it->second;
85 84
86 #ifdef DEBUG 85 #ifdef DEBUG
87 G4cout << it->first->GetName() << " " << i 86 G4cout << it->first->GetName() << " " << it->second
88 << " D : " << it->first->GetDiffusi 87 << " D : " << it->first->GetDiffusionCoefficient()
89 << " LengthY : " << LengthY << " Pr 88 << " LengthY : " << LengthY << " PropensityFunction : " << alpha
90 << G4endl; 89 << G4endl;
91 #endif 90 #endif
92 } 91 }
93 return alpha; 92 return alpha;
94 } 93 }
95 94
96 G4double G4DNAGillespieDirectMethod::Propensit << 95 G4double G4DNAGillespieDirectMethod::PropensityFunction(const Index& index,
97 96 ReactionData* data)
98 { 97 {
99 G4double value; 98 G4double value;
100 auto ConfA = data->GetReactant 99 auto ConfA = data->GetReactant1();
101 auto ConfB = data->GetReactant 100 auto ConfB = data->GetReactant2();
102 G4double scavengerNumber = 0; 101 G4double scavengerNumber = 0;
103 auto typeANumber = FindScavenging(vo << 102 auto typeANumber = FindScavenging(index, ConfA, scavengerNumber)
104 ? scavengerNumb << 103 ? scavengerNumber
105 : ComputeNumber << 104 : ComputeNumberInNode(index, ConfA);
106 105
107 auto typeBNumber = FindScavenging(voxel, Con << 106 auto typeBNumber = FindScavenging(index, ConfB, scavengerNumber)
108 ? scavengerNumber 107 ? scavengerNumber
109 : ComputeNumberInNode(v << 108 : ComputeNumberInNode(index, ConfB);
110 109
111 if(typeANumber == 0 || typeBNumber == 0) 110 if(typeANumber == 0 || typeBNumber == 0)
112 { 111 {
113 return 0; 112 return 0;
114 } 113 }
115 114
116 auto k = 115 auto k =
117 data->GetObservedReactionRateConstant() / << 116 data->GetObservedReactionRateConstant() / (Avogadro * VolumeOfNode(index));
118 if(ConfA == ConfB) 117 if(ConfA == ConfB)
119 { 118 {
120 value = typeANumber * (typeBNumber - 1) * 119 value = typeANumber * (typeBNumber - 1) * k;
121 } 120 }
122 else 121 else
123 { 122 {
124 value = typeANumber * typeBNumber * k; 123 value = typeANumber * typeBNumber * k;
125 } 124 }
126 125
127 if(value < 0) 126 if(value < 0)
128 { 127 {
129 G4ExceptionDescription exceptionDescriptio << 128 G4cout << "G4DNAGillespieDirectMethod::PropensityFunction for : "
130 exceptionDescription << 129 << ConfA->GetName() << "(" << typeANumber << ") + "
131 << "G4DNAGillespieDirectMethod::Propensi << 130 << ConfB->GetName() << "(" << typeBNumber
132 << ConfA->GetName() << "(" << typeANumbe << 131 << ") : propensity : " << value
133 << "(" << typeBNumber << ") : propensity << 132 << " GetObservedReactionRateConstant : "
134 << " GetObservedReactionRateConstant : " << 133 << data->GetObservedReactionRateConstant()
135 << data->GetObservedReactionRateConstant << 134 << " GetEffectiveReactionRadius : "
136 << " GetEffectiveReactionRadius : " << 135 << G4BestUnit(data->GetEffectiveReactionRadius(), "Length")
137 << G4BestUnit(data->GetEffectiveReaction << 136 << " k : " << k << " volume : " << VolumeOfNode(index)
138 << " k : " << k << " volume : " << Volum << 137 << " Index : " << index << G4endl;
139 G4Exception("G4DNAGillespieDirectMethod::P << 138 assert(false);
140 "G4DNAGillespieDirectMethod013 <<
141 exceptionDescription); <<
142 } 139 }
143 140
144 #ifdef DEBUG 141 #ifdef DEBUG
145 if(value > 0) 142 if(value > 0)
146 G4cout << "G4DNAGillespieDirectMethod::Pro 143 G4cout << "G4DNAGillespieDirectMethod::PropensityFunction for : "
147 << ConfA->GetName() << "(" << typeA 144 << ConfA->GetName() << "(" << typeANumber << ") + "
148 << ConfB->GetName() << "(" << typeB 145 << ConfB->GetName() << "(" << typeBNumber
149 << ") : propensity : " << value 146 << ") : propensity : " << value
150 << " Time to Reaction : " << G4Bes 147 << " Time to Reaction : " << G4BestUnit(timeToReaction, "Time")
151 << " k : " << k << " Index : " << i 148 << " k : " << k << " Index : " << index << G4endl;
152 #endif 149 #endif
153 150
154 return value; 151 return value;
155 } 152 }
156 153
157 void G4DNAGillespieDirectMethod::Initialize() 154 void G4DNAGillespieDirectMethod::Initialize()
158 { 155 {
159 // for Scavenger 156 // for Scavenger
160 fpScavengerMaterial = dynamic_cast<G4DNAScav 157 fpScavengerMaterial = dynamic_cast<G4DNAScavengerMaterial*>(
161 G4Scheduler::Instance()->GetScavengerMater 158 G4Scheduler::Instance()->GetScavengerMaterial());
162 fEquilibriumProcesses.emplace( <<
163 std::make_pair(6, std::make_unique<G4ChemE <<
164 fEquilibriumProcesses.emplace( <<
165 std::make_pair(7, std::make_unique<G4ChemE <<
166 fEquilibriumProcesses.emplace( <<
167 std::make_pair(8, std::make_unique<G4ChemE <<
168 for(auto& it : fEquilibriumProcesses) <<
169 { <<
170 it.second->Initialize(); <<
171 it.second->SetVerbose(fVerbose); <<
172 } <<
173 } <<
174 159
175 void G4DNAGillespieDirectMethod::CreateEvents( <<
176 { <<
177 auto begin = fpMesh->begin(); 160 auto begin = fpMesh->begin();
178 auto end = fpMesh->end(); 161 auto end = fpMesh->end();
179 for(; begin != end; begin++) 162 for(; begin != end; begin++)
180 { 163 {
181 auto index = std::get<0>(*begin); << 164 auto key = begin->first;
182 #ifdef DEBUG 165 #ifdef DEBUG
183 fpMesh->PrintVoxel(index); << 166 fpMesh->PrintVoxel(fpMesh->GetIndex(key));
184 #endif 167 #endif
185 CreateEvent(index); << 168 CreateEvent(key);
186 } 169 }
187 } 170 }
188 171
189 void G4DNAGillespieDirectMethod::SetTimeStep(c 172 void G4DNAGillespieDirectMethod::SetTimeStep(const G4double& stepTime)
190 { 173 {
191 fTimeStep = stepTime; 174 fTimeStep = stepTime;
192 } 175 }
193 void G4DNAGillespieDirectMethod::CreateEvent(c << 176 void G4DNAGillespieDirectMethod::CreateEvent(unsigned int key)
194 { 177 {
195 const auto& voxel = fpMesh->GetVoxel(index); <<
196 if(std::get<2>(voxel).empty()) <<
197 { <<
198 G4ExceptionDescription exceptionDescriptio <<
199 exceptionDescription << "This voxel : " << <<
200 << " is not ready to <<
201 G4Exception("G4DNAGillespieDirectMethod::C <<
202 "G4DNAGillespieDirectMethod05" <<
203 exceptionDescription); <<
204 } <<
205 G4double r1 = G4UniformRand(); 178 G4double r1 = G4UniformRand();
206 G4double r2 = G4UniformRand(); 179 G4double r2 = G4UniformRand();
207 G4double dAlpha0 = DiffusiveJumping(voxel << 180 auto index = fpMesh->GetIndex(key);
208 G4double rAlpha0 = Reaction(voxel); << 181 G4double dAlpha0 = DiffusiveJumping(index);
>> 182 G4double rAlpha0 = Reaction(index);
209 G4double alphaTotal = dAlpha0 + rAlpha0; 183 G4double alphaTotal = dAlpha0 + rAlpha0;
210 184
211 if(alphaTotal == 0) 185 if(alphaTotal == 0)
212 { 186 {
213 return; 187 return;
214 } 188 }
215 auto timeStep = ((1.0 / (alphaTotal)) * std: 189 auto timeStep = ((1.0 / (alphaTotal)) * std::log(1.0 / r1)) + fTimeStep;
216 190
217 #ifdef DEBUG 191 #ifdef DEBUG
218 G4cout << "r2 : " << r2 << " rAlpha0 : " << 192 G4cout << "r2 : " << r2 << " rAlpha0 : " << rAlpha0
219 << " dAlpha0 : " << dAlpha0 << " r 193 << " dAlpha0 : " << dAlpha0 << " rAlpha0 / (dAlpha0 + rAlpha0) : "
220 << rAlpha0 / (dAlpha0 + rAlpha0) << G 194 << rAlpha0 / (dAlpha0 + rAlpha0) << G4endl;
221 #endif 195 #endif
222 if(r2 < rAlpha0 / alphaTotal) 196 if(r2 < rAlpha0 / alphaTotal)
223 { 197 {
224 if(fVerbose > 1) 198 if(fVerbose > 1)
225 { 199 {
226 G4cout << "=>>>>reaction at : " << timeS 200 G4cout << "=>>>>reaction at : " << timeStep << " timeStep : "
227 << G4BestUnit(((1.0 / alphaTotal) 201 << G4BestUnit(((1.0 / alphaTotal) * std::log(1.0 / r1)), "Time")
228 << G4endl; 202 << G4endl;
229 } 203 }
230 auto rSelectedIter = fReactionDataMap.uppe 204 auto rSelectedIter = fReactionDataMap.upper_bound(r2 * alphaTotal);
231 fpEventSet->CreateEvent(timeStep, index, r << 205 fpEventSet->CreateEvent(timeStep, key, rSelectedIter->second);
232 } 206 }
233 else if(dAlpha0 > 0) 207 else if(dAlpha0 > 0)
234 { 208 {
235 if(fVerbose > 1) 209 if(fVerbose > 1)
236 { 210 {
237 G4cout << "=>>>>jumping at : " << timeSt 211 G4cout << "=>>>>jumping at : " << timeStep << " timeStep : "
238 << G4BestUnit(((1.0 / alphaTotal) 212 << G4BestUnit(((1.0 / alphaTotal) * std::log(1.0 / r1)), "Time")
239 << G4endl; 213 << G4endl;
240 } 214 }
241 215
242 auto dSelectedIter = fJumpingDataMap.upper 216 auto dSelectedIter = fJumpingDataMap.upper_bound(r2 * alphaTotal - rAlpha0);
243 auto pDSelected = 217 auto pDSelected =
244 std::make_unique<std::pair<MolType, Inde 218 std::make_unique<std::pair<MolType, Index>>(dSelectedIter->second);
245 fpEventSet->CreateEvent(timeStep, index, s << 219 fpEventSet->CreateEvent(timeStep, key, std::move(pDSelected));
246 } 220 }
247 #ifdef DEBUG 221 #ifdef DEBUG
248 G4cout << G4endl; 222 G4cout << G4endl;
249 #endif 223 #endif
250 } 224 }
251 225
252 G4double G4DNAGillespieDirectMethod::Reaction( << 226 G4double G4DNAGillespieDirectMethod::Reaction(const Index& index)
253 { 227 {
254 fReactionDataMap.clear(); 228 fReactionDataMap.clear();
255 G4double alpha0 = 0; 229 G4double alpha0 = 0;
256 const auto& dataList = << 230 auto dataList = fMolecularReactions->GetVectorOfReactionData();
257 fMolecularReactions->GetVectorOfReactionDa <<
258 if(dataList.empty()) 231 if(dataList.empty())
259 { 232 {
260 G4ExceptionDescription exceptionDescriptio << 233 G4cout << "MolecularReactionTable empty" << G4endl;
261 exceptionDescription << "MolecularReaction << 234 assert(false);
262 G4Exception("G4DNAGillespieDirectMethod::R <<
263 "G4DNAGillespieDirectMethod01" <<
264 exceptionDescription); <<
265 } 235 }
266 <<
267 for(const auto& it : dataList) 236 for(const auto& it : dataList)
268 { 237 {
269 if(!IsEquilibrium(it->GetReactionType())) << 238 auto propensity = PropensityFunction(index, it);
270 { <<
271 continue; <<
272 } <<
273 auto propensity = PropensityFunction(voxel <<
274 if(propensity == 0) 239 if(propensity == 0)
275 { 240 {
276 continue; 241 continue;
277 } 242 }
278 alpha0 += propensity; 243 alpha0 += propensity;
279 fReactionDataMap[alpha0] = it; 244 fReactionDataMap[alpha0] = it;
280 } 245 }
281 #ifdef DEBUG 246 #ifdef DEBUG
282 G4cout << "Reaction :alpha0 : " << alpha0 < 247 G4cout << "Reaction :alpha0 : " << alpha0 << G4endl;
283 #endif 248 #endif
284 return alpha0; 249 return alpha0;
285 } 250 }
286 251
287 G4double G4DNAGillespieDirectMethod::Diffusive << 252 G4double G4DNAGillespieDirectMethod::DiffusiveJumping(const Index& index)
288 { 253 {
289 fJumpingDataMap.clear(); 254 fJumpingDataMap.clear();
290 G4double alpha0 = 0; 255 G4double alpha0 = 0;
291 auto index = std::get<0>(voxel); <<
292 auto NeighboringVoxels = fpMesh->FindNeighbo 256 auto NeighboringVoxels = fpMesh->FindNeighboringVoxels(index);
293 if(NeighboringVoxels.empty()) 257 if(NeighboringVoxels.empty())
294 { 258 {
295 return 0; 259 return 0;
296 } 260 }
297 auto iter = G4MoleculeTable::Instance()->Get 261 auto iter = G4MoleculeTable::Instance()->GetConfigurationIterator();
298 while(iter()) 262 while(iter())
299 { 263 {
300 const auto* conf = iter.value(); << 264 const auto conf = iter.value();
301 auto propensity = PropensityFunction(voxe << 265 auto propensity = PropensityFunction(index, conf);
302 if(propensity == 0) 266 if(propensity == 0)
303 { 267 {
304 continue; 268 continue;
305 } 269 }
306 for(const auto& it_Neighbor : NeighboringV 270 for(const auto& it_Neighbor : NeighboringVoxels)
307 { 271 {
308 alpha0 += propensity; 272 alpha0 += propensity;
309 fJumpingDataMap[alpha0] = std::make_pair 273 fJumpingDataMap[alpha0] = std::make_pair(conf, it_Neighbor);
310 #ifdef DEBUG 274 #ifdef DEBUG
311 G4cout << "mole : " << conf->GetName() 275 G4cout << "mole : " << conf->GetName()
312 << " number : " << ComputeNumberI 276 << " number : " << ComputeNumberInNode(index, conf)
313 << " propensity : " << propensity 277 << " propensity : " << propensity << " alpha0 : " << alpha0
314 << G4endl; 278 << G4endl;
315 #endif 279 #endif
316 } 280 }
317 } 281 }
318 #ifdef DEBUG 282 #ifdef DEBUG
319 G4cout << "DiffusiveJumping :alpha0 : " << 283 G4cout << "DiffusiveJumping :alpha0 : " << alpha0 << G4endl;
320 #endif 284 #endif
321 return alpha0; 285 return alpha0;
322 } 286 }
323 287
324 G4double G4DNAGillespieDirectMethod::ComputeNu 288 G4double G4DNAGillespieDirectMethod::ComputeNumberInNode(
325 const Voxel& voxel, MolType type) // depend << 289 const Index& index, MolType type) // depend node ?
326 { 290 {
327 if(type->GetDiffusionCoefficient() != 0) 291 if(type->GetDiffusionCoefficient() != 0)
328 { 292 {
329 const auto& node = std::get<2>(voxel); << 293 const auto& node = fpMesh->GetVoxelMapList(index);
330 const auto& it = node.find(type); 294 const auto& it = node.find(type);
331 return (it != node.end()) ? (it->second) : 295 return (it != node.end()) ? (it->second) : 0;
332 } 296 }
333 return 0; << 297 else
>> 298 {
>> 299 return 0;
>> 300 }
334 } 301 }
335 302
336 G4bool G4DNAGillespieDirectMethod::FindScaveng << 303 G4bool G4DNAGillespieDirectMethod::FindScavenging(const Index& index,
337 304 MolType moletype,
338 305 G4double& numberOfScavenger)
339 { 306 {
340 numberOfScavenger = 0; 307 numberOfScavenger = 0;
341 if(fpScavengerMaterial == nullptr) 308 if(fpScavengerMaterial == nullptr)
342 { 309 {
343 return false; 310 return false;
344 } 311 }
345 auto volumeOfNode = VolumeOfNode(voxel); << 312 auto volumeOfNode = VolumeOfNode(index);
346 if(G4MoleculeTable::Instance()->GetConfigura 313 if(G4MoleculeTable::Instance()->GetConfiguration("H2O") == moletype)
347 { 314 {
348 auto factor = Avogadro * volumeOfNod 315 auto factor = Avogadro * volumeOfNode;
349 numberOfScavenger = factor; 316 numberOfScavenger = factor;
350 return true; 317 return true;
351 } 318 }
352 319
353 G4double totalNumber = 320 G4double totalNumber =
354 fpScavengerMaterial->GetNumberMoleculePerV 321 fpScavengerMaterial->GetNumberMoleculePerVolumeUnitForMaterialConf(
355 moletype); 322 moletype);
356 if(totalNumber == 0) 323 if(totalNumber == 0)
357 { 324 {
358 return false; 325 return false;
359 } 326 }
360 << 327 else
361 G4double numberInDouble = volumeOfNode * std <<
362 fpMesh->GetBoundin <<
363 auto numberInInterg = (int64_t) (std::floor( <<
364 G4double change = numberInDouble - numbe <<
365 G4UniformRand() > change ? numberOfScavenger <<
366 : numberOfScave <<
367 return true; <<
368 } <<
369 <<
370 G4bool G4DNAGillespieDirectMethod::IsEquilibri <<
371 { <<
372 auto reaction = fEquilibriumProcesses.find(r <<
373 if(reaction == fEquilibriumProcesses.end()) <<
374 { 328 {
>> 329 G4double numberInDouble =
>> 330 volumeOfNode * std::floor(totalNumber) / fpMesh->GetBoundingBox().Volume();
>> 331 auto numberInInterg = (int) (std::floor(numberInDouble));
>> 332 G4double ram = G4UniformRand();
>> 333 G4double change = numberInDouble - numberInInterg;
>> 334 if(ram > change)
>> 335 {
>> 336 numberOfScavenger = numberInInterg;
>> 337 }
>> 338 else
>> 339 {
>> 340 numberOfScavenger = numberInInterg + 1;
>> 341 }
375 return true; 342 return true;
376 }else <<
377 { <<
378 return (reaction->second->GetEquilibriumSt <<
379 } 343 }
380 <<
381 } 344 }
382 <<
383 G4bool G4DNAGillespieDirectMethod::SetEquilibr <<
384 { <<
385 for(auto& it : fEquilibriumProcesses) <<
386 { <<
387 it.second->SetGlobalTime(fTimeStep); <<
388 it.second->SetEquilibrium(pReaction); <<
389 if(it.second->IsStatusChanged()) return tr <<
390 } <<
391 return false; <<
392 } <<
393 <<
394 void G4DNAGillespieDirectMethod::ResetEquilibr <<
395 { <<
396 for(auto& it : fEquilibriumProcesses) <<
397 { <<
398 it.second->Reset(); <<
399 } <<
400 } <<
401 <<
402 <<
403 <<
404 345