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