| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 #include <memory> 27 #include <memory>
28 #include "G4DNAEventScheduler.hh" 28 #include "G4DNAEventScheduler.hh"
29 #include "G4DNAGillespieDirectMethod.hh" 29 #include "G4DNAGillespieDirectMethod.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 #include "G4UnitsTable.hh" 31 #include "G4UnitsTable.hh"
32 #include "G4DNAUpdateSystemModel.hh" 32 #include "G4DNAUpdateSystemModel.hh"
33 #include "G4DNAMolecularReactionTable.hh" 33 #include "G4DNAMolecularReactionTable.hh"
34 #include "G4Timer.hh" 34 #include "G4Timer.hh"
35 #include "G4Scheduler.hh" 35 #include "G4Scheduler.hh"
36 #include "G4UserMeshAction.hh" 36 #include "G4UserMeshAction.hh"
37 #include "G4MoleculeCounter.hh" 37 #include "G4MoleculeCounter.hh"
38 #include "G4DNAScavengerMaterial.hh" 38 #include "G4DNAScavengerMaterial.hh"
39 #include "G4Molecule.hh" 39 #include "G4Molecule.hh"
40 40
41 G4DNAEventScheduler::G4DNAEventScheduler() << 41 G4DNAEventScheduler::G4DNAEventScheduler(const G4DNABoundingBox& boundingBox,
42 : fpGillespieReaction(new G4DNAGillespieDire << 42 G4int pixel)
>> 43 : IEventScheduler()
>> 44 , fPixel(pixel)
>> 45 , fInitialPixels(fPixel)
>> 46 , fpMesh(new G4DNAMesh(boundingBox, fPixel))
>> 47 , fpGillespieReaction(new G4DNAGillespieDirectMethod())
43 , fpEventSet(new G4DNAEventSet()) 48 , fpEventSet(new G4DNAEventSet())
44 , fpUpdateSystem(new G4DNAUpdateSystemModel( 49 , fpUpdateSystem(new G4DNAUpdateSystemModel())
45 {} << 50 {
>> 51 if(!CheckingReactionRadius(fpMesh->GetResolution()))
>> 52 {
>> 53 G4String WarMessage = "resolution is not good : " +
>> 54 std::to_string(fpMesh->GetResolution() / nm);
>> 55 G4Exception("G4DNAEventScheduler::InitializeInMesh()", "WrongResolution",
>> 56 JustWarning, WarMessage);
>> 57 }
>> 58 }
46 59
47 void G4DNAEventScheduler::ClearAndReChargeCoun 60 void G4DNAEventScheduler::ClearAndReChargeCounter()
48 { 61 {
49 fCounterMap.clear(); 62 fCounterMap.clear();
50 if(fTimeToRecord.empty()) 63 if(fTimeToRecord.empty())
51 { 64 {
52 G4String WarMessage = "fTimeToRecord is em 65 G4String WarMessage = "fTimeToRecord is empty ";
53 G4Exception("G4DNAEventScheduler::ClearAnd 66 G4Exception("G4DNAEventScheduler::ClearAndReChargeCounter()",
54 "TimeToRecord is empty", JustW 67 "TimeToRecord is empty", JustWarning, WarMessage);
55 } 68 }
56 fLastRecoredTime = fTimeToRecord.begin(); 69 fLastRecoredTime = fTimeToRecord.begin();
57 70
58 if(G4VMoleculeCounter::Instance()->InUse()) 71 if(G4VMoleculeCounter::Instance()->InUse()) // copy from MoleculeCounter
59 { 72 {
60 G4MoleculeCounter::RecordedMolecules speci 73 G4MoleculeCounter::RecordedMolecules species;
61 species = G4MoleculeCounter::Instance()->G 74 species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
62 if(species.get() == nullptr) 75 if(species.get() == nullptr)
63 { 76 {
64 return; 77 return;
65 } 78 }
66 if(species->empty()) << 79 else if(species->empty())
67 { 80 {
68 G4MoleculeCounter::Instance()->ResetCoun 81 G4MoleculeCounter::Instance()->ResetCounter();
69 return; 82 return;
70 } 83 }
71 for(auto time_mol : fTimeToRecord) 84 for(auto time_mol : fTimeToRecord)
72 { 85 {
73 if(time_mol > fStartTime) 86 if(time_mol > fStartTime)
74 { 87 {
75 continue; 88 continue;
76 } 89 }
77 90
78 for(auto molecule : *species) 91 for(auto molecule : *species)
79 { 92 {
80 G4int n_mol = G4MoleculeCounter::Insta 93 G4int n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(
81 molecule, time_mol); 94 molecule, time_mol);
82 95
83 if(n_mol < 0) 96 if(n_mol < 0)
84 { 97 {
85 G4cerr << "G4DNAEventScheduler::Clea 98 G4cerr << "G4DNAEventScheduler::ClearAndReChargeCounter() ::N "
86 "molecules not valid < 0 " 99 "molecules not valid < 0 "
87 << G4endl; 100 << G4endl;
88 G4Exception("", "N<0", FatalExceptio 101 G4Exception("", "N<0", FatalException, "");
89 } 102 }
90 fCounterMap[time_mol][molecule] = n_mo 103 fCounterMap[time_mol][molecule] = n_mol;
91 } 104 }
92 fLastRecoredTime++; 105 fLastRecoredTime++;
93 } 106 }
94 G4MoleculeCounter::Instance()->ResetCounte 107 G4MoleculeCounter::Instance()->ResetCounter(); // reset
95 G4MoleculeCounter::Instance()->Use(false); 108 G4MoleculeCounter::Instance()->Use(false); // no more used
96 } 109 }
97 else 110 else
98 { 111 {
99 G4ExceptionDescription exceptionDescriptio 112 G4ExceptionDescription exceptionDescription;
100 exceptionDescription << "G4VMoleculeCounte 113 exceptionDescription << "G4VMoleculeCounter is not used";
101 G4Exception("G4DNAEventScheduler::ClearAnd 114 G4Exception("G4DNAEventScheduler::ClearAndReChargeCounter()",
102 "G4DNAEventScheduler010", Just 115 "G4DNAEventScheduler010", JustWarning, exceptionDescription);
103 } 116 }
104 } 117 }
105 118
106 [[maybe_unused]] void G4DNAEventScheduler::Add 119 [[maybe_unused]] void G4DNAEventScheduler::AddTimeToRecord(const G4double& time)
107 { 120 {
108 if(fTimeToRecord.find(time) == fTimeToRecord 121 if(fTimeToRecord.find(time) == fTimeToRecord.end())
109 { 122 {
110 fTimeToRecord.insert(time); 123 fTimeToRecord.insert(time);
111 } 124 }
112 } 125 }
113 126
114 G4DNAEventScheduler::~G4DNAEventScheduler() = 127 G4DNAEventScheduler::~G4DNAEventScheduler() = default;
115 128
116 void G4DNAEventScheduler::Voxelizing() 129 void G4DNAEventScheduler::Voxelizing()
117 { 130 {
118 auto pMainList = G4ITTrackHolder::Instance() 131 auto pMainList = G4ITTrackHolder::Instance()->GetMainList();
119 std::map<G4VDNAMesh::Index, MapList> TrackKe 132 std::map<G4VDNAMesh::Index, MapList> TrackKeyMap;
120 for(auto track : *pMainList) 133 for(auto track : *pMainList)
121 { 134 {
122 auto molType = GetMolecule(track)->GetMole 135 auto molType = GetMolecule(track)->GetMolecularConfiguration();
123 136
124 auto pScavengerMaterial = dynamic_cast<G4D 137 auto pScavengerMaterial = dynamic_cast<G4DNAScavengerMaterial*>(
125 G4Scheduler::Instance()->GetScavengerMat 138 G4Scheduler::Instance()->GetScavengerMaterial());
126 if(pScavengerMaterial != nullptr && 139 if(pScavengerMaterial != nullptr &&
127 pScavengerMaterial->find(molType)) // 140 pScavengerMaterial->find(molType)) // avoid voxelize the scavenger
128 { 141 {
129 continue; 142 continue;
130 } 143 }
131 144
132 auto key = fpMesh->GetIndex(track->GetPosi 145 auto key = fpMesh->GetIndex(track->GetPosition());
133 if(TrackKeyMap.find(key) != TrackKeyMap.en 146 if(TrackKeyMap.find(key) != TrackKeyMap.end())
134 { 147 {
135 std::map<MolType, size_t>& TrackTypeMap 148 std::map<MolType, size_t>& TrackTypeMap = TrackKeyMap[key];
136 if(TrackTypeMap.find(molType) != TrackTy 149 if(TrackTypeMap.find(molType) != TrackTypeMap.end())
137 { 150 {
138 TrackTypeMap[molType]++; 151 TrackTypeMap[molType]++;
139 } 152 }
140 else 153 else
141 { 154 {
142 TrackTypeMap[molType] = 1; 155 TrackTypeMap[molType] = 1;
143 } 156 }
144 } 157 }
145 else 158 else
146 { 159 {
147 TrackKeyMap[key][molType] = 1; 160 TrackKeyMap[key][molType] = 1;
148 } 161 }
149 } 162 }
150 163
151 for(auto& it : TrackKeyMap) 164 for(auto& it : TrackKeyMap)
152 { 165 {
153 fpMesh->InitializeVoxel(it.first, std::mov 166 fpMesh->InitializeVoxel(it.first, std::move(it.second));
154 } 167 }
155 } 168 }
156 169
157 void G4DNAEventScheduler::ReVoxelizing(G4int p 170 void G4DNAEventScheduler::ReVoxelizing(G4int pixel)
158 { 171 {
159 fPixel = pixel; 172 fPixel = pixel;
160 auto newMesh = new G4DNAMesh(fpMesh->GetBoun 173 auto newMesh = new G4DNAMesh(fpMesh->GetBoundingBox(), fPixel);
161 174
162 auto begin = fpMesh->begin(); 175 auto begin = fpMesh->begin();
163 auto end = fpMesh->end(); 176 auto end = fpMesh->end();
164 std::map<G4VDNAMesh::Index, MapList> TrackKe 177 std::map<G4VDNAMesh::Index, MapList> TrackKeyMap;
165 for(; begin != end; begin++) 178 for(; begin != end; begin++)
166 { 179 {
167 auto index = std::get<0>(*begin); 180 auto index = std::get<0>(*begin);
168 auto newIndex = fpMesh->ConvertIndex(index 181 auto newIndex = fpMesh->ConvertIndex(index, fPixel);
169 if(TrackKeyMap.find(newIndex) == TrackKeyM 182 if(TrackKeyMap.find(newIndex) == TrackKeyMap.end())
170 { 183 {
171 TrackKeyMap[newIndex] = std::get<2>(*beg 184 TrackKeyMap[newIndex] = std::get<2>(*begin);
172 } 185 }
173 else 186 else
174 { 187 {
175 for(const auto& it : std::get<2>(*begin) 188 for(const auto& it : std::get<2>(*begin))
176 { 189 {
177 TrackKeyMap[newIndex][it.first] += it. 190 TrackKeyMap[newIndex][it.first] += it.second;
178 } 191 }
179 if(fVerbose > 1) 192 if(fVerbose > 1)
180 { 193 {
181 G4cout << " ReVoxelizing:: Old index : 194 G4cout << " ReVoxelizing:: Old index : " << index
182 << " new index : " << fpMesh->C 195 << " new index : " << fpMesh->ConvertIndex(index, fPixel)
183 << " number: " << std::get<2>(* 196 << " number: " << std::get<2>(*begin).begin()->second << G4endl;
184 } 197 }
185 } 198 }
186 } 199 }
187 fpMesh.reset(newMesh); 200 fpMesh.reset(newMesh);
188 201
189 for(auto& it : TrackKeyMap) 202 for(auto& it : TrackKeyMap)
190 { 203 {
191 fpMesh->InitializeVoxel(it.first, std::mov 204 fpMesh->InitializeVoxel(it.first, std::move(it.second));
192 } 205 }
193 } 206 }
194 void G4DNAEventScheduler::Reset() 207 void G4DNAEventScheduler::Reset()
195 { 208 {
196 // find another solultion 209 // find another solultion
197 fGlobalTime = fEndTime; 210 fGlobalTime = fEndTime;
198 211
199 // 212 //
200 LastRegisterForCounter();//Last register for << 213 // RecordTime();//Last register for counter
201 214
202 if(fVerbose > 0) 215 if(fVerbose > 0)
203 { 216 {
204 G4cout << "End Processing and reset Gird, 217 G4cout << "End Processing and reset Gird, ScavengerTable, EventSet for new "
205 "simulation!!!!" 218 "simulation!!!!"
206 << G4endl; 219 << G4endl;
207 } 220 }
208 fInitialized = false; 221 fInitialized = false;
209 fTimeStep = 0; 222 fTimeStep = 0;
210 fStepNumber = 0; 223 fStepNumber = 0;
211 fGlobalTime = fStartTime; 224 fGlobalTime = fStartTime;
212 fRunning = true; 225 fRunning = true;
213 fReactionNumber = 0; 226 fReactionNumber = 0;
214 fJumpingNumber = 0; 227 fJumpingNumber = 0;
215 228
216 fpEventSet->RemoveEventSet(); 229 fpEventSet->RemoveEventSet();
217 fpMesh->Reset(); 230 fpMesh->Reset();
218 fpGillespieReaction->ResetEquilibrium(); <<
219 } 231 }
220 232
221 void G4DNAEventScheduler::Initialize(const G4D << 233 void G4DNAEventScheduler::Initialize()
222 G4int pix <<
223 { 234 {
224 if(!fInitialized) 235 if(!fInitialized)
225 { 236 {
226 fPixel = pixel; << 237 fPixel = fInitialPixels;
227 fpMesh = std::make_unique<G4DNAMesh>(bound << 238 fpMesh = std::make_unique<G4DNAMesh>(fpMesh->GetBoundingBox(), fPixel);
228 <<
229 if(!CheckingReactionRadius(fpMesh->GetReso <<
230 { <<
231 G4String WarMessage = "resolution is not <<
232 std::to_string(fpM <<
233 G4Exception("G4DNAEventScheduler::Initia <<
234 JustWarning, WarMessage); <<
235 } <<
236 239
237 // Scavenger(); 240 // Scavenger();
238 241
239 auto pScavengerMaterial = dynamic_cast<G4D 242 auto pScavengerMaterial = dynamic_cast<G4DNAScavengerMaterial*>(
240 G4Scheduler::Instance()->GetScavengerMat 243 G4Scheduler::Instance()->GetScavengerMaterial());
241 if(pScavengerMaterial == nullptr) 244 if(pScavengerMaterial == nullptr)
242 { 245 {
243 G4cout << "There is no scavenger" << G4e 246 G4cout << "There is no scavenger" << G4endl;
244 } 247 }
245 else 248 else
246 { 249 {
247 if(fVerbose > 1) 250 if(fVerbose > 1)
248 { 251 {
249 pScavengerMaterial->PrintInfo(); 252 pScavengerMaterial->PrintInfo();
250 } 253 }
251 } 254 }
252 255
253 Voxelizing(); 256 Voxelizing();
254 fpGillespieReaction->SetVoxelMesh(*fpMesh) 257 fpGillespieReaction->SetVoxelMesh(*fpMesh);
255 fpGillespieReaction->SetEventSet(fpEventSe 258 fpGillespieReaction->SetEventSet(fpEventSet.get());
256 fpGillespieReaction->SetTimeStep(0);// res << 259 fpGillespieReaction->SetTimeStep(
>> 260 0); // reset fTimeStep = 0 in fpGillespieReaction
257 fpGillespieReaction->Initialize(); 261 fpGillespieReaction->Initialize();
258 fpGillespieReaction->CreateEvents(); <<
259 fpUpdateSystem->SetMesh(fpMesh.get()); 262 fpUpdateSystem->SetMesh(fpMesh.get());
260 ClearAndReChargeCounter(); 263 ClearAndReChargeCounter();
261 fInitialized = true; 264 fInitialized = true;
262 } 265 }
263 266
264 if(fVerbose > 0) 267 if(fVerbose > 0)
265 { 268 {
266 fpUpdateSystem->SetVerbose(1); 269 fpUpdateSystem->SetVerbose(1);
267 } 270 }
268 271
269 if(fVerbose > 2) 272 if(fVerbose > 2)
270 { 273 {
271 fpMesh->PrintMesh(); 274 fpMesh->PrintMesh();
272 } 275 }
273 } 276 }
274 void G4DNAEventScheduler::InitializeInMesh() 277 void G4DNAEventScheduler::InitializeInMesh()
275 { 278 {
276 if(fPixel <= 1) 279 if(fPixel <= 1)
277 { 280 {
278 fRunning = false; 281 fRunning = false;
279 return; 282 return;
280 } 283 }
281 // TEST /3 284 // TEST /3
282 ReVoxelizing(fPixel / 2); // 285 ReVoxelizing(fPixel / 2); //
283 // ReVoxelizing(fPixel/3);// 286 // ReVoxelizing(fPixel/3);//
284 287
285 fpGillespieReaction->SetVoxelMesh(*fpMesh); 288 fpGillespieReaction->SetVoxelMesh(*fpMesh);
286 fpUpdateSystem->SetMesh(fpMesh.get()); 289 fpUpdateSystem->SetMesh(fpMesh.get());
287 fpGillespieReaction->CreateEvents(); << 290 fpGillespieReaction->Initialize();
288 } 291 }
289 292
290 void G4DNAEventScheduler::ResetInMesh() 293 void G4DNAEventScheduler::ResetInMesh()
291 { 294 {
292 if(fVerbose > 0) 295 if(fVerbose > 0)
293 { 296 {
294 G4cout 297 G4cout
295 << "*** End Processing In Mesh and reset 298 << "*** End Processing In Mesh and reset Mesh, EventSet for new Mesh!!!!"
296 << G4endl; 299 << G4endl;
297 } 300 }
298 fpEventSet->RemoveEventSet(); 301 fpEventSet->RemoveEventSet();
299 fInitialized = false; 302 fInitialized = false;
300 fIsChangeMesh = false; 303 fIsChangeMesh = false;
301 fReactionNumber = 0; 304 fReactionNumber = 0;
302 fJumpingNumber = 0; 305 fJumpingNumber = 0;
303 fStepNumberInMesh = 0; 306 fStepNumberInMesh = 0;
304 } 307 }
305 308
306 G4double G4DNAEventScheduler::GetStartTime() c 309 G4double G4DNAEventScheduler::GetStartTime() const { return fStartTime; }
307 310
308 G4double G4DNAEventScheduler::GetEndTime() con 311 G4double G4DNAEventScheduler::GetEndTime() const { return fEndTime; }
309 312
310 [[maybe_unused]] G4double G4DNAEventScheduler: 313 [[maybe_unused]] G4double G4DNAEventScheduler::GetTimeStep() const
311 { 314 {
312 return fTimeStep; 315 return fTimeStep;
313 } 316 }
314 317
315 G4int G4DNAEventScheduler::GetVerbose() const 318 G4int G4DNAEventScheduler::GetVerbose() const { return fVerbose; }
316 319
317 [[maybe_unused]] void G4DNAEventScheduler::Set 320 [[maybe_unused]] void G4DNAEventScheduler::SetMaxNbSteps(G4int max)
318 { 321 {
319 fMaxStep = max; 322 fMaxStep = max;
320 } 323 }
321 324
322 [[maybe_unused]] void G4DNAEventScheduler::Set 325 [[maybe_unused]] void G4DNAEventScheduler::SetStartTime(G4double time)
323 { 326 {
324 fStartTime = time; 327 fStartTime = time;
325 fGlobalTime = fStartTime; 328 fGlobalTime = fStartTime;
326 } 329 }
327 330
328 void G4DNAEventScheduler::Stop() { fRunning = 331 void G4DNAEventScheduler::Stop() { fRunning = false; }
329 void G4DNAEventScheduler::Run() 332 void G4DNAEventScheduler::Run()
330 { 333 {
331 G4Timer localtimer; 334 G4Timer localtimer;
332 if(fVerbose > 2) 335 if(fVerbose > 2)
333 { 336 {
334 localtimer.Start(); 337 localtimer.Start();
335 G4cout << "***G4DNAEventScheduler::Run*** 338 G4cout << "***G4DNAEventScheduler::Run*** for Pixel : " << fPixel << G4endl;
336 } 339 }
337 while(fEndTime > fGlobalTime && fRunning) 340 while(fEndTime > fGlobalTime && fRunning)
338 { 341 {
339 RunInMesh(); 342 RunInMesh();
340 } 343 }
341 if(fVerbose > 2) 344 if(fVerbose > 2)
342 { 345 {
343 if(!fRunning) 346 if(!fRunning)
344 { 347 {
345 G4cout << " StepNumber(" << fStepNumber 348 G4cout << " StepNumber(" << fStepNumber << ") = MaxStep(" << fMaxStep
346 << ")" << G4endl; 349 << ")" << G4endl;
347 } 350 }
348 else if(fEndTime <= fGlobalTime) 351 else if(fEndTime <= fGlobalTime)
349 { 352 {
350 G4cout << " GlobalTime(" << fGlobalTime 353 G4cout << " GlobalTime(" << fGlobalTime << ") > EndTime(" << fEndTime
351 << ")" 354 << ")"
352 << " StepNumber : " << fStepNumbe 355 << " StepNumber : " << fStepNumber << G4endl;
353 } 356 }
354 localtimer.Stop(); 357 localtimer.Stop();
355 G4cout << "***G4DNAEventScheduler::Ending: 358 G4cout << "***G4DNAEventScheduler::Ending::"
356 << G4BestUnit(fGlobalTime, "Time") 359 << G4BestUnit(fGlobalTime, "Time")
357 << " Events left : " << fpEventSet- 360 << " Events left : " << fpEventSet->size() << G4endl;
358 if(fVerbose > 1) 361 if(fVerbose > 1)
359 { 362 {
360 fpMesh->PrintMesh(); 363 fpMesh->PrintMesh();
361 } 364 }
362 G4cout << " Computing Time : " << localtim 365 G4cout << " Computing Time : " << localtimer << G4endl;
363 } 366 }
364 Reset(); 367 Reset();
365 } 368 }
366 369
367 void G4DNAEventScheduler::RunInMesh() 370 void G4DNAEventScheduler::RunInMesh()
368 { 371 {
369 if(!fInitialized) 372 if(!fInitialized)
370 { 373 {
371 InitializeInMesh(); 374 InitializeInMesh();
372 } 375 }
373 if(fVerbose > 0) 376 if(fVerbose > 0)
374 { 377 {
375 G4double resolution = fpMesh->GetResolutio 378 G4double resolution = fpMesh->GetResolution();
376 G4cout << "At Time : " << std::setw(7) << 379 G4cout << "At Time : " << std::setw(7) << G4BestUnit(fGlobalTime, "Time")
377 << " the Mesh has " << fPixel << " 380 << " the Mesh has " << fPixel << " x " << fPixel << " x " << fPixel
378 << " voxels with Resolution " << G4 381 << " voxels with Resolution " << G4BestUnit(resolution, "Length")
379 << " during next " 382 << " during next "
380 << G4BestUnit(resolution * resoluti 383 << G4BestUnit(resolution * resolution * C / (6 * D), "Time")
381 << G4endl; 384 << G4endl;
382 } 385 }
383 386
384 if(fVerbose > 2) 387 if(fVerbose > 2)
385 { 388 {
386 fpMesh->PrintMesh(); 389 fpMesh->PrintMesh();
387 } 390 }
388 391
389 if(fpUserMeshAction != nullptr) 392 if(fpUserMeshAction != nullptr)
390 { 393 {
391 fpUserMeshAction->BeginOfMesh(fpMesh.get() 394 fpUserMeshAction->BeginOfMesh(fpMesh.get(), fGlobalTime);
392 } 395 }
393 396
394 // if diffusive jumping is avaiable, EventSe 397 // if diffusive jumping is avaiable, EventSet is never empty
395 398
396 while(!fpEventSet->Empty() && !fIsChangeMesh 399 while(!fpEventSet->Empty() && !fIsChangeMesh && fEndTime > fGlobalTime)
397 { 400 {
398 Stepping(); 401 Stepping();
399 fGlobalTime = fTimeStep + fStartTime; 402 fGlobalTime = fTimeStep + fStartTime;
400 403
401 if(fpUserMeshAction != nullptr) 404 if(fpUserMeshAction != nullptr)
402 { 405 {
403 fpUserMeshAction->InMesh(fpMesh.get(), f 406 fpUserMeshAction->InMesh(fpMesh.get(), fGlobalTime);
404 } 407 }
405 408
406 if(fVerbose > 2) 409 if(fVerbose > 2)
407 { 410 {
408 G4cout << "fGlobalTime : " << G4BestUnit 411 G4cout << "fGlobalTime : " << G4BestUnit(fGlobalTime, "Time")
409 << " fTimeStep : " << G4BestUnit( 412 << " fTimeStep : " << G4BestUnit(fTimeStep, "Time") << G4endl;
410 } 413 }
411 G4double resolution = fpMesh->GetResolutio 414 G4double resolution = fpMesh->GetResolution();
412 fTransferTime = resolution * resolut 415 fTransferTime = resolution * resolution * C / (6 * D);
413 if(fTransferTime == 0) 416 if(fTransferTime == 0)
414 { 417 {
415 G4ExceptionDescription exceptionDescript 418 G4ExceptionDescription exceptionDescription;
416 exceptionDescription << "fTransferTime = 419 exceptionDescription << "fTransferTime == 0";
417 G4Exception("G4DNAEventScheduler::RunInM 420 G4Exception("G4DNAEventScheduler::RunInMesh", "G4DNAEventScheduler001",
418 FatalErrorInArgument, except 421 FatalErrorInArgument, exceptionDescription);
419 } 422 }
420 if(fTransferTime < fTimeStep && 423 if(fTransferTime < fTimeStep &&
421 fPixel != 1) // dont change Mesh if fP 424 fPixel != 1) // dont change Mesh if fPixel == 1
422 { 425 {
423 if(fSetChangeMesh) 426 if(fSetChangeMesh)
424 { 427 {
425 if(fVerbose > 1) 428 if(fVerbose > 1)
426 { 429 {
427 G4cout << " Pixels : " << fPixel << 430 G4cout << " Pixels : " << fPixel << " resolution : "
428 << G4BestUnit(fpMesh->GetReso 431 << G4BestUnit(fpMesh->GetResolution(), "Length")
429 << " fStepNumberInMesh : " < 432 << " fStepNumberInMesh : " << fStepNumberInMesh
430 << " at fGlobalTime : " << G4 433 << " at fGlobalTime : " << G4BestUnit(fGlobalTime, "Time")
431 << " at fTimeStep : " << G4Be 434 << " at fTimeStep : " << G4BestUnit(fTimeStep, "Time")
432 << " fReactionNumber : " << 435 << " fReactionNumber : " << fReactionNumber
433 << " transferTime : " << G4Be 436 << " transferTime : " << G4BestUnit(fTransferTime, "Time")
434 << G4endl; 437 << G4endl;
435 } 438 }
436 fIsChangeMesh = true; 439 fIsChangeMesh = true;
437 } 440 }
438 } 441 }
439 } 442 }
440 443
441 if(fVerbose > 1) 444 if(fVerbose > 1)
442 { 445 {
443 G4cout << "***G4DNAEventScheduler::Ending: 446 G4cout << "***G4DNAEventScheduler::Ending::"
444 << G4BestUnit(fGlobalTime, "Time") 447 << G4BestUnit(fGlobalTime, "Time")
445 << " Event left : " << fpEventSet-> 448 << " Event left : " << fpEventSet->size() << G4endl;
446 G4cout << " Due to : "; 449 G4cout << " Due to : ";
447 if(fpEventSet->Empty()) 450 if(fpEventSet->Empty())
448 { 451 {
449 G4cout << "EventSet is Empty" << G4endl; 452 G4cout << "EventSet is Empty" << G4endl;
450 } 453 }
451 else if(fIsChangeMesh) 454 else if(fIsChangeMesh)
452 { 455 {
453 G4cout << "Changing Mesh from : " << fPi 456 G4cout << "Changing Mesh from : " << fPixel
454 << " pixels to : " << fPixel / 2 457 << " pixels to : " << fPixel / 2 << " pixels" << G4endl;
455 G4cout << "Info : ReactionNumber : " << 458 G4cout << "Info : ReactionNumber : " << fReactionNumber
456 << " JumpingNumber : " << fJump 459 << " JumpingNumber : " << fJumpingNumber << G4endl;
457 } 460 }
458 else if(fEndTime > fGlobalTime) 461 else if(fEndTime > fGlobalTime)
459 { 462 {
460 G4cout << " GlobalTime(" << fGlobalTime 463 G4cout << " GlobalTime(" << fGlobalTime << ") > EndTime(" << fEndTime
461 << ")" 464 << ")"
462 << " StepNumber : " << fStepNumbe 465 << " StepNumber : " << fStepNumber << G4endl;
463 } 466 }
464 if(fVerbose > 2) 467 if(fVerbose > 2)
465 { 468 {
466 fpMesh->PrintMesh(); 469 fpMesh->PrintMesh();
467 } 470 }
468 G4cout << G4endl; 471 G4cout << G4endl;
469 } 472 }
470 473
471 if(fpUserMeshAction != nullptr) 474 if(fpUserMeshAction != nullptr)
472 { 475 {
473 fpUserMeshAction->EndOfMesh(fpMesh.get(), 476 fpUserMeshAction->EndOfMesh(fpMesh.get(), fGlobalTime);
474 } 477 }
475 ResetInMesh(); 478 ResetInMesh();
476 } 479 }
477 480
478 void G4DNAEventScheduler::Stepping() // this 481 void G4DNAEventScheduler::Stepping() // this event loop
479 { 482 {
480 fStepNumber < fMaxStep ? fStepNumber++ : sta << 483 fStepNumber < fMaxStep ? fStepNumber++ : fRunning = false;
481 if(fpEventSet->size() > fpMesh->size()) 484 if(fpEventSet->size() > fpMesh->size())
482 { 485 {
483 G4ExceptionDescription exceptionDescriptio 486 G4ExceptionDescription exceptionDescription;
484 exceptionDescription 487 exceptionDescription
485 << "impossible that fpEventSet->size() > 488 << "impossible that fpEventSet->size() > fpMesh->size()";
486 G4Exception("G4DNAEventScheduler::Stepping 489 G4Exception("G4DNAEventScheduler::Stepping", "G4DNAEventScheduler002",
487 FatalErrorInArgument, exceptio 490 FatalErrorInArgument, exceptionDescription);
488 } 491 }
489 492
490 auto selected = fpEventSet->begin(); 493 auto selected = fpEventSet->begin();
491 auto index = (*selected)->GetIndex(); 494 auto index = (*selected)->GetIndex();
492 495
493 if(fVerbose > 1) 496 if(fVerbose > 1)
494 { 497 {
495 G4cout << "G4DNAEventScheduler::Stepping() 498 G4cout << "G4DNAEventScheduler::Stepping()*********************************"
496 "*******" 499 "*******"
497 << G4endl; 500 << G4endl;
498 (*selected)->PrintEvent(); 501 (*selected)->PrintEvent();
499 } 502 }
500 503
501 // get selected time step 504 // get selected time step
502 fTimeStep = (*selected)->GetTime(); 505 fTimeStep = (*selected)->GetTime();
503 506
504 // selected data 507 // selected data
505 auto pJumping = (*selected)->GetJumpingData 508 auto pJumping = (*selected)->GetJumpingData();
506 auto pReaction = (*selected)->GetReactionDat 509 auto pReaction = (*selected)->GetReactionData();
507 510
508 fpUpdateSystem->SetGlobalTime(fTimeStep + 511 fpUpdateSystem->SetGlobalTime(fTimeStep +
509 fStartTime); 512 fStartTime); // this is just for printing
510 fpGillespieReaction->SetTimeStep(fTimeStep); 513 fpGillespieReaction->SetTimeStep(fTimeStep);
511 if(pJumping == nullptr && pReaction != nullp 514 if(pJumping == nullptr && pReaction != nullptr)
512 { 515 {
513 fpUpdateSystem->UpdateSystem(index, *pReac 516 fpUpdateSystem->UpdateSystem(index, *pReaction);
514 fpEventSet->RemoveEvent(selected); 517 fpEventSet->RemoveEvent(selected);
515 <<
516 //hoang's exp <<
517 if(fpGillespieReaction->SetEquilibrium(pRe <<
518 { <<
519 ResetEventSet(); <<
520 } <<
521 //end Hoang's exp <<
522 <<
523 // create new event 518 // create new event
524 fpGillespieReaction->CreateEvent(index); 519 fpGillespieReaction->CreateEvent(index);
525 fReactionNumber++; 520 fReactionNumber++;
526 // recordTime in reaction 521 // recordTime in reaction
527 RecordTime(); 522 RecordTime();
528 } 523 }
529 else if(pJumping != nullptr && pReaction == 524 else if(pJumping != nullptr && pReaction == nullptr)
530 { 525 {
531 // dont change this 526 // dont change this
532 fpUpdateSystem->UpdateSystem(index, *pJump 527 fpUpdateSystem->UpdateSystem(index, *pJumping);
533 // save jumping Index before delete select 528 // save jumping Index before delete selected event
534 auto jumpingIndex = pJumping->second; 529 auto jumpingIndex = pJumping->second;
535 fpEventSet->RemoveEvent(selected); 530 fpEventSet->RemoveEvent(selected);
536 // create new event 531 // create new event
537 // should create Jumping before key 532 // should create Jumping before key
538 fpGillespieReaction->CreateEvent(jumpingIn 533 fpGillespieReaction->CreateEvent(jumpingIndex);
539 fpGillespieReaction->CreateEvent(index); 534 fpGillespieReaction->CreateEvent(index);
540 fJumpingNumber++; 535 fJumpingNumber++;
541 } 536 }
542 else 537 else
543 { 538 {
544 G4ExceptionDescription exceptionDescriptio 539 G4ExceptionDescription exceptionDescription;
545 exceptionDescription << "pJumping == nullp 540 exceptionDescription << "pJumping == nullptr && pReaction == nullptr";
546 G4Exception("G4DNAEventScheduler::Stepping 541 G4Exception("G4DNAEventScheduler::Stepping", "G4DNAEventScheduler003",
547 FatalErrorInArgument, exceptio 542 FatalErrorInArgument, exceptionDescription);
548 } 543 }
549 544
550 if(fVerbose > 1) 545 if(fVerbose > 1)
551 { 546 {
552 G4cout << "G4DNAEventScheduler::Stepping:: 547 G4cout << "G4DNAEventScheduler::Stepping::end "
553 "Print************************** 548 "Print***********************************"
554 << G4endl; 549 << G4endl;
555 G4cout << G4endl; 550 G4cout << G4endl;
556 } 551 }
557 fStepNumberInMesh++; 552 fStepNumberInMesh++;
558 } 553 }
559 554
560 void G4DNAEventScheduler::SetEndTime(const G4d 555 void G4DNAEventScheduler::SetEndTime(const G4double& endTime)
561 { 556 {
562 fEndTime = endTime; 557 fEndTime = endTime;
563 } 558 }
564 559
565 void G4DNAEventScheduler::RecordTime() 560 void G4DNAEventScheduler::RecordTime()
566 { 561 {
567 auto recordTime = *fLastRecoredTime; 562 auto recordTime = *fLastRecoredTime;
568 if(fGlobalTime >= recordTime && fCounterMap[ 563 if(fGlobalTime >= recordTime && fCounterMap[recordTime].empty())
569 { 564 {
570 auto begin = fpMesh->begin(); 565 auto begin = fpMesh->begin();
571 auto end = fpMesh->end(); 566 auto end = fpMesh->end();
572 for(; begin != end; begin++) 567 for(; begin != end; begin++)
573 { 568 {
574 const auto& mapData = std::get<2>(*begin 569 const auto& mapData = std::get<2>(*begin);
575 if(mapData.empty()) 570 if(mapData.empty())
576 { 571 {
577 continue; 572 continue;
578 } 573 }
579 for(const auto& it : mapData) 574 for(const auto& it : mapData)
580 { 575 {
581 fCounterMap[recordTime][it.first] += i 576 fCounterMap[recordTime][it.first] += it.second;
582 } 577 }
583 } 578 }
584 fLastRecoredTime++; 579 fLastRecoredTime++;
585 } 580 }
586 } 581 }
587 582
588 void G4DNAEventScheduler::PrintRecordTime() 583 void G4DNAEventScheduler::PrintRecordTime()
589 { 584 {
590 G4cout << "CounterMap.size : " << fCounterMa 585 G4cout << "CounterMap.size : " << fCounterMap.size() << G4endl;
591 for(const auto& i : fCounterMap) 586 for(const auto& i : fCounterMap)
592 { 587 {
593 auto map = i.second; 588 auto map = i.second;
594 auto begin = map.begin(); // 589 auto begin = map.begin(); //
595 auto end = map.end(); // 590 auto end = map.end(); //
596 for(; begin != end; begin++) 591 for(; begin != end; begin++)
597 { 592 {
598 auto molecule = begin->first; 593 auto molecule = begin->first;
599 auto number = begin->second; 594 auto number = begin->second;
600 if(number == 0) 595 if(number == 0)
601 { 596 {
602 continue; 597 continue;
603 } 598 }
604 G4cout << "molecule : " << molecule->Get 599 G4cout << "molecule : " << molecule->GetName() << " number : " << number
605 << G4endl; 600 << G4endl;
606 } 601 }
607 } 602 }
608 } 603 }
609 604
610 std::map<G4double /*time*/, G4DNAEventSchedule 605 std::map<G4double /*time*/, G4DNAEventScheduler::MapCounter>
611 G4DNAEventScheduler::GetCounterMap() const 606 G4DNAEventScheduler::GetCounterMap() const
612 { 607 {
613 return fCounterMap; 608 return fCounterMap;
614 } 609 }
615 610
616 void G4DNAEventScheduler::SetUserMeshAction( 611 void G4DNAEventScheduler::SetUserMeshAction(
617 std::unique_ptr<G4UserMeshAction> pUserMeshA 612 std::unique_ptr<G4UserMeshAction> pUserMeshAction)
618 { 613 {
619 fpUserMeshAction = std::move(pUserMeshAction 614 fpUserMeshAction = std::move(pUserMeshAction);
620 } 615 }
621 616
622 G4DNAMesh* G4DNAEventScheduler::GetMesh() cons 617 G4DNAMesh* G4DNAEventScheduler::GetMesh() const { return fpMesh.get(); }
623 618
624 G4int G4DNAEventScheduler::GetPixels() const { 619 G4int G4DNAEventScheduler::GetPixels() const { return fPixel; }
625 620
626 G4bool G4DNAEventScheduler::CheckingReactionRa 621 G4bool G4DNAEventScheduler::CheckingReactionRadius(G4double resolution)
627 { 622 {
628 auto pMolecularReactionTable = G4DNAMolecula 623 auto pMolecularReactionTable = G4DNAMolecularReactionTable::Instance();
629 auto reactionDataList = pMolecularReactionTa 624 auto reactionDataList = pMolecularReactionTable->GetVectorOfReactionData();
630 if(reactionDataList.empty()) 625 if(reactionDataList.empty())
631 { 626 {
632 G4cout << "reactionDataList.empty()" << G4 627 G4cout << "reactionDataList.empty()" << G4endl;
633 return true; 628 return true;
634 } 629 }
635 << 630 else
636 for(auto it : reactionDataList) <<
637 { <<
638 if(it->GetEffectiveReactionRadius() >= res <<
639 { <<
640 G4cout << it->GetReactant1()->GetName() <<
641 << it->GetReactant2()->GetName() <<
642 G4cout << "G4DNAEventScheduler::Reaction <<
643 << G4BestUnit(it->GetEffectiveRea <<
644 << G4endl; <<
645 G4cout << "resolution : " << G4BestUnit( <<
646 return false; <<
647 } <<
648 } <<
649 return true; <<
650 } <<
651 <<
652 void G4DNAEventScheduler::ResetEventSet() <<
653 { <<
654 fpEventSet->RemoveEventSet(); <<
655 fpGillespieReaction->CreateEvents(); <<
656 } <<
657 <<
658 void G4DNAEventScheduler::LastRegisterForCount <<
659 { <<
660 if(fLastRecoredTime == fTimeToRecord.end()) <<
661 { <<
662 //fully recorded -> happy ending <<
663 return; <<
664 }else <<
665 { 631 {
666 auto lastRecorded = --fLastRecoredTime;//f << 632 for(auto it : reactionDataList)
667 while (fLastRecoredTime != fTimeToRecord.e <<
668 { 633 {
669 fCounterMap[*fLastRecoredTime] = fCounte << 634 if(it->GetEffectiveReactionRadius() >= resolution / CLHEP::pi)
670 fLastRecoredTime++; << 635 {
>> 636 G4cout << it->GetReactant1()->GetName() << " + "
>> 637 << it->GetReactant2()->GetName() << G4endl;
>> 638 G4cout << "G4DNAEventScheduler::ReactionRadius : "
>> 639 << G4BestUnit(it->GetEffectiveReactionRadius(), "Length")
>> 640 << G4endl;
>> 641 G4cout << "resolution : " << G4BestUnit(resolution, "Length") << G4endl;
>> 642 return false;
>> 643 }
671 } 644 }
>> 645 return true;
672 } 646 }
673 <<
674 } 647 }