| 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
28 #include "G4MoleculeCounter.hh"
29
30 #include "G4MolecularConfiguration.hh"
31 #include "G4MoleculeDefinition.hh"
32 #include "G4MoleculeTable.hh"
33 #include "G4Scheduler.hh" // TODO: remove this dependency
34 #include "G4SystemOfUnits.hh"
35 #include "G4UIcommand.hh"
36 #include "G4UnitsTable.hh"
37
38 #include <iomanip>
39 #include <memory>
40
41 using namespace std;
42
43 namespace G4{
44 namespace MoleculeCounter {
45
46 bool TimePrecision::operator()(const double& a, const double& b) const
47 {
48 if (std::fabs(a - b) < fPrecision)
49 {
50 return false;
51 }
52
53 return a < b;
54 }
55
56 G4ThreadLocal double TimePrecision::fPrecision = 0.5 * picosecond;
57 }
58 }
59
60 //------------------------------------------------------------------------------
61 G4MoleculeCounter* G4MoleculeCounter::Instance()
62 {
63 if (fpInstance == nullptr)
64 {
65 fpInstance = new G4MoleculeCounter();
66 }
67 return dynamic_cast<G4MoleculeCounter*>(fpInstance);
68 }
69
70 //------------------------------------------------------------------------------
71
72 G4MoleculeCounter::G4MoleculeCounter()
73 {
74 fVerbose = 0;
75 fCheckTimeIsConsistentWithScheduler = true;
76 }
77
78 //------------------------------------------------------------------------------
79
80 G4MoleculeCounter::~G4MoleculeCounter() = default;
81
82 //------------------------------------------------------------------------------
83
84 void G4MoleculeCounter::Initialize()
85 {
86 auto mol_iterator = G4MoleculeTable::Instance()->GetConfigurationIterator();
87 while ((mol_iterator)())
88 {
89 if (!IsRegistered(mol_iterator.value()->GetDefinition()))
90 {
91 continue;
92 }
93
94 fCounterMap[mol_iterator.value()]; // initialize the second map
95 }
96 }
97
98 //------------------------------------------------------------------------------
99
100 void G4MoleculeCounter::SetTimeSlice(double timeSlice)
101 {
102 G4::MoleculeCounter::TimePrecision::fPrecision = timeSlice;
103 }
104
105 //------------------------------------------------------------------------------
106
107 G4bool G4MoleculeCounter::SearchTimeMap(Reactant* molecule)
108 {
109 if (fpLastSearch == nullptr)
110 {
111 fpLastSearch = std::make_unique<Search>();
112 }
113 else
114 {
115 if (fpLastSearch->fLowerBoundSet &&
116 fpLastSearch->fLastMoleculeSearched->first == molecule)
117 {
118 return true;
119 }
120 }
121
122 auto mol_it = fCounterMap.find(molecule);
123 fpLastSearch->fLastMoleculeSearched = mol_it;
124
125 if (mol_it != fCounterMap.end())
126 {
127 fpLastSearch->fLowerBoundTime = fpLastSearch->fLastMoleculeSearched->second
128 .end();
129 fpLastSearch->fLowerBoundSet = true;
130 }
131 else
132 {
133 fpLastSearch->fLowerBoundSet = false;
134 }
135
136 return false;
137 }
138
139 //------------------------------------------------------------------------------
140
141 int G4MoleculeCounter::SearchUpperBoundTime(double time,
142 bool sameTypeOfMolecule)
143 {
144 auto mol_it = fpLastSearch->fLastMoleculeSearched;
145 if (mol_it == fCounterMap.end())
146 {
147 return 0;
148 }
149
150 NbMoleculeAgainstTime& timeMap = mol_it->second;
151 if (timeMap.empty())
152 {
153 return 0;
154 }
155
156 if (sameTypeOfMolecule)
157 {
158 if (fpLastSearch->fLowerBoundSet && fpLastSearch->fLowerBoundTime != timeMap.end())
159 {
160 if (fpLastSearch->fLowerBoundTime->first < time)
161 {
162 auto upperToLast = fpLastSearch->fLowerBoundTime;
163 upperToLast++;
164
165 if (upperToLast == timeMap.end())
166 {
167 return fpLastSearch->fLowerBoundTime->second;
168 }
169
170 if (upperToLast->first > time)
171 {
172 return fpLastSearch->fLowerBoundTime->second;
173 }
174 }
175 }
176 }
177
178 auto up_time_it = timeMap.upper_bound(time);
179
180 if (up_time_it == timeMap.end())
181 {
182 auto last_time = timeMap.rbegin();
183 return last_time->second;
184 }
185 if (up_time_it == timeMap.begin())
186 {
187 return 0;
188 }
189
190 up_time_it--;
191
192 fpLastSearch->fLowerBoundTime = up_time_it;
193 fpLastSearch->fLowerBoundSet = true;
194
195 return fpLastSearch->fLowerBoundTime->second;
196 }
197
198 //------------------------------------------------------------------------------
199
200 int G4MoleculeCounter::GetNMoleculesAtTime(Reactant* molecule,
201 double time)
202 {
203 G4bool sameTypeOfMolecule = SearchTimeMap(molecule);
204 return SearchUpperBoundTime(time, sameTypeOfMolecule);
205 }
206
207 //------------------------------------------------------------------------------
208
209 void G4MoleculeCounter::AddAMoleculeAtTime(Reactant* molecule,
210 G4double time,
211 const G4ThreeVector* /*position*/,
212 int number)
213 {
214 if (fDontRegister[molecule->GetDefinition()])
215 {
216 return;
217 }
218
219 if (fVerbose != 0)
220 {
221 G4cout << "G4MoleculeCounter::AddAMoleculeAtTime : " << molecule->GetName()
222 << " at time : " << G4BestUnit(time, "Time") << G4endl;
223 }
224
225 auto counterMap_i = fCounterMap.find(molecule);
226
227 if (counterMap_i == fCounterMap.end())
228 {
229 fCounterMap[molecule][time] = number;
230 }
231 else if (counterMap_i->second.empty())
232 {
233 counterMap_i->second[time] = number;
234 }
235 else
236 {
237 auto end = counterMap_i->second.rbegin();
238
239 if (end->first <= time ||
240 fabs(end->first - time) <= G4::MoleculeCounter::TimePrecision::fPrecision)
241 // Case 1 = new time comes after last recorded data
242 // Case 2 = new time is about the same as the last recorded one
243 {
244 double newValue = end->second + number;
245 counterMap_i->second[time] = newValue;
246 }
247 else
248 {
249 // if(fabs(time - G4Scheduler::Instance()->GetGlobalTime()) >
250 // G4Scheduler::Instance()->GetTimeTolerance())
251 {
252 G4ExceptionDescription errMsg;
253 errMsg << "Time of species "
254 << molecule->GetName() << " is "
255 << G4BestUnit(time, "Time") << " while "
256 << " global time is "
257 << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
258 << G4endl;
259 G4Exception("G4MoleculeCounter::AddAMoleculeAtTime",
260 "TIME_DONT_MATCH",
261 FatalException, errMsg);
262 }
263 }
264 }
265 }
266
267 //------------------------------------------------------------------------------
268
269 void G4MoleculeCounter::RemoveAMoleculeAtTime(const G4MolecularConfiguration* pMolecule,
270 G4double time,
271 const G4ThreeVector* /*position*/,
272 int number)
273 {
274 if (fDontRegister[pMolecule->GetDefinition()])
275 {
276 return;
277 }
278
279 if (fVerbose != 0)
280 {
281 G4cout << "G4MoleculeCounter::RemoveAMoleculeAtTime : "
282 << pMolecule->GetName() << " at time : " << G4BestUnit(time, "Time")
283 << G4endl;
284 }
285
286 if (fCheckTimeIsConsistentWithScheduler)
287 {
288 if (fabs(time - G4Scheduler::Instance()->GetGlobalTime()) >
289 G4Scheduler::Instance()->GetTimeTolerance())
290 {
291 G4ExceptionDescription errMsg;
292 errMsg << "Time of species "
293 << pMolecule->GetName() << " is "
294 << G4BestUnit(time, "Time") << " while "
295 << " global time is "
296 << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
297 << G4endl;
298 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
299 "TIME_DONT_MATCH",
300 FatalException, errMsg);
301 }
302 }
303
304 NbMoleculeAgainstTime& nbMolPerTime = fCounterMap[pMolecule];
305
306 if (nbMolPerTime.empty())
307 {
308 pMolecule->PrintState();
309 Dump();
310 G4String errMsg =
311 "You are trying to remove molecule " + pMolecule->GetName() +
312 " from the counter while this kind of molecules has not been registered yet";
313 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime", "",
314 FatalErrorInArgument, errMsg);
315
316 return;
317 }
318
319 auto it = nbMolPerTime.rbegin();
320
321 if (it == nbMolPerTime.rend())
322 {
323 it--;
324
325 G4String errMsg =
326 "There was no " + pMolecule->GetName() + " recorded at the time or even before the time asked";
327 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime", "",
328 FatalErrorInArgument, errMsg);
329 }
330
331 if (time - it->first < -G4::MoleculeCounter::TimePrecision::fPrecision)
332 {
333 Dump();
334 G4ExceptionDescription errMsg;
335 errMsg << "Is time going back?? " << pMolecule->GetName()
336 << " is being removed at time " << G4BestUnit(time, "Time")
337 << " while last recorded time was "
338 << G4BestUnit(it->first, "Time") << ".";
339 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
340 "RETURN_TO_THE_FUTUR",
341 FatalErrorInArgument,
342 errMsg);
343 }
344
345 double finalN = it->second - number;
346
347 if (finalN < 0)
348 {
349 Dump();
350 G4ExceptionDescription errMsg;
351 errMsg << "After removal of " << number << " species of "
352 << pMolecule->GetName() << " the final number at time "
353 << G4BestUnit(time, "Time") << " is less than zero and so not valid."
354 << " Global time is "
355 << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
356 << ". Previous selected time is "
357 << G4BestUnit(it->first, "Time")
358 << G4endl;
359 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
360 "N_INF_0",
361 FatalException, errMsg);
362 }
363
364 nbMolPerTime[time] = finalN;
365 }
366
367 //------------------------------------------------------------------------------
368
369 G4MoleculeCounter::RecordedMolecules G4MoleculeCounter::GetRecordedMolecules()
370 {
371 if (fVerbose > 1)
372 {
373 G4cout << "Entering in G4MoleculeCounter::RecordMolecules" << G4endl;
374 }
375
376 RecordedMolecules output(new ReactantList());
377
378 for (const auto & it : fCounterMap)
379 {
380 output->push_back(it.first);
381 }
382 return output;
383 }
384
385 //------------------------------------------------------------------------------
386
387 RecordedTimes G4MoleculeCounter::GetRecordedTimes()
388 {
389 RecordedTimes output(new std::set<G4double>);
390
391 for(const auto& it : fCounterMap)
392 {
393 for(const auto& it2 : it.second)
394 {
395 //time = it2->first;
396 output->insert(it2.first);
397 }
398 }
399
400 return output;
401 }
402
403 //------------------------------------------------------------------------------
404
405 // >>DEV<<
406 //void G4MoleculeCounter::SignalReceiver(G4SpeciesInCM* /*speciesInCM*/,
407 // size_t moleculeID,
408 // int /*number*/,
409 // G4SpeciesInCM::SpeciesChange speciesChange,
410 // int diff)
411 //{
412 // switch(speciesChange)
413 // {
414 // case G4SpeciesInCM::eAdd:
415 // AddAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
416 // G4Scheduler::Instance()->GetGlobalTime(),
417 // diff);
418 // break;
419 // case G4SpeciesInCM::eRemove:
420 // RemoveAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
421 // G4Scheduler::Instance()->GetGlobalTime(),
422 // diff);
423 // break;
424 // }
425 //}
426
427 //------------------------------------------------------------------------------
428
429 void G4MoleculeCounter::Dump()
430 {
431 for (const auto& it : fCounterMap)
432 {
433 auto pReactant = it.first;
434
435 G4cout << " --- > For " << pReactant->GetName() << G4endl;
436
437 for (const auto& it2 : it.second)
438 {
439 G4cout << " " << G4BestUnit(it2.first, "Time")
440 << " " << it2.second << G4endl;
441 }
442 }
443 }
444
445 //------------------------------------------------------------------------------
446
447 void G4MoleculeCounter::ResetCounter()
448 {
449 if (fVerbose != 0)
450 {
451 G4cout << " ---> G4MoleculeCounter::ResetCounter" << G4endl;
452 }
453 fCounterMap.clear();
454 fpLastSearch.reset(nullptr);
455 }
456
457 //------------------------------------------------------------------------------
458
459 const NbMoleculeAgainstTime& G4MoleculeCounter::GetNbMoleculeAgainstTime(Reactant* molecule)
460 {
461 return fCounterMap[molecule];
462 }
463
464 //------------------------------------------------------------------------------
465
466 void G4MoleculeCounter::SetVerbose(G4int level)
467 {
468 fVerbose = level;
469 }
470
471 //------------------------------------------------------------------------------
472
473 G4int G4MoleculeCounter::GetVerbose()
474 {
475 return fVerbose;
476 }
477
478 //------------------------------------------------------------------------------
479
480 void G4MoleculeCounter::DontRegister(const G4MoleculeDefinition* molDef)
481 {
482 fDontRegister[molDef] = true;
483 }
484
485 //------------------------------------------------------------------------------
486
487 bool G4MoleculeCounter::IsRegistered(const G4MoleculeDefinition* molDef)
488 {
489 return fDontRegister.find(molDef) == fDontRegister.end();
490 }
491
492 //------------------------------------------------------------------------------
493
494 void G4MoleculeCounter::RegisterAll()
495 {
496 fDontRegister.clear();
497 }
498
499 G4bool G4MoleculeCounter::IsTimeCheckedForConsistency() const
500 {
501 return fCheckTimeIsConsistentWithScheduler;
502 }
503
504 void G4MoleculeCounter::CheckTimeForConsistency(G4bool flag)
505 {
506 fCheckTimeIsConsistentWithScheduler = flag;
507 }
508
509