Geant4  10.02.p01
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 94423 2015-11-16 08:25:40Z gcosmo $
27 //
28 #include "G4MoleculeCounter.hh"
29 #include "G4MoleculeTable.hh"
30 #include "G4UIcommand.hh"
31 #include "G4UnitsTable.hh"
33 #include "G4MoleculeDefinition.hh"
34 #include <iomanip>
35 #include "G4SystemOfUnits.hh"
36 #include "G4Scheduler.hh" // TODO: remove this dependency
37 
38 using namespace std;
39 
43 
45 {
46  fUse=flag;
47 }
48 
50 {
51  return fUse;
52 }
53 
54 
56 {
57  fVerbose = 0;
58  fCheckTimeIsConsistentWithScheduler = true;
60  {
62  }
63 }
64 
66 {
67 }
68 
70 {
71  if (!fpInstance) fpInstance = new G4MoleculeCounter();
72 
73  return fpInstance;
74 }
75 
77 {
78  if (!fpInstance) fpInstance = new G4MoleculeCounter();
79 
80  return fpInstance;
81 }
82 
84 {
85  if (fpInstance)
86  {
87  delete fpInstance;
88  fpInstance = 0;
89  }
90 }
91 
93 {
94  if(fpInstance) fpInstance->Initialize();
95 }
96 
98 {
101  while ((mol_iterator)())
102  {
103  // G4cout << "G4MoleculeCounter::Initialize" << G4endl;
104  // G4cout << mol_iterator->value()->GetName() << G4endl;
105  fCounterMap[mol_iterator.value()]; // initialize the second map
106  }
107 }
108 
109 void G4MoleculeCounter::SetTimeSlice(double timeSlice)
110 {
112 }
113 
115 {
116  if (fpLastSearch.get() == 0)
117  {
118  fpLastSearch.reset(new Search());
119  }
120  else
121  {
122  if (fpLastSearch->fLowerBoundSet &&
123  fpLastSearch->fLastMoleculeSearched->first == molecule)
124  return true;
125  }
126 
127  CounterMapType::iterator mol_it = fCounterMap.find(molecule);
128  fpLastSearch->fLastMoleculeSearched = mol_it;
129 
130  if (mol_it != fCounterMap.end()) // TODO
131  {
132  fpLastSearch->fLowerBoundTime = fpLastSearch->fLastMoleculeSearched->second
133  .end();
134  fpLastSearch->fLowerBoundSet = true;
135  }
136  else
137  {
138  fpLastSearch->fLowerBoundSet = false;
139  }
140 
141  return false;
142 }
143 
145  bool sameTypeOfMolecule)
146 {
147  CounterMapType::iterator mol_it = fpLastSearch->fLastMoleculeSearched;
148  if (mol_it == fCounterMap.end()) return 0; // RETURN
149 
150  NbMoleculeAgainstTime& timeMap = mol_it->second;
151  if (timeMap.empty()) return 0;
152 
153  NbMoleculeAgainstTime::iterator end_time = timeMap.end();
154 
155  if (sameTypeOfMolecule == true)
156  {
157  //G4cout << "SAME MOLECULE" << G4endl;
158  if (fpLastSearch->fLowerBoundSet && fpLastSearch->fLowerBoundTime
159  != end_time)
160  {
161  //G4cout << fpLastSearch->fLowerBoundTime->first << G4endl;
162 // G4cout << "fpLastSearch->fLowerBoundTime != timeMap.end() " << time << G4endl;
163  if (fpLastSearch->fLowerBoundTime->first < time)
164  {
165  NbMoleculeAgainstTime::iterator upperToLast = fpLastSearch
166  ->fLowerBoundTime;
167  upperToLast++;
168 
169  if (upperToLast == end_time)
170  {
171  return fpLastSearch->fLowerBoundTime->second;
172  }
173 
174  if (upperToLast->first > time)
175  {
176  return fpLastSearch->fLowerBoundTime->second;
177  }
178  }
179  }
180  }
181  /*
182  else
183  {
184  G4cout << "\n" << G4endl;
185  G4cout << "Molecule has changed" << G4endl;
186  G4cout << "\n" << G4endl;
187  }
188  */
189  //G4cout << "Searching" << G4endl;
190  // With upper bound
191  NbMoleculeAgainstTime::iterator up_time_it = timeMap.upper_bound(time);
192 
193  if (up_time_it == end_time)
194  {
195  NbMoleculeAgainstTime::reverse_iterator last_time = timeMap.rbegin();
196 
197 // if (last_time->first <= time)
198 // {
199  //G4cout << "RETURN LAST : " << G4BestUnit(time, "Time") << G4endl;
200  return last_time->second;
201 // }
202 
203 // G4cout << "RETURN 0 (1)" << G4endl;
204 // return 0; // RETURN
205  }
206  if (up_time_it == timeMap.begin())
207  {
208 // G4cout << "RETURN 0 (2)" << G4endl;
209  return 0; // RETURN
210  }
211 
212  //G4cout << "Going back : " << up_time_it->first << "-->";
213 
214  up_time_it--;
215 
216 // G4cout << up_time_it->first << G4endl;
217 
218  fpLastSearch->fLowerBoundTime = up_time_it;
219  fpLastSearch->fLowerBoundSet = true;
220 
221 // G4cout << "returning : " << fpLastSearch->fLowerBoundTime->second << G4endl;
222 
223  return fpLastSearch->fLowerBoundTime->second;
224 }
225 
227  double time)
228 {
229  G4bool sameTypeOfMolecule = SearchTimeMap(molecule);
230  return SearchUpperBoundTime(time, sameTypeOfMolecule);
231  // NbMoleculeAgainstTime::iterator low_time_it = timeMap.lower_bound(time);
232 }
233 
235  G4double time, int number)
236 {
237  if (fDontRegister[molecule->GetDefinition()]) return;
238 
239  if (fVerbose)
240  {
241  G4cout << "G4MoleculeCounter::AddAMoleculeAtTime : " << molecule->GetName()
242  << " at time : " << G4BestUnit(time, "Time") << G4endl;
243  }
244 
245  CounterMapType::iterator counterMap_i =
246  fCounterMap.find(molecule);
247 
248  if (counterMap_i == fCounterMap.end())
249  {
250  // DEBUG
251  // if(fVerbose) G4cout << " !! ***** Map is empty " << G4endl;
252  fCounterMap[molecule][time] = number;
253  }
254  else if (counterMap_i->second.empty())
255  {
256  // DEBUG
257  // if(fVerbose) G4cout << " !! ***** Map is empty " << G4endl;
258  counterMap_i->second[time] = number;
259 // G4cout << " !! New N = " << number << G4endl;
260  }
261  else
262  {
263  NbMoleculeAgainstTime::reverse_iterator end = counterMap_i->second.rbegin();
264 // end--;
265 
266  // DEBUG
267  // if(fVerbose)
268  // G4cout<<"!! End Time = "<< G4BestUnit(end->first, "Time") <<G4endl;
269 
270  if (end->first <= time ||
271  fabs(end->first - time) <= compDoubleWithPrecision::fPrecision)
272  // Case 1 = new time comes after last recorded data
273  // Case 2 = new time is about the same as the last recorded one
274  {
275  double newValue = end->second + number;
276 // G4cout << " !! New N = " << newValue << G4endl;
277  counterMap_i->second[time] = newValue;
278 // double newValue = end->second + number;
279 // G4cout << " !! New N = " << end->second + number << G4endl;
280 // counterMap_i->second[time] = end->second + number;
281  }
282  else
283  {
284 // if(fabs(time - G4Scheduler::Instance()->GetGlobalTime()) >
285 // G4Scheduler::Instance()->GetTimeTolerance())
286  {
287  G4ExceptionDescription errMsg;
288  errMsg << "Time of species "
289  << molecule->GetName() << " is "
290  << G4BestUnit(time, "Time") << " while "
291  << " global time is "
292  << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
293  << G4endl;
294  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
295  "TIME_DONT_MATCH",
296  FatalException, errMsg);
297  }
298 
299 
300 // NbMoleculeAgainstTime::iterator it = counterMap_i->second.lower_bound(
301 // time);
302 //
303 // while (it->first > time && it != counterMap_i->second.begin())
304 // {
305 // // DEBUG
306 // // if(fVerbose)
307 // // G4cout<<"!! ********** Is going back!!!!"<<G4endl;
308 // it--;
309 // }
310 //
311 // if (it == counterMap_i->second.begin() && it->first > time)
312 // {
313 // // DEBUG
314 // // if(fVerbose)
315 // // G4cout<<"!! ********** Illegal !!!!"<<G4endl;
316 // return;
317 // }
318 //
319 // // DEBUG
320 // // if(fVerbose)
321 // // {
322 // // G4cout<<"!! PREVIOUS NB = "<< it->second <<G4endl;
323 // // G4cout<<"!! PREVIOUS TIME = "<< G4BestUnit(it->first,"Time") <<G4endl;
324 // // }
325 // counterMap_i->second[time] = it->second + number;
326  }
327  }
328 
329  // DEBUG
330  // if(fVerbose)
331  // G4cout<<"!! NB = "<< fCounterMap[molecule][time]<<G4endl;
332 }
333 
335  G4double time, int number)
336 {
337  if (fDontRegister[molecule->GetDefinition()]) return;
338 
339  if (fVerbose)
340  {
341  G4cout << "G4MoleculeCounter::RemoveAMoleculeAtTime : "
342  << molecule->GetName() << " at time : " << G4BestUnit(time, "Time")
343  << G4endl;
344  }
345 
346  if(fCheckTimeIsConsistentWithScheduler)
347  {
348  if(fabs(time - G4Scheduler::Instance()->GetGlobalTime()) >
349  G4Scheduler::Instance()->GetTimeTolerance())
350  {
351  G4ExceptionDescription errMsg;
352  errMsg << "Time of species "
353  << molecule->GetName() << " is "
354  << G4BestUnit(time, "Time") << " while "
355  << " global time is "
356  << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
357  << G4endl;
358  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
359  "TIME_DONT_MATCH",
360  FatalException, errMsg);
361  }
362  }
363 
364  NbMoleculeAgainstTime& nbMolPerTime = fCounterMap[molecule];
365 
366  if (nbMolPerTime.empty())
367  {
368  molecule->PrintState();
369  Dump();
370  G4String errMsg =
371  "You are trying to remove molecule " + molecule->GetName()
372  + " from the counter while this kind of molecules has not been registered yet";
373  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime", "",
374  FatalErrorInArgument, errMsg);
375 
376  return;
377  }
378  else
379  {
380  NbMoleculeAgainstTime::reverse_iterator it = nbMolPerTime.rbegin();
381 
382 // if (nbMolPerTime.size() == 1)
383 // {
384 // it = nbMolPerTime.begin();
385 // // DEBUG
386 // // if(fVerbose)
387 // // G4cout << "!! fCounterMap[molecule].size() == 1" << G4endl;
388 // }
389 // else
390 // {
392 // it = nbMolPerTime.end();
393 // --it;
394 // }
395 
396  if (it == nbMolPerTime.rend())
397  {
398  // DEBUG
399  // if(fVerbose)
400  // G4cout << " ********** NO ITERATOR !!!!!!!!! " << G4endl;
401  it--;
402 
403  //if (time < it->first)
404  {
405  G4String errMsg = "There was no " + molecule->GetName()
406  + " recorded at the time or even before the time asked";
407  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime", "",
408  FatalErrorInArgument, errMsg);
409  }
410  }
411 
412  // DEBUG
413  // if(fVerbose)
414  // {
416  // G4cout<<"!! Molecule = " << molecule.GetName() << G4endl;
417  // G4cout<<"!! At Time = "<< G4BestUnit(time,"Time") <<G4endl;
418  // G4cout<<"!! PREVIOUS TIME = "<< G4BestUnit(it->first,"Time")<<G4endl;
419  // G4cout<<"!! PREVIOUS Nb = "<< it->second <<G4endl;
420  // }
421 
422  // If valgrind problem on the line below, it means that the pointer "it"
423  // points nowhere
424 // if (nbMolPerTime.value_comp()(*it, *nbMolPerTime.begin()))
425 // {
426 // // DEBUG
427 // // if(fVerbose)
428 // // G4cout<<"!! ***** In value_comp ... " << G4endl;
429 // it++;
430 // if (time < it->first)
431 // {
432 // G4String errMsg = "There was no " + molecule->GetName()
433 // + " record at the time or even before the time asked";
434 // G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime", "",
435 // FatalErrorInArgument, errMsg);
436 // }
437 // }
438 //
439 // while (it->first - time > compDoubleWithPrecision::fPrecision
440 // && it != nbMolPerTime.begin())
441 // {
442 // // DEBUG
443 // // if(fVerbose)
444 // // {
445 // // G4cout<<"!! ***** Is going back!!!!"<<G4endl;
446 // // G4cout<<"!! PREVIOUS TIME = "<< G4BestUnit(it-> first,"Time") <<G4endl;
447 // // }
448 // it--;
449 // }
450 
451  if (time - it->first < -compDoubleWithPrecision::fPrecision)
452  {
453  Dump();
454  G4ExceptionDescription errMsg;
455  errMsg << "Is time going back?? " << molecule->GetName()
456  << " is being removed at time " << G4BestUnit(time, "Time")
457  << " while last recorded time was "
458  << G4BestUnit(it->first, "Time") << ".";
459  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
460  "RETURN_TO_THE_FUTUR",
462  errMsg);
463  }
464 /*
465  if (it == nbMolPerTime.begin() && it->first > time)
466  {
467  // DEBUG
468  // if(fVerbose)
469  // G4cout<<"!! ********** Illegal !!!!"<<G4endl;
470  return;
471  }
472 */
473  // DEBUG
474  // if(fVerbose)
475  // {
476  // G4cout<<"!! PREVIOUS NB = "<< (*it).second <<G4endl;
477  // G4cout<<"!! PREVIOUS TIME = "<< G4BestUnit(it->first,"Time")<<G4endl;
478  // }
479 
480  double finalN = it->second - number;
481 
482  if(finalN < 0)
483  {
484  Dump();
485  G4ExceptionDescription errMsg;
486  errMsg << "After removal of " << number << " species of "
487  << molecule->GetName() << " the final number at time "
488  << G4BestUnit(time, "Time") << " is less than zero and so not valid."
489  << " Global time is "
490  << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
491  << ". Previous selected time is "
492  << G4BestUnit(it->first, "Time")
493  << G4endl;
494  G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
495  "N_INF_0",
496  FatalException, errMsg);
497  }
498 
499  nbMolPerTime[time] = finalN;
500  }
501 
502  // DEBUG
503  // if(fVerbose)
504  // {
505  // G4cout<<"!! NB = "<< nbMolPerTime[time]<<G4endl;
506  // }
507 }
508 
510 {
511  if (fVerbose > 1)
512  {
513  G4cout << "Entering in G4MoleculeCounter::RecordMolecules" << G4endl;
514  }
515 
516  CounterMapType::iterator it;
517  RecordedMolecules output (new vector<G4MolecularConfiguration*>);
518 
519  for(it = fCounterMap.begin(); it != fCounterMap.end(); it++)
520  {
521  output->push_back(it->first);
522  }
523  return output;
524 }
525 
527 {
528  RecordedTimes output(new std::set<G4double>);
529 
530  //G4double time;
531 
532  CounterMapType::iterator it;
533  CounterMapType::const_iterator ite;
534 
535  NbMoleculeAgainstTime::iterator it2;
536  NbMoleculeAgainstTime::const_iterator ite2;
537 
538  // iterate on each molecule
539  for (it = fCounterMap.begin(), ite = fCounterMap.end(); it != ite; ++it)
540  {
541  // iterate on each time
542  for (it2 = (it->second).begin(), ite2 = (it->second).end(); it2 != ite2;
543  ++it2)
544  {
545  //time = it2->first;
546  output->insert(it2->first);
547  }
548  }
549 
550  return output;
551 }
552 
553 // >>DEV<<
554 //void G4MoleculeCounter::SignalReceiver(G4SpeciesInCM* /*speciesInCM*/,
555 // size_t moleculeID,
556 // int /*number*/,
557 // G4SpeciesInCM::SpeciesChange speciesChange,
558 // int diff)
559 //{
560 // switch(speciesChange)
561 // {
562 // case G4SpeciesInCM::eAdd:
563 // AddAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
564 // G4Scheduler::Instance()->GetGlobalTime(),
565 // diff);
566 // break;
567 // case G4SpeciesInCM::eRemove:
568 // RemoveAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
569 // G4Scheduler::Instance()->GetGlobalTime(),
570 // diff);
571 // break;
572 // }
573 //}
574 
575 
577 {
578  CounterMapType::iterator it = fCounterMap.begin();
579  CounterMapType::iterator end = fCounterMap.end();
580 
581  for(;it!=end;++it)
582  {
583  G4MolecularConfiguration* molConf = it->first;
584 
585  G4cout << " --- > For " << molConf->GetName() << G4endl;
586  NbMoleculeAgainstTime::iterator it2 = it->second.begin();
587  NbMoleculeAgainstTime::iterator end2 = it->second.end();
588 
589  for(;it2!=end2;++it2)
590  {
591  G4cout << " "<< G4BestUnit(it2->first, "Time") << " " << it2->second << G4endl;
592  }
593  }
594 }
The pointer G4MolecularConfiguration will be shared by all the molecules having the same molecule def...
void PrintState() const
Display the electronic state of the molecule.
int SearchUpperBoundTime(double time, bool sameTypeOfMolecule)
const G4String & GetName() const
Returns the name of the molecule.
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
static G4MoleculeCounter * GetMoleculeCounter()
virtual void RemoveAMoleculeAtTime(G4MolecularConfiguration *, G4double time, int number=1)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
#define G4ThreadLocal
Definition: tls.hh:89
stdunique_ptr< std::vector< G4MolecularConfiguration * > > RecordedMolecules
static G4Scheduler * Instance()
Definition: G4Scheduler.cc:102
G4bool SearchTimeMap(G4MolecularConfiguration *molecule)
int GetNMoleculesAtTime(G4MolecularConfiguration *molecule, double time)
RecordedMolecules GetRecordedMolecules()
G4GLOB_DLL std::ostream G4cout
static G4bool InUse()
static G4MoleculeCounter * Instance()
bool G4bool
Definition: G4Types.hh:79
static G4MoleculeTable * Instance()
stdunique_ptr< std::set< G4double > > RecordedTimes
static void Use(G4bool flag=true)
virtual void AddAMoleculeAtTime(G4MolecularConfiguration *, G4double time, int number=1)
G4ConfigurationIterator GetConfigurationIterator()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4ThreadLocal double fPrecision
static G4ThreadLocal G4MoleculeCounter * fpInstance
const G4MoleculeDefinition * GetDefinition() const
std::map< G4double, G4int, compDoubleWithPrecision > NbMoleculeAgainstTime
static void InitializeInstance()
#define G4endl
Definition: G4ios.hh:61
RecordedTimes GetRecordedTimes()
static const double picosecond
Definition: G4SIunits.hh:160
double G4double
Definition: G4Types.hh:76
virtual ~G4MoleculeCounter()
void SetTimeSlice(double)
static void DeleteInstance()