| 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
28 #include "G4MoleculeCounter.hh"
29
30 #include "G4MolecularConfiguration.hh"
31 #include "G4MoleculeDefinition.hh"
32 #include "G4MoleculeTable.hh"
33 #include "G4Scheduler.hh" // TODO: remove this
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
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
57 }
58 }
59
60 //--------------------------------------------
61 G4MoleculeCounter* G4MoleculeCounter::Instance
62 {
63 if (fpInstance == nullptr)
64 {
65 fpInstance = new G4MoleculeCounter();
66 }
67 return dynamic_cast<G4MoleculeCounter*>(fp
68 }
69
70 //--------------------------------------------
71
72 G4MoleculeCounter::G4MoleculeCounter()
73 {
74 fVerbose = 0;
75 fCheckTimeIsConsistentWithScheduler = true
76 }
77
78 //--------------------------------------------
79
80 G4MoleculeCounter::~G4MoleculeCounter() = defa
81
82 //--------------------------------------------
83
84 void G4MoleculeCounter::Initialize()
85 {
86 auto mol_iterator = G4MoleculeTable::Insta
87 while ((mol_iterator)())
88 {
89 if (!IsRegistered(mol_iterator.value()
90 {
91 continue;
92 }
93
94 fCounterMap[mol_iterator.value()]; //
95 }
96 }
97
98 //--------------------------------------------
99
100 void G4MoleculeCounter::SetTimeSlice(double ti
101 {
102 G4::MoleculeCounter::TimePrecision::fPreci
103 }
104
105 //--------------------------------------------
106
107 G4bool G4MoleculeCounter::SearchTimeMap(Reacta
108 {
109 if (fpLastSearch == nullptr)
110 {
111 fpLastSearch = std::make_unique<Search
112 }
113 else
114 {
115 if (fpLastSearch->fLowerBoundSet &&
116 fpLastSearch->fLastMoleculeSearche
117 {
118 return true;
119 }
120 }
121
122 auto mol_it = fCounterMap.find(molecule);
123 fpLastSearch->fLastMoleculeSearched = mol_
124
125 if (mol_it != fCounterMap.end())
126 {
127 fpLastSearch->fLowerBoundTime = fpLast
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(do
142 bo
143 {
144 auto mol_it = fpLastSearch->fLastMoleculeS
145 if (mol_it == fCounterMap.end())
146 {
147 return 0;
148 }
149
150 NbMoleculeAgainstTime& timeMap = mol_it->s
151 if (timeMap.empty())
152 {
153 return 0;
154 }
155
156 if (sameTypeOfMolecule)
157 {
158 if (fpLastSearch->fLowerBoundSet && fp
159 {
160 if (fpLastSearch->fLowerBoundTime-
161 {
162 auto upperToLast = fpLastSearc
163 upperToLast++;
164
165 if (upperToLast == timeMap.end
166 {
167 return fpLastSearch->fLowe
168 }
169
170 if (upperToLast->first > time)
171 {
172 return fpLastSearch->fLowe
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->seco
196 }
197
198 //--------------------------------------------
199
200 int G4MoleculeCounter::GetNMoleculesAtTime(Rea
201 dou
202 {
203 G4bool sameTypeOfMolecule = SearchTimeMap(
204 return SearchUpperBoundTime(time, sameType
205 }
206
207 //--------------------------------------------
208
209 void G4MoleculeCounter::AddAMoleculeAtTime(Rea
210 G4d
211 con
212 int
213 {
214 if (fDontRegister[molecule->GetDefinition(
215 {
216 return;
217 }
218
219 if (fVerbose != 0)
220 {
221 G4cout << "G4MoleculeCounter::AddAMole
222 << " at time : " << G4BestUnit(
223 }
224
225 auto counterMap_i = fCounterMap.find(molec
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::Mol
241 // Case 1 = new time comes after l
242 // Case 2 = new time is about the
243 {
244 double newValue = end->second + nu
245 counterMap_i->second[time] = newVa
246 }
247 else
248 {
249 // if(fabs(time - G4Scheduler
250 // G4Scheduler::Instance()
251 {
252 G4ExceptionDescription errMsg;
253 errMsg << "Time of species "
254 << molecule->GetName()
255 << G4BestUnit(time, "Ti
256 << " global time is "
257 << G4BestUnit(G4Schedul
258 << G4endl;
259 G4Exception("G4MoleculeCounter
260 "TIME_DONT_MATCH",
261 FatalException, er
262 }
263 }
264 }
265 }
266
267 //--------------------------------------------
268
269 void G4MoleculeCounter::RemoveAMoleculeAtTime(
270
271
272
273 {
274 if (fDontRegister[pMolecule->GetDefinition
275 {
276 return;
277 }
278
279 if (fVerbose != 0)
280 {
281 G4cout << "G4MoleculeCounter::RemoveAM
282 << pMolecule->GetName() << " at
283 << G4endl;
284 }
285
286 if (fCheckTimeIsConsistentWithScheduler)
287 {
288 if (fabs(time - G4Scheduler::Instance(
289 G4Scheduler::Instance()->GetTimeTo
290 {
291 G4ExceptionDescription errMsg;
292 errMsg << "Time of species "
293 << pMolecule->GetName() <<
294 << G4BestUnit(time, "Time")
295 << " global time is "
296 << G4BestUnit(G4Scheduler::
297 << G4endl;
298 G4Exception("G4MoleculeCounter::Re
299 "TIME_DONT_MATCH",
300 FatalException, errMsg
301 }
302 }
303
304 NbMoleculeAgainstTime& nbMolPerTime = fCou
305
306 if (nbMolPerTime.empty())
307 {
308 pMolecule->PrintState();
309 Dump();
310 G4String errMsg =
311 "You are trying to remove mole
312 " from the counter while this
313 G4Exception("G4MoleculeCounter::Remove
314 FatalErrorInArgument, errM
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->G
327 G4Exception("G4MoleculeCounter::Remove
328 FatalErrorInArgument, errM
329 }
330
331 if (time - it->first < -G4::MoleculeCounte
332 {
333 Dump();
334 G4ExceptionDescription errMsg;
335 errMsg << "Is time going back?? " << p
336 << " is being removed at time "
337 << " while last recorded time w
338 << G4BestUnit(it->first, "Time"
339 G4Exception("G4MoleculeCounter::Remove
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 " << numbe
352 << pMolecule->GetName() << " th
353 << G4BestUnit(time, "Time") <<
354 << " Global time is "
355 << G4BestUnit(G4Scheduler::Inst
356 << ". Previous selected time is
357 << G4BestUnit(it->first, "Time"
358 << G4endl;
359 G4Exception("G4MoleculeCounter::Remove
360 "N_INF_0",
361 FatalException, errMsg);
362 }
363
364 nbMolPerTime[time] = finalN;
365 }
366
367 //--------------------------------------------
368
369 G4MoleculeCounter::RecordedMolecules G4Molecul
370 {
371 if (fVerbose > 1)
372 {
373 G4cout << "Entering in G4MoleculeCount
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::GetRecordedTi
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(G4Spe
407 // size_
408 // int /
409 // G4Spe
410 // int d
411 //{
412 // switch(speciesChange)
413 // {
414 // case G4SpeciesInCM::eAdd:
415 // AddAMoleculeAtTime(G4MoleculeTable::In
416 // G4Scheduler::Instan
417 // diff);
418 // break;
419 // case G4SpeciesInCM::eRemove:
420 // RemoveAMoleculeAtTime(G4MoleculeTable:
421 // G4Scheduler::Instan
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->
436
437 for (const auto& it2 : it.second)
438 {
439 G4cout << " " << G4BestUnit(it2.fi
440 << " " << it2.second <<
441 }
442 }
443 }
444
445 //--------------------------------------------
446
447 void G4MoleculeCounter::ResetCounter()
448 {
449 if (fVerbose != 0)
450 {
451 G4cout << " ---> G4MoleculeCounter::Re
452 }
453 fCounterMap.clear();
454 fpLastSearch.reset(nullptr);
455 }
456
457 //--------------------------------------------
458
459 const NbMoleculeAgainstTime& G4MoleculeCounter
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 G4M
481 {
482 fDontRegister[molDef] = true;
483 }
484
485 //--------------------------------------------
486
487 bool G4MoleculeCounter::IsRegistered(const G4M
488 {
489 return fDontRegister.find(molDef) == fDont
490 }
491
492 //--------------------------------------------
493
494 void G4MoleculeCounter::RegisterAll()
495 {
496 fDontRegister.clear();
497 }
498
499 G4bool G4MoleculeCounter::IsTimeCheckedForCons
500 {
501 return fCheckTimeIsConsistentWithScheduler
502 }
503
504 void G4MoleculeCounter::CheckTimeForConsistenc
505 {
506 fCheckTimeIsConsistentWithScheduler = flag
507 }
508
509