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