| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4MoleculeCounter.cc 94423 2015-11-16 08:25:40Z gcosmo $
26 // 27 //
27 <<
28 #include "G4MoleculeCounter.hh" 28 #include "G4MoleculeCounter.hh"
29 <<
30 #include "G4MolecularConfiguration.hh" <<
31 #include "G4MoleculeDefinition.hh" <<
32 #include "G4MoleculeTable.hh" 29 #include "G4MoleculeTable.hh"
33 #include "G4Scheduler.hh" // TODO: remove this <<
34 #include "G4SystemOfUnits.hh" <<
35 #include "G4UIcommand.hh" 30 #include "G4UIcommand.hh"
36 #include "G4UnitsTable.hh" 31 #include "G4UnitsTable.hh"
37 << 32 #include "G4MolecularConfiguration.hh"
>> 33 #include "G4MoleculeDefinition.hh"
38 #include <iomanip> 34 #include <iomanip>
39 #include <memory> << 35 #include "G4SystemOfUnits.hh"
>> 36 #include "G4Scheduler.hh" // TODO: remove this dependency
40 37
41 using namespace std; 38 using namespace std;
42 39
43 namespace G4{ << 40 G4ThreadLocal double compDoubleWithPrecision::fPrecision = 0;
44 namespace MoleculeCounter { << 41 G4bool G4MoleculeCounter::fUse = false;
>> 42 G4ThreadLocal G4MoleculeCounter* G4MoleculeCounter::fpInstance = 0;
45 43
46 bool TimePrecision::operator()(const double& a << 44 void G4MoleculeCounter::Use(G4bool flag)
47 { 45 {
48 if (std::fabs(a - b) < fPrecision) << 46 fUse=flag;
49 { <<
50 return false; <<
51 } <<
52 <<
53 return a < b; <<
54 } <<
55 <<
56 G4ThreadLocal double TimePrecision::fPrecision <<
57 } <<
58 } 47 }
59 48
60 //-------------------------------------------- << 49 G4bool G4MoleculeCounter::InUse()
61 G4MoleculeCounter* G4MoleculeCounter::Instance <<
62 { 50 {
63 if (fpInstance == nullptr) << 51 return fUse;
64 { <<
65 fpInstance = new G4MoleculeCounter(); <<
66 } <<
67 return dynamic_cast<G4MoleculeCounter*>(fp <<
68 } 52 }
69 53
70 //-------------------------------------------- <<
71 54
72 G4MoleculeCounter::G4MoleculeCounter() 55 G4MoleculeCounter::G4MoleculeCounter()
73 { 56 {
74 fVerbose = 0; << 57 fVerbose = 0;
75 fCheckTimeIsConsistentWithScheduler = true << 58 fCheckTimeIsConsistentWithScheduler = true;
>> 59 if(compDoubleWithPrecision::fPrecision == 0)
>> 60 {
>> 61 compDoubleWithPrecision::fPrecision = 0.5*picosecond;
>> 62 }
76 } 63 }
77 64
78 //-------------------------------------------- << 65 G4MoleculeCounter::~G4MoleculeCounter()
>> 66 {
>> 67 }
79 68
80 G4MoleculeCounter::~G4MoleculeCounter() = defa << 69 G4MoleculeCounter* G4MoleculeCounter::GetMoleculeCounter()
>> 70 {
>> 71 if (!fpInstance) fpInstance = new G4MoleculeCounter();
81 72
82 //-------------------------------------------- << 73 return fpInstance;
>> 74 }
83 75
84 void G4MoleculeCounter::Initialize() << 76 G4MoleculeCounter* G4MoleculeCounter::Instance()
85 { 77 {
86 auto mol_iterator = G4MoleculeTable::Insta << 78 if (!fpInstance) fpInstance = new G4MoleculeCounter();
87 while ((mol_iterator)()) <<
88 { <<
89 if (!IsRegistered(mol_iterator.value() <<
90 { <<
91 continue; <<
92 } <<
93 79
94 fCounterMap[mol_iterator.value()]; // << 80 return fpInstance;
95 } <<
96 } 81 }
97 82
98 //-------------------------------------------- << 83 void G4MoleculeCounter::DeleteInstance()
99 <<
100 void G4MoleculeCounter::SetTimeSlice(double ti <<
101 { 84 {
102 G4::MoleculeCounter::TimePrecision::fPreci << 85 if (fpInstance)
>> 86 {
>> 87 delete fpInstance;
>> 88 fpInstance = 0;
>> 89 }
103 } 90 }
104 91
105 //-------------------------------------------- << 92 void G4MoleculeCounter::InitializeInstance()
106 <<
107 G4bool G4MoleculeCounter::SearchTimeMap(Reacta <<
108 { 93 {
109 if (fpLastSearch == nullptr) << 94 if(fpInstance) fpInstance->Initialize();
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 } 95 }
138 96
139 //-------------------------------------------- << 97 void G4MoleculeCounter::Initialize()
140 <<
141 int G4MoleculeCounter::SearchUpperBoundTime(do <<
142 bo <<
143 { 98 {
144 auto mol_it = fpLastSearch->fLastMoleculeS << 99 G4ConfigurationIterator mol_iterator = G4MoleculeTable::Instance()
145 if (mol_it == fCounterMap.end()) << 100 ->GetConfigurationIterator();
146 { << 101 while ((mol_iterator)())
147 return 0; << 102 {
148 } << 103 // G4cout << "G4MoleculeCounter::Initialize" << G4endl;
149 << 104 // G4cout << mol_iterator->value()->GetName() << G4endl;
150 NbMoleculeAgainstTime& timeMap = mol_it->s << 105 fCounterMap[mol_iterator.value()]; // initialize the second map
151 if (timeMap.empty()) << 106 }
152 { << 107 }
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 108
190 up_time_it--; << 109 void G4MoleculeCounter::SetTimeSlice(double timeSlice)
>> 110 {
>> 111 compDoubleWithPrecision::fPrecision = timeSlice;
>> 112 }
191 113
192 fpLastSearch->fLowerBoundTime = up_time_it << 114 G4bool G4MoleculeCounter::SearchTimeMap(G4MolecularConfiguration* molecule)
>> 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();
193 fpLastSearch->fLowerBoundSet = true; 134 fpLastSearch->fLowerBoundSet = true;
>> 135 }
>> 136 else
>> 137 {
>> 138 fpLastSearch->fLowerBoundSet = false;
>> 139 }
194 140
195 return fpLastSearch->fLowerBoundTime->seco << 141 return false;
196 } 142 }
197 143
198 //-------------------------------------------- << 144 int G4MoleculeCounter::SearchUpperBoundTime(double time,
199 << 145 bool sameTypeOfMolecule)
200 int G4MoleculeCounter::GetNMoleculesAtTime(Rea <<
201 dou <<
202 { 146 {
203 G4bool sameTypeOfMolecule = SearchTimeMap( << 147 CounterMapType::iterator mol_it = fpLastSearch->fLastMoleculeSearched;
204 return SearchUpperBoundTime(time, sameType << 148 if (mol_it == fCounterMap.end()) return 0; // RETURN
205 } <<
206 <<
207 //-------------------------------------------- <<
208 149
209 void G4MoleculeCounter::AddAMoleculeAtTime(Rea << 150 NbMoleculeAgainstTime& timeMap = mol_it->second;
210 G4d << 151 if (timeMap.empty()) return 0;
211 con <<
212 int <<
213 { <<
214 if (fDontRegister[molecule->GetDefinition( <<
215 { <<
216 return; <<
217 } <<
218 152
219 if (fVerbose != 0) << 153 NbMoleculeAgainstTime::iterator end_time = timeMap.end();
220 { <<
221 G4cout << "G4MoleculeCounter::AddAMole <<
222 << " at time : " << G4BestUnit( <<
223 } <<
224 154
225 auto counterMap_i = fCounterMap.find(molec << 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++;
226 168
227 if (counterMap_i == fCounterMap.end()) << 169 if (upperToLast == end_time)
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 { 170 {
244 double newValue = end->second + nu << 171 return fpLastSearch->fLowerBoundTime->second;
245 counterMap_i->second[time] = newVa <<
246 } 172 }
247 else << 173
>> 174 if (upperToLast->first > time)
248 { 175 {
249 // if(fabs(time - G4Scheduler << 176 return fpLastSearch->fLowerBoundTime->second;
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 } 177 }
>> 178 }
264 } 179 }
265 } << 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;
266 220
267 //-------------------------------------------- << 221 // G4cout << "returning : " << fpLastSearch->fLowerBoundTime->second << G4endl;
268 222
269 void G4MoleculeCounter::RemoveAMoleculeAtTime( << 223 return fpLastSearch->fLowerBoundTime->second;
270 << 224 }
271 <<
272 <<
273 { <<
274 if (fDontRegister[pMolecule->GetDefinition <<
275 { <<
276 return; <<
277 } <<
278 225
279 if (fVerbose != 0) << 226 int G4MoleculeCounter::GetNMoleculesAtTime(G4MolecularConfiguration* molecule,
280 { << 227 double time)
281 G4cout << "G4MoleculeCounter::RemoveAM << 228 {
282 << pMolecule->GetName() << " at << 229 G4bool sameTypeOfMolecule = SearchTimeMap(molecule);
283 << G4endl; << 230 return SearchUpperBoundTime(time, sameTypeOfMolecule);
>> 231 // NbMoleculeAgainstTime::iterator low_time_it = timeMap.lower_bound(time);
>> 232 }
>> 233
>> 234 void G4MoleculeCounter::AddAMoleculeAtTime(G4MolecularConfiguration* molecule,
>> 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;
284 } 281 }
285 << 282 else
286 if (fCheckTimeIsConsistentWithScheduler) <<
287 { 283 {
288 if (fabs(time - G4Scheduler::Instance( << 284 // if(fabs(time - G4Scheduler::Instance()->GetGlobalTime()) >
289 G4Scheduler::Instance()->GetTimeTo << 285 // G4Scheduler::Instance()->GetTimeTolerance())
290 { << 286 {
291 G4ExceptionDescription errMsg; << 287 G4ExceptionDescription errMsg;
292 errMsg << "Time of species " << 288 errMsg << "Time of species "
293 << pMolecule->GetName() << << 289 << molecule->GetName() << " is "
294 << G4BestUnit(time, "Time") << 290 << G4BestUnit(time, "Time") << " while "
295 << " global time is " << 291 << " global time is "
296 << G4BestUnit(G4Scheduler:: << 292 << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
297 << G4endl; << 293 << G4endl;
298 G4Exception("G4MoleculeCounter::Re << 294 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
299 "TIME_DONT_MATCH", << 295 "TIME_DONT_MATCH",
300 FatalException, errMsg << 296 FatalException, errMsg);
301 } << 297 }
302 } <<
303 <<
304 NbMoleculeAgainstTime& nbMolPerTime = fCou <<
305 298
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 299
316 return; << 300 // NbMoleculeAgainstTime::iterator it = counterMap_i->second.lower_bound(
317 } << 301 // time);
318 << 302 //
319 auto it = nbMolPerTime.rbegin(); << 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
>> 334 void G4MoleculeCounter::RemoveAMoleculeAtTime(G4MolecularConfiguration* molecule,
>> 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 // {
>> 391 //// it = nbMolPerTime.lower_bound(time);
>> 392 // it = nbMolPerTime.end();
>> 393 // --it;
>> 394 // }
320 395
321 if (it == nbMolPerTime.rend()) 396 if (it == nbMolPerTime.rend())
322 { 397 {
323 it--; << 398 // DEBUG
324 << 399 // if(fVerbose)
325 G4String errMsg = << 400 // G4cout << " ********** NO ITERATOR !!!!!!!!! " << G4endl;
326 "There was no " + pMolecule->G << 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";
327 G4Exception("G4MoleculeCounter::Remove 407 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime", "",
328 FatalErrorInArgument, errM 408 FatalErrorInArgument, errMsg);
>> 409 }
329 } 410 }
330 411
331 if (time - it->first < -G4::MoleculeCounte << 412 // DEBUG
332 { << 413 // if(fVerbose)
333 Dump(); << 414 // {
334 G4ExceptionDescription errMsg; << 415 //// G4cout << "G4MoleculeCounter::RemoveAMoleculeAtTime " << G4endl;
335 errMsg << "Is time going back?? " << p << 416 // G4cout<<"!! Molecule = " << molecule.GetName() << G4endl;
336 << " is being removed at time " << 417 // G4cout<<"!! At Time = "<< G4BestUnit(time,"Time") <<G4endl;
337 << " while last recorded time w << 418 // G4cout<<"!! PREVIOUS TIME = "<< G4BestUnit(it->first,"Time")<<G4endl;
338 << G4BestUnit(it->first, "Time" << 419 // G4cout<<"!! PREVIOUS Nb = "<< it->second <<G4endl;
339 G4Exception("G4MoleculeCounter::Remove << 420 // }
340 "RETURN_TO_THE_FUTUR", << 421
341 FatalErrorInArgument, << 422 // If valgrind problem on the line below, it means that the pointer "it"
342 errMsg); << 423 // points nowhere
343 } << 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",
>> 461 FatalErrorInArgument,
>> 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 // }
344 479
345 double finalN = it->second - number; 480 double finalN = it->second - number;
346 481
347 if (finalN < 0) << 482 if(finalN < 0)
348 { 483 {
349 Dump(); << 484 Dump();
350 G4ExceptionDescription errMsg; << 485 G4ExceptionDescription errMsg;
351 errMsg << "After removal of " << numbe << 486 errMsg << "After removal of " << number << " species of "
352 << pMolecule->GetName() << " th << 487 << molecule->GetName() << " the final number at time "
353 << G4BestUnit(time, "Time") << << 488 << G4BestUnit(time, "Time") << " is less than zero and so not valid."
354 << " Global time is " << 489 << " Global time is "
355 << G4BestUnit(G4Scheduler::Inst << 490 << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
356 << ". Previous selected time is << 491 << ". Previous selected time is "
357 << G4BestUnit(it->first, "Time" << 492 << G4BestUnit(it->first, "Time")
358 << G4endl; << 493 << G4endl;
359 G4Exception("G4MoleculeCounter::Remove << 494 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime",
360 "N_INF_0", << 495 "N_INF_0",
361 FatalException, errMsg); << 496 FatalException, errMsg);
362 } 497 }
363 498
364 nbMolPerTime[time] = finalN; 499 nbMolPerTime[time] = finalN;
365 } << 500 }
366 501
367 //-------------------------------------------- << 502 // DEBUG
>> 503 // if(fVerbose)
>> 504 // {
>> 505 // G4cout<<"!! NB = "<< nbMolPerTime[time]<<G4endl;
>> 506 // }
>> 507 }
368 508
369 G4MoleculeCounter::RecordedMolecules G4Molecul 509 G4MoleculeCounter::RecordedMolecules G4MoleculeCounter::GetRecordedMolecules()
370 { 510 {
371 if (fVerbose > 1) << 511 if (fVerbose > 1)
372 { << 512 {
373 G4cout << "Entering in G4MoleculeCount << 513 G4cout << "Entering in G4MoleculeCounter::RecordMolecules" << G4endl;
374 } << 514 }
375 << 515
376 RecordedMolecules output(new ReactantList( << 516 CounterMapType::iterator it;
377 << 517 RecordedMolecules output (new vector<G4MolecularConfiguration*>);
378 for (const auto & it : fCounterMap) << 518
379 { << 519 for(it = fCounterMap.begin(); it != fCounterMap.end(); it++)
380 output->push_back(it.first); << 520 {
381 } << 521 output->push_back(it->first);
382 return output; << 522 }
>> 523 return output;
383 } 524 }
384 525
385 //-------------------------------------------- <<
386 <<
387 RecordedTimes G4MoleculeCounter::GetRecordedTi 526 RecordedTimes G4MoleculeCounter::GetRecordedTimes()
388 { 527 {
389 RecordedTimes output(new std::set<G4double << 528 RecordedTimes output(new std::set<G4double>);
390 529
391 for(const auto& it : fCounterMap) << 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)
392 { 544 {
393 for(const auto& it2 : it.second) << 545 //time = it2->first;
394 { << 546 output->insert(it2->first);
395 //time = it2->first; <<
396 output->insert(it2.first); <<
397 } <<
398 } 547 }
>> 548 }
399 549
400 return output; << 550 return output;
401 } 551 }
402 552
403 //-------------------------------------------- <<
404 <<
405 // >>DEV<< 553 // >>DEV<<
406 //void G4MoleculeCounter::SignalReceiver(G4Spe 554 //void G4MoleculeCounter::SignalReceiver(G4SpeciesInCM* /*speciesInCM*/,
407 // size_ 555 // size_t moleculeID,
408 // int / 556 // int /*number*/,
409 // G4Spe 557 // G4SpeciesInCM::SpeciesChange speciesChange,
410 // int d 558 // int diff)
411 //{ 559 //{
412 // switch(speciesChange) 560 // switch(speciesChange)
413 // { 561 // {
414 // case G4SpeciesInCM::eAdd: 562 // case G4SpeciesInCM::eAdd:
415 // AddAMoleculeAtTime(G4MoleculeTable::In 563 // AddAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
416 // G4Scheduler::Instan 564 // G4Scheduler::Instance()->GetGlobalTime(),
417 // diff); 565 // diff);
418 // break; 566 // break;
419 // case G4SpeciesInCM::eRemove: 567 // case G4SpeciesInCM::eRemove:
420 // RemoveAMoleculeAtTime(G4MoleculeTable: 568 // RemoveAMoleculeAtTime(G4MoleculeTable::Instance()->GetConfiguration((int)moleculeID),
421 // G4Scheduler::Instan 569 // G4Scheduler::Instance()->GetGlobalTime(),
422 // diff); 570 // diff);
423 // break; 571 // break;
424 // } 572 // }
425 //} 573 //}
426 574
427 //-------------------------------------------- <<
428 575
429 void G4MoleculeCounter::Dump() 576 void G4MoleculeCounter::Dump()
430 { 577 {
431 for (const auto& it : fCounterMap) << 578 CounterMapType::iterator it = fCounterMap.begin();
432 { << 579 CounterMapType::iterator end = fCounterMap.end();
433 auto pReactant = it.first; <<
434 <<
435 G4cout << " --- > For " << pReactant-> <<
436 580
437 for (const auto& it2 : it.second) << 581 for(;it!=end;++it)
438 { << 582 {
439 G4cout << " " << G4BestUnit(it2.fi << 583 G4MolecularConfiguration* molConf = it->first;
440 << " " << it2.second << <<
441 } <<
442 } <<
443 } <<
444 584
445 //-------------------------------------------- << 585 G4cout << " --- > For " << molConf->GetName() << G4endl;
>> 586 NbMoleculeAgainstTime::iterator it2 = it->second.begin();
>> 587 NbMoleculeAgainstTime::iterator end2 = it->second.end();
446 588
447 void G4MoleculeCounter::ResetCounter() << 589 for(;it2!=end2;++it2)
448 { <<
449 if (fVerbose != 0) <<
450 { 590 {
451 G4cout << " ---> G4MoleculeCounter::Re << 591 G4cout << " "<< G4BestUnit(it2->first, "Time") << " " << it2->second << G4endl;
452 } 592 }
453 fCounterMap.clear(); << 593 }
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 } 594 }
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 595