| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
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 G4DNAGillespieDirectMethod())
43 , fpEventSet(new G4DNAEventSet())
44 , fpUpdateSystem(new G4DNAUpdateSystemModel())
45 {}
46
47 void G4DNAEventScheduler::ClearAndReChargeCounter()
48 {
49 fCounterMap.clear();
50 if(fTimeToRecord.empty())
51 {
52 G4String WarMessage = "fTimeToRecord is empty ";
53 G4Exception("G4DNAEventScheduler::ClearAndReChargeCounter()",
54 "TimeToRecord is empty", JustWarning, WarMessage);
55 }
56 fLastRecoredTime = fTimeToRecord.begin();
57
58 if(G4VMoleculeCounter::Instance()->InUse()) // copy from MoleculeCounter
59 {
60 G4MoleculeCounter::RecordedMolecules species;
61 species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
62 if(species.get() == nullptr)
63 {
64 return;
65 }
66 if(species->empty())
67 {
68 G4MoleculeCounter::Instance()->ResetCounter();
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::Instance()->GetNMoleculesAtTime(
81 molecule, time_mol);
82
83 if(n_mol < 0)
84 {
85 G4cerr << "G4DNAEventScheduler::ClearAndReChargeCounter() ::N "
86 "molecules not valid < 0 "
87 << G4endl;
88 G4Exception("", "N<0", FatalException, "");
89 }
90 fCounterMap[time_mol][molecule] = n_mol;
91 }
92 fLastRecoredTime++;
93 }
94 G4MoleculeCounter::Instance()->ResetCounter(); // reset
95 G4MoleculeCounter::Instance()->Use(false); // no more used
96 }
97 else
98 {
99 G4ExceptionDescription exceptionDescription;
100 exceptionDescription << "G4VMoleculeCounter is not used";
101 G4Exception("G4DNAEventScheduler::ClearAndReChargeCounter()",
102 "G4DNAEventScheduler010", JustWarning, exceptionDescription);
103 }
104 }
105
106 [[maybe_unused]] void G4DNAEventScheduler::AddTimeToRecord(const G4double& time)
107 {
108 if(fTimeToRecord.find(time) == fTimeToRecord.end())
109 {
110 fTimeToRecord.insert(time);
111 }
112 }
113
114 G4DNAEventScheduler::~G4DNAEventScheduler() = default;
115
116 void G4DNAEventScheduler::Voxelizing()
117 {
118 auto pMainList = G4ITTrackHolder::Instance()->GetMainList();
119 std::map<G4VDNAMesh::Index, MapList> TrackKeyMap;
120 for(auto track : *pMainList)
121 {
122 auto molType = GetMolecule(track)->GetMolecularConfiguration();
123
124 auto pScavengerMaterial = dynamic_cast<G4DNAScavengerMaterial*>(
125 G4Scheduler::Instance()->GetScavengerMaterial());
126 if(pScavengerMaterial != nullptr &&
127 pScavengerMaterial->find(molType)) // avoid voxelize the scavenger
128 {
129 continue;
130 }
131
132 auto key = fpMesh->GetIndex(track->GetPosition());
133 if(TrackKeyMap.find(key) != TrackKeyMap.end())
134 {
135 std::map<MolType, size_t>& TrackTypeMap = TrackKeyMap[key];
136 if(TrackTypeMap.find(molType) != TrackTypeMap.end())
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::move(it.second));
154 }
155 }
156
157 void G4DNAEventScheduler::ReVoxelizing(G4int pixel)
158 {
159 fPixel = pixel;
160 auto newMesh = new G4DNAMesh(fpMesh->GetBoundingBox(), fPixel);
161
162 auto begin = fpMesh->begin();
163 auto end = fpMesh->end();
164 std::map<G4VDNAMesh::Index, MapList> TrackKeyMap;
165 for(; begin != end; begin++)
166 {
167 auto index = std::get<0>(*begin);
168 auto newIndex = fpMesh->ConvertIndex(index, fPixel);
169 if(TrackKeyMap.find(newIndex) == TrackKeyMap.end())
170 {
171 TrackKeyMap[newIndex] = std::get<2>(*begin);
172 }
173 else
174 {
175 for(const auto& it : std::get<2>(*begin))
176 {
177 TrackKeyMap[newIndex][it.first] += it.second;
178 }
179 if(fVerbose > 1)
180 {
181 G4cout << " ReVoxelizing:: Old index : " << index
182 << " new index : " << fpMesh->ConvertIndex(index, fPixel)
183 << " number: " << std::get<2>(*begin).begin()->second << G4endl;
184 }
185 }
186 }
187 fpMesh.reset(newMesh);
188
189 for(auto& it : TrackKeyMap)
190 {
191 fpMesh->InitializeVoxel(it.first, std::move(it.second));
192 }
193 }
194 void G4DNAEventScheduler::Reset()
195 {
196 // find another solultion
197 fGlobalTime = fEndTime;
198
199 //
200 LastRegisterForCounter();//Last register for counter
201
202 if(fVerbose > 0)
203 {
204 G4cout << "End Processing and reset Gird, ScavengerTable, EventSet for new "
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 G4DNABoundingBox& boundingBox,
222 G4int pixel)
223 {
224 if(!fInitialized)
225 {
226 fPixel = pixel;
227 fpMesh = std::make_unique<G4DNAMesh>(boundingBox, pixel);
228
229 if(!CheckingReactionRadius(fpMesh->GetResolution()))
230 {
231 G4String WarMessage = "resolution is not good : " +
232 std::to_string(fpMesh->GetResolution() / nm);
233 G4Exception("G4DNAEventScheduler::InitializeInMesh()", "WrongResolution",
234 JustWarning, WarMessage);
235 }
236
237 // Scavenger();
238
239 auto pScavengerMaterial = dynamic_cast<G4DNAScavengerMaterial*>(
240 G4Scheduler::Instance()->GetScavengerMaterial());
241 if(pScavengerMaterial == nullptr)
242 {
243 G4cout << "There is no scavenger" << G4endl;
244 }
245 else
246 {
247 if(fVerbose > 1)
248 {
249 pScavengerMaterial->PrintInfo();
250 }
251 }
252
253 Voxelizing();
254 fpGillespieReaction->SetVoxelMesh(*fpMesh);
255 fpGillespieReaction->SetEventSet(fpEventSet.get());
256 fpGillespieReaction->SetTimeStep(0);// reset fTimeStep = 0 in fpGillespieReaction
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 Mesh, EventSet for new Mesh!!!!"
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() const { return fStartTime; }
307
308 G4double G4DNAEventScheduler::GetEndTime() const { return fEndTime; }
309
310 [[maybe_unused]] G4double G4DNAEventScheduler::GetTimeStep() const
311 {
312 return fTimeStep;
313 }
314
315 G4int G4DNAEventScheduler::GetVerbose() const { return fVerbose; }
316
317 [[maybe_unused]] void G4DNAEventScheduler::SetMaxNbSteps(G4int max)
318 {
319 fMaxStep = max;
320 }
321
322 [[maybe_unused]] void G4DNAEventScheduler::SetStartTime(G4double time)
323 {
324 fStartTime = time;
325 fGlobalTime = fStartTime;
326 }
327
328 void G4DNAEventScheduler::Stop() { fRunning = false; }
329 void G4DNAEventScheduler::Run()
330 {
331 G4Timer localtimer;
332 if(fVerbose > 2)
333 {
334 localtimer.Start();
335 G4cout << "***G4DNAEventScheduler::Run*** for Pixel : " << fPixel << G4endl;
336 }
337 while(fEndTime > fGlobalTime && fRunning)
338 {
339 RunInMesh();
340 }
341 if(fVerbose > 2)
342 {
343 if(!fRunning)
344 {
345 G4cout << " StepNumber(" << fStepNumber << ") = MaxStep(" << fMaxStep
346 << ")" << G4endl;
347 }
348 else if(fEndTime <= fGlobalTime)
349 {
350 G4cout << " GlobalTime(" << fGlobalTime << ") > EndTime(" << fEndTime
351 << ")"
352 << " StepNumber : " << fStepNumber << G4endl;
353 }
354 localtimer.Stop();
355 G4cout << "***G4DNAEventScheduler::Ending::"
356 << G4BestUnit(fGlobalTime, "Time")
357 << " Events left : " << fpEventSet->size() << G4endl;
358 if(fVerbose > 1)
359 {
360 fpMesh->PrintMesh();
361 }
362 G4cout << " Computing Time : " << localtimer << G4endl;
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->GetResolution();
376 G4cout << "At Time : " << std::setw(7) << G4BestUnit(fGlobalTime, "Time")
377 << " the Mesh has " << fPixel << " x " << fPixel << " x " << fPixel
378 << " voxels with Resolution " << G4BestUnit(resolution, "Length")
379 << " during next "
380 << G4BestUnit(resolution * resolution * C / (6 * D), "Time")
381 << G4endl;
382 }
383
384 if(fVerbose > 2)
385 {
386 fpMesh->PrintMesh();
387 }
388
389 if(fpUserMeshAction != nullptr)
390 {
391 fpUserMeshAction->BeginOfMesh(fpMesh.get(), fGlobalTime);
392 }
393
394 // if diffusive jumping is avaiable, EventSet is never empty
395
396 while(!fpEventSet->Empty() && !fIsChangeMesh && fEndTime > fGlobalTime)
397 {
398 Stepping();
399 fGlobalTime = fTimeStep + fStartTime;
400
401 if(fpUserMeshAction != nullptr)
402 {
403 fpUserMeshAction->InMesh(fpMesh.get(), fGlobalTime);
404 }
405
406 if(fVerbose > 2)
407 {
408 G4cout << "fGlobalTime : " << G4BestUnit(fGlobalTime, "Time")
409 << " fTimeStep : " << G4BestUnit(fTimeStep, "Time") << G4endl;
410 }
411 G4double resolution = fpMesh->GetResolution();
412 fTransferTime = resolution * resolution * C / (6 * D);
413 if(fTransferTime == 0)
414 {
415 G4ExceptionDescription exceptionDescription;
416 exceptionDescription << "fTransferTime == 0";
417 G4Exception("G4DNAEventScheduler::RunInMesh", "G4DNAEventScheduler001",
418 FatalErrorInArgument, exceptionDescription);
419 }
420 if(fTransferTime < fTimeStep &&
421 fPixel != 1) // dont change Mesh if fPixel == 1
422 {
423 if(fSetChangeMesh)
424 {
425 if(fVerbose > 1)
426 {
427 G4cout << " Pixels : " << fPixel << " resolution : "
428 << G4BestUnit(fpMesh->GetResolution(), "Length")
429 << " fStepNumberInMesh : " << fStepNumberInMesh
430 << " at fGlobalTime : " << G4BestUnit(fGlobalTime, "Time")
431 << " at fTimeStep : " << G4BestUnit(fTimeStep, "Time")
432 << " fReactionNumber : " << fReactionNumber
433 << " transferTime : " << G4BestUnit(fTransferTime, "Time")
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->size() << G4endl;
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 : " << fPixel
454 << " pixels to : " << fPixel / 2 << " pixels" << G4endl;
455 G4cout << "Info : ReactionNumber : " << fReactionNumber
456 << " JumpingNumber : " << fJumpingNumber << G4endl;
457 }
458 else if(fEndTime > fGlobalTime)
459 {
460 G4cout << " GlobalTime(" << fGlobalTime << ") > EndTime(" << fEndTime
461 << ")"
462 << " StepNumber : " << fStepNumber << G4endl;
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(), fGlobalTime);
474 }
475 ResetInMesh();
476 }
477
478 void G4DNAEventScheduler::Stepping() // this event loop
479 {
480 fStepNumber < fMaxStep ? fStepNumber++ : static_cast<int>(fRunning = false);
481 if(fpEventSet->size() > fpMesh->size())
482 {
483 G4ExceptionDescription exceptionDescription;
484 exceptionDescription
485 << "impossible that fpEventSet->size() > fpMesh->size()";
486 G4Exception("G4DNAEventScheduler::Stepping", "G4DNAEventScheduler002",
487 FatalErrorInArgument, exceptionDescription);
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)->GetReactionData();
507
508 fpUpdateSystem->SetGlobalTime(fTimeStep +
509 fStartTime); // this is just for printing
510 fpGillespieReaction->SetTimeStep(fTimeStep);
511 if(pJumping == nullptr && pReaction != nullptr)
512 {
513 fpUpdateSystem->UpdateSystem(index, *pReaction);
514 fpEventSet->RemoveEvent(selected);
515
516 //hoang's exp
517 if(fpGillespieReaction->SetEquilibrium(pReaction))
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 == nullptr)
530 {
531 // dont change this
532 fpUpdateSystem->UpdateSystem(index, *pJumping);
533 // save jumping Index before delete selected event
534 auto jumpingIndex = pJumping->second;
535 fpEventSet->RemoveEvent(selected);
536 // create new event
537 // should create Jumping before key
538 fpGillespieReaction->CreateEvent(jumpingIndex);
539 fpGillespieReaction->CreateEvent(index);
540 fJumpingNumber++;
541 }
542 else
543 {
544 G4ExceptionDescription exceptionDescription;
545 exceptionDescription << "pJumping == nullptr && pReaction == nullptr";
546 G4Exception("G4DNAEventScheduler::Stepping", "G4DNAEventScheduler003",
547 FatalErrorInArgument, exceptionDescription);
548 }
549
550 if(fVerbose > 1)
551 {
552 G4cout << "G4DNAEventScheduler::Stepping::end "
553 "Print***********************************"
554 << G4endl;
555 G4cout << G4endl;
556 }
557 fStepNumberInMesh++;
558 }
559
560 void G4DNAEventScheduler::SetEndTime(const G4double& endTime)
561 {
562 fEndTime = endTime;
563 }
564
565 void G4DNAEventScheduler::RecordTime()
566 {
567 auto recordTime = *fLastRecoredTime;
568 if(fGlobalTime >= recordTime && fCounterMap[recordTime].empty())
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] += it.second;
582 }
583 }
584 fLastRecoredTime++;
585 }
586 }
587
588 void G4DNAEventScheduler::PrintRecordTime()
589 {
590 G4cout << "CounterMap.size : " << fCounterMap.size() << G4endl;
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->GetName() << " number : " << number
605 << G4endl;
606 }
607 }
608 }
609
610 std::map<G4double /*time*/, G4DNAEventScheduler::MapCounter>
611 G4DNAEventScheduler::GetCounterMap() const
612 {
613 return fCounterMap;
614 }
615
616 void G4DNAEventScheduler::SetUserMeshAction(
617 std::unique_ptr<G4UserMeshAction> pUserMeshAction)
618 {
619 fpUserMeshAction = std::move(pUserMeshAction);
620 }
621
622 G4DNAMesh* G4DNAEventScheduler::GetMesh() const { return fpMesh.get(); }
623
624 G4int G4DNAEventScheduler::GetPixels() const { return fPixel; }
625
626 G4bool G4DNAEventScheduler::CheckingReactionRadius(G4double resolution)
627 {
628 auto pMolecularReactionTable = G4DNAMolecularReactionTable::Instance();
629 auto reactionDataList = pMolecularReactionTable->GetVectorOfReactionData();
630 if(reactionDataList.empty())
631 {
632 G4cout << "reactionDataList.empty()" << G4endl;
633 return true;
634 }
635
636 for(auto it : reactionDataList)
637 {
638 if(it->GetEffectiveReactionRadius() >= resolution / CLHEP::pi)
639 {
640 G4cout << it->GetReactant1()->GetName() << " + "
641 << it->GetReactant2()->GetName() << G4endl;
642 G4cout << "G4DNAEventScheduler::ReactionRadius : "
643 << G4BestUnit(it->GetEffectiveReactionRadius(), "Length")
644 << G4endl;
645 G4cout << "resolution : " << G4BestUnit(resolution, "Length") << G4endl;
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::LastRegisterForCounter()
659 {
660 if(fLastRecoredTime == fTimeToRecord.end())
661 {
662 //fully recorded -> happy ending
663 return;
664 }else
665 {
666 auto lastRecorded = --fLastRecoredTime;//fixed
667 while (fLastRecoredTime != fTimeToRecord.end())
668 {
669 fCounterMap[*fLastRecoredTime] = fCounterMap[*lastRecorded];
670 fLastRecoredTime++;
671 }
672 }
673
674 }