Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MoleculeCounter.cc
Go to the documentation of this file.
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 // $Id: G4MoleculeCounter.cc 101354 2016-11-15 08:27:51Z gcosmo $
27 //
28 
29 #include <iomanip>
30 #include "G4MoleculeCounter.hh"
31 #include "G4MoleculeTable.hh"
32 #include "G4UIcommand.hh"
33 #include "G4UnitsTable.hh"
35 #include "G4MoleculeDefinition.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4Scheduler.hh" // TODO: remove this dependency
38 
39 using namespace std;
40 
42 
43 //------------------------------------------------------------------------------
45  if (!fpInstance) fpInstance = new G4MoleculeCounter();
46  return dynamic_cast<G4MoleculeCounter*>(fpInstance);
47 }
48 
49 //------------------------------------------------------------------------------
50 
52 {
53  fVerbose = 0;
54  fCheckTimeIsConsistentWithScheduler = true;
56  {
58  }
59 }
60 
61 //------------------------------------------------------------------------------
62 
64 {
65 }
66 
67 //------------------------------------------------------------------------------
68 
70 {
71 // G4cout << "G4MoleculeCounter::Initialize" << G4endl;
72 
75  while ((mol_iterator)())
76  {
77  if(IsRegistered(mol_iterator.value()->GetDefinition()) == false)
78  {
79  continue;
80  }
81 
82  // G4cout << "G4MoleculeCounter::Initialize" << G4endl;
83  // G4cout << mol_iterator->value()->GetName() << G4endl;
84  fCounterMap[mol_iterator.value()]; // initialize the second map
85  }
86 }
87 
88 //------------------------------------------------------------------------------
89 
90 void G4MoleculeCounter::SetTimeSlice(double timeSlice)
91 {
93 }
94 
95 //------------------------------------------------------------------------------
96 
98 {
99  if (fpLastSearch.get() == 0)
100  {
101  fpLastSearch.reset(new Search());
102  }
103  else
104  {
105  if (fpLastSearch->fLowerBoundSet &&
106  fpLastSearch->fLastMoleculeSearched->first == molecule)
107  return true;
108  }
109 
110  CounterMapType::iterator mol_it = fCounterMap.find(molecule);
111  fpLastSearch->fLastMoleculeSearched = mol_it;
112 
113  if (mol_it != fCounterMap.end()) // TODO
114  {
115  fpLastSearch->fLowerBoundTime = fpLastSearch->fLastMoleculeSearched->second
116  .end();
117  fpLastSearch->fLowerBoundSet = true;
118  }
119  else
120  {
121  fpLastSearch->fLowerBoundSet = false;
122  }
123 
124  return false;
125 }
126 
127 //------------------------------------------------------------------------------
128 
130  bool sameTypeOfMolecule)
131 {
132  CounterMapType::iterator mol_it = fpLastSearch->fLastMoleculeSearched;
133  if (mol_it == fCounterMap.end()) return 0; // RETURN
134 
135  NbMoleculeAgainstTime& timeMap = mol_it->second;
136  if (timeMap.empty()) return 0;
137 
138  NbMoleculeAgainstTime::iterator end_time = timeMap.end();
139 
140  if (sameTypeOfMolecule == true)
141  {
142  //G4cout << "SAME MOLECULE" << G4endl;
143  if (fpLastSearch->fLowerBoundSet && fpLastSearch->fLowerBoundTime
144  != end_time)
145  {
146  if (fpLastSearch->fLowerBoundTime->first < time)
147  {
148  NbMoleculeAgainstTime::iterator upperToLast = fpLastSearch
149  ->fLowerBoundTime;
150  upperToLast++;
151 
152  if (upperToLast == end_time)
153  {
154  return fpLastSearch->fLowerBoundTime->second;
155  }
156 
157  if (upperToLast->first > time)
158  {
159  return fpLastSearch->fLowerBoundTime->second;
160  }
161  }
162  }
163  }
164  /*
165  else
166  {
167  G4cout << "--> Molecule has changed" << G4endl;
168  }
169  */
170  //G4cout << "Searching" << G4endl;
171  // With upper bound
172  NbMoleculeAgainstTime::iterator up_time_it = timeMap.upper_bound(time);
173 
174  if (up_time_it == end_time)
175  {
176  NbMoleculeAgainstTime::reverse_iterator last_time = timeMap.rbegin();
177 
178 // {
179  //G4cout << "RETURN LAST : " << G4BestUnit(time, "Time") << G4endl;
180  return last_time->second;
181 // }
182 
183 // G4cout << "RETURN 0 (1)" << G4endl;
184 // return 0; // RETURN
185  }
186  if (up_time_it == timeMap.begin())
187  {
188 // G4cout << "RETURN 0 (2)" << G4endl;
189  return 0; // RETURN
190  }
191 
192  //G4cout << "Going back : " << up_time_it->first << "-->";
193 
194  up_time_it--;
195 
196 // G4cout << up_time_it->first << G4endl;
197 
198  fpLastSearch->fLowerBoundTime = up_time_it;
199  fpLastSearch->fLowerBoundSet = true;
200 
201 // G4cout << "returning : " << fpLastSearch->fLowerBoundTime->second << G4endl;
202 
203  return fpLastSearch->fLowerBoundTime->second;
204 }
205 
206 //------------------------------------------------------------------------------
207 
209  double time)
210 {
211  G4bool sameTypeOfMolecule = SearchTimeMap(molecule);
212  return SearchUpperBoundTime(time, sameTypeOfMolecule);
213 }
214 
215 //------------------------------------------------------------------------------
216 
218  G4double time,
219  const G4ThreeVector* /*position*/,
220  int number)
221 {
222  if (fDontRegister[molecule->GetDefinition()]) return;
223 
224  if (fVerbose){
225  G4cout << "G4MoleculeCounter::AddAMoleculeAtTime : " << molecule->GetName()
226  << " at time : " << G4BestUnit(time, "Time") << G4endl;
227  }
228 
229  CounterMapType::iterator counterMap_i =
230  fCounterMap.find(molecule);
231 
232  if (counterMap_i == fCounterMap.end()){
233  fCounterMap[molecule][time] = number;
234  }
235  else if (counterMap_i->second.empty()){
236  counterMap_i->second[time] = number;
237  }
238  else{
239  NbMoleculeAgainstTime::reverse_iterator end = counterMap_i->second.rbegin();
240 
241  if (end->first <= time ||
242  fabs(end->first - time) <= compDoubleWithPrecision::fPrecision)
243  // Case 1 = new time comes after last recorded data
244  // Case 2 = new time is about the same as the last recorded one
245  {
246  double newValue = end->second + number;
247  counterMap_i->second[time] = newValue;
248  }
249  else
250  {
251 // if(fabs(time - G4Scheduler::Instance()->GetGlobalTime()) >
252 // G4Scheduler::Instance()->GetTimeTolerance())
253  {
254  G4ExceptionDescription errMsg;
255  errMsg << "Time of species "
256  << molecule->GetName() << " is "
257  << G4BestUnit(time, "Time") << " while "
258  << " global time is "
259  << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
260  << G4endl;
261  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
262  "TIME_DONT_MATCH",
263  FatalException, errMsg);
264  }
265  }
266  }
267 }
268 
269 //------------------------------------------------------------------------------
270 
271 void
273  G4double time,
274  const G4ThreeVector* /*position*/,
275  int number)
276 {
277  if (fDontRegister[molecule->GetDefinition()]) return;
278 
279  if (fVerbose)
280  {
281  G4cout << "G4MoleculeCounter::RemoveAMoleculeAtTime : "
282  << molecule->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  << molecule->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[molecule];
305 
306  if (nbMolPerTime.empty())
307  {
308  molecule->PrintState();
309  Dump();
310  G4String errMsg =
311  "You are trying to remove molecule " + molecule->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  else
319  {
320  NbMoleculeAgainstTime::reverse_iterator it = nbMolPerTime.rbegin();
321 
322  if (it == nbMolPerTime.rend()){
323  it--;
324 
325  G4String errMsg = "There was no " + molecule->GetName()
326  + " recorded at the time or even before the time asked";
327  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime", "",
328  FatalErrorInArgument, errMsg);
329  }
330 
331  if (time - it->first < -compDoubleWithPrecision::fPrecision){
332  Dump();
333  G4ExceptionDescription errMsg;
334  errMsg << "Is time going back?? " << molecule->GetName()
335  << " is being removed at time " << G4BestUnit(time, "Time")
336  << " while last recorded time was "
337  << G4BestUnit(it->first, "Time") << ".";
338  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
339  "RETURN_TO_THE_FUTUR",
341  errMsg);
342  }
343 
344  double finalN = it->second - number;
345 
346  if(finalN < 0){
347  Dump();
348  G4ExceptionDescription errMsg;
349  errMsg << "After removal of " << number << " species of "
350  << molecule->GetName() << " the final number at time "
351  << G4BestUnit(time, "Time") << " is less than zero and so not valid."
352  << " Global time is "
353  << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
354  << ". Previous selected time is "
355  << G4BestUnit(it->first, "Time")
356  << G4endl;
357  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
358  "N_INF_0",
359  FatalException, errMsg);
360  }
361 
362  nbMolPerTime[time] = finalN;
363  }
364 }
365 
366 //------------------------------------------------------------------------------
367 
369 {
370  if (fVerbose > 1)
371  {
372  G4cout << "Entering in G4MoleculeCounter::RecordMolecules" << G4endl;
373  }
374 
375  CounterMapType::iterator it;
376  RecordedMolecules output (new vector<G4MolecularConfiguration*>);
377 
378  for(it = fCounterMap.begin(); it != fCounterMap.end(); it++)
379  {
380  output->push_back(it->first);
381  }
382  return output;
383 }
384 
385 //------------------------------------------------------------------------------
386 
388 {
389  RecordedTimes output(new std::set<G4double>);
390 
391  //G4double time;
392 
393  CounterMapType::iterator it;
394  CounterMapType::const_iterator ite;
395 
396  NbMoleculeAgainstTime::iterator it2;
397  NbMoleculeAgainstTime::const_iterator ite2;
398 
399  // iterate on each molecule
400  for (it = fCounterMap.begin(), ite = fCounterMap.end(); it != ite; ++it)
401  {
402  // iterate on each time
403  for (it2 = (it->second).begin(), ite2 = (it->second).end(); it2 != ite2;
404  ++it2)
405  {
406  //time = it2->first;
407  output->insert(it2->first);
408  }
409  }
410 
411  return output;
412 }
413 
414 //------------------------------------------------------------------------------
415 
416 // >>DEV<<
417 //void G4MoleculeCounter::SignalReceiver(G4SpeciesInCM* /*speciesInCM*/,
418 // size_t moleculeID,
419 // int /*number*/,
420 // G4SpeciesInCM::SpeciesChange speciesChange,
421 // int diff)
422 //{
423 // switch(speciesChange)
424 // {
425 // case G4SpeciesInCM::eAdd:
426 // AddAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
427 // G4Scheduler::Instance()->GetGlobalTime(),
428 // diff);
429 // break;
430 // case G4SpeciesInCM::eRemove:
431 // RemoveAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
432 // G4Scheduler::Instance()->GetGlobalTime(),
433 // diff);
434 // break;
435 // }
436 //}
437 
438 //------------------------------------------------------------------------------
439 
441 {
442  CounterMapType::iterator it = fCounterMap.begin();
443  CounterMapType::iterator end = fCounterMap.end();
444 
445  for(;it!=end;++it)
446  {
447  G4MolecularConfiguration* molConf = it->first;
448 
449  G4cout << " --- > For " << molConf->GetName() << G4endl;
450  NbMoleculeAgainstTime::iterator it2 = it->second.begin();
451  NbMoleculeAgainstTime::iterator end2 = it->second.end();
452 
453  for(;it2!=end2;++it2)
454  {
455  G4cout << " " << G4BestUnit(it2->first, "Time")
456  << " " << it2->second << G4endl;
457  }
458  }
459 }
void AddAMoleculeAtTime(G4MolecularConfiguration *, G4double time, const G4ThreeVector *position=nullptr, int number=1) override
int SearchUpperBoundTime(double time, bool sameTypeOfMolecule)
const G4String & GetName() const
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
#define G4ThreadLocal
Definition: tls.hh:89
std::unique_ptr< std::vector< G4MolecularConfiguration * > > RecordedMolecules
void RemoveAMoleculeAtTime(G4MolecularConfiguration *, G4double time, const G4ThreeVector *position=nullptr, int number=1) override
static G4Scheduler * Instance()
Definition: G4Scheduler.cc:102
static constexpr double picosecond
Definition: G4SIunits.hh:161
G4bool SearchTimeMap(G4MolecularConfiguration *molecule)
int GetNMoleculesAtTime(G4MolecularConfiguration *molecule, double time)
RecordedMolecules GetRecordedMolecules()
G4GLOB_DLL std::ostream G4cout
static G4MoleculeCounter * Instance()
bool G4bool
Definition: G4Types.hh:79
void Initialize() override
static G4MoleculeTable * Instance()
G4ConfigurationIterator GetConfigurationIterator()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::unique_ptr< std::set< G4double > > RecordedTimes
static G4ThreadLocal double fPrecision
const G4MoleculeDefinition * GetDefinition() const
std::map< G4double, G4int, compDoubleWithPrecision > NbMoleculeAgainstTime
#define G4endl
Definition: G4ios.hh:61
RecordedTimes GetRecordedTimes()
double G4double
Definition: G4Types.hh:76
virtual ~G4MoleculeCounter()
void SetTimeSlice(double)