| 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 // 26 //
>> 27 #include <memory>
27 #include "G4DNAScavengerMaterial.hh" 28 #include "G4DNAScavengerMaterial.hh"
28 << 29 #include "G4StateManager.hh"
29 #include "G4DNABoundingBox.hh" <<
30 #include "G4DNAMolecularMaterial.hh" 30 #include "G4DNAMolecularMaterial.hh"
31 #include "G4MolecularConfiguration.hh" <<
32 #include "G4PhysicalConstants.hh" 31 #include "G4PhysicalConstants.hh"
33 #include "G4Scheduler.hh" << 32 #include "G4MolecularConfiguration.hh"
34 #include "G4StateManager.hh" <<
35 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh"
36 #include "G4UnitsTable.hh" << 34 #include "G4DNABoundingBox.hh"
37 #include "G4VChemistryWorld.hh" 35 #include "G4VChemistryWorld.hh"
38 << 36 #include "G4Scheduler.hh"
39 #include <memory> << 37 #include "G4UnitsTable.hh"
>> 38 #include "G4MoleculeTable.hh"
40 39
41 using namespace std; 40 using namespace std;
42 41
43 //-------------------------------------------- 42 //------------------------------------------------------------------------------
44 43
45 G4DNAScavengerMaterial::G4DNAScavengerMaterial << 44 G4DNAScavengerMaterial::G4DNAScavengerMaterial(
46 : fpChemistryInfo(pChemistryInfo), fIsInitia << 45 G4VChemistryWorld* pChemistryInfo)
>> 46 : G4VScavengerMaterial()
>> 47 , fpChemistryInfo(pChemistryInfo)
>> 48 , fIsInitialized(false)
>> 49 , fCounterAgainstTime(false)
>> 50 , fVerbose(0)
47 { 51 {
48 Initialize(); 52 Initialize();
49 } 53 }
50 54
51 //-------------------------------------------- 55 //------------------------------------------------------------------------------
52 56
53 void G4DNAScavengerMaterial::Initialize() 57 void G4DNAScavengerMaterial::Initialize()
54 { 58 {
55 if (fIsInitialized) { << 59 if(fIsInitialized)
>> 60 {
56 return; 61 return;
57 } 62 }
58 63
59 if (fpChemistryInfo->size() == 0) { << 64 if(fpChemistryInfo->size() == 0)
>> 65 {
60 G4cout << "G4DNAScavengerMaterial existed 66 G4cout << "G4DNAScavengerMaterial existed but empty" << G4endl;
61 } 67 }
62 Reset(); 68 Reset();
63 fIsInitialized = true; 69 fIsInitialized = true;
64 } 70 }
65 71
66 G4double << 72 G4double G4DNAScavengerMaterial::GetNumberMoleculePerVolumeUnitForMaterialConf(
67 G4DNAScavengerMaterial::GetNumberMoleculePerVo << 73 MolType matConf) const
68 { 74 {
69 // no change these molecules 75 // no change these molecules
70 if (fH2O == matConf) { << 76 if(G4MoleculeTable::Instance()->GetConfiguration("H2O") == matConf)
>> 77 {
71 G4ExceptionDescription exceptionDescriptio 78 G4ExceptionDescription exceptionDescription;
72 exceptionDescription << "matConf : " << ma << 79 exceptionDescription << "matConf : "<<matConf->GetName();
73 G4Exception("G4DNAScavengerMaterial::GetNu << 80 G4Exception("G4DNAScavengerMaterial::GetNumberMoleculePerVolumeUnitForMaterialConf", "G4DNAScavengerMaterial001",
74 "G4DNAScavengerMaterial001", F << 81 FatalErrorInArgument, exceptionDescription);
75 } 82 }
76 83
77 auto iter = fScavengerTable.find(matConf); 84 auto iter = fScavengerTable.find(matConf);
78 if (iter == fScavengerTable.end()) { << 85 if(iter == fScavengerTable.end())
>> 86 {
>> 87 // G4cout<<matConf->GetName()<<G4endl;
>> 88 // throw std::runtime_error("this material is not existed");
79 return 0; 89 return 0;
80 } 90 }
81 << 91 else
82 if (iter->second >= 1) { << 92 {
83 return (floor)(iter->second); << 93 if(iter->second >= 1)
>> 94 {
>> 95 return (floor)(iter->second);
>> 96 }
>> 97 else
>> 98 {
>> 99 return 0;
>> 100 }
84 } 101 }
85 <<
86 return 0; <<
87 } 102 }
88 103
89 void G4DNAScavengerMaterial::ReduceNumberMolec << 104 void G4DNAScavengerMaterial::ReduceNumberMoleculePerVolumeUnitForMaterialConf(
90 << 105 MolType matConf, G4double time)
91 { 106 {
92 // no change these molecules 107 // no change these molecules
93 if (fH2O == matConf || fH3Op == matConf || << 108 if(G4MoleculeTable::Instance()->GetConfiguration("H2O") == matConf ||
94 fHOm == matConf) << 109 G4MoleculeTable::Instance()->GetConfiguration("H3Op(B)") ==
>> 110 matConf || // suppose that pH is not changed during simu
>> 111 G4MoleculeTable::Instance()->GetConfiguration("OHm(B)") == matConf)
95 { 112 {
96 // G4cout<<"moletype : "<<matConf->GetName 113 // G4cout<<"moletype : "<<matConf->GetName()<<G4endl;
97 // kobs is already counted these molecule 114 // kobs is already counted these molecule concentrations
98 return; 115 return;
99 } 116 }
100 if (!find(matConf)) // matConf must greater << 117 if(!find(matConf)) // matConf must greater than 0
101 { 118 {
102 return; 119 return;
103 } 120 }
104 fScavengerTable[matConf]--; 121 fScavengerTable[matConf]--;
105 if (fScavengerTable[matConf] < 0) // projec << 122 if(fScavengerTable[matConf] < 0) // projection
106 { 123 {
107 assert(false); 124 assert(false);
108 } 125 }
109 126
110 if (fCounterAgainstTime) { << 127 if(fCounterAgainstTime)
>> 128 {
111 RemoveAMoleculeAtTime(matConf, time); 129 RemoveAMoleculeAtTime(matConf, time);
112 } 130 }
113 } 131 }
114 132
115 void G4DNAScavengerMaterial::AddNumberMolecule << 133 void G4DNAScavengerMaterial::AddNumberMoleculePerVolumeUnitForMaterialConf(
116 << 134 MolType matConf, G4double time)
117 { 135 {
118 // no change these molecules 136 // no change these molecules
119 137
120 if (fH2O == matConf || fH3Op == matConf || << 138 if(G4MoleculeTable::Instance()->GetConfiguration("H2O") == matConf ||
121 G4MoleculeTable::Instance()->GetConfigur << 139 G4MoleculeTable::Instance()->GetConfiguration("H3Op(B)") ==
>> 140 matConf || // pH has no change
>> 141 G4MoleculeTable::Instance()->GetConfiguration("OHm(B)") == matConf)
122 { 142 {
123 // G4cout<<"moletype : "<<matConf->GetName 143 // G4cout<<"moletype : "<<matConf->GetName()<<G4endl;
124 // kobs is already counted these molecule 144 // kobs is already counted these molecule concentrations
125 return; 145 return;
126 } 146 }
127 147
128 auto it = fScavengerTable.find(matConf); 148 auto it = fScavengerTable.find(matConf);
129 if (it == fScavengerTable.end()) // matConf << 149 if(it == fScavengerTable.end()) // matConf must be in fScavengerTable
130 { 150 {
131 return; 151 return;
132 } 152 }
133 fScavengerTable[matConf]++; 153 fScavengerTable[matConf]++;
134 154
135 if (fCounterAgainstTime) { << 155 if(fCounterAgainstTime)
>> 156 {
136 AddAMoleculeAtTime(matConf, time); 157 AddAMoleculeAtTime(matConf, time);
137 } 158 }
138 } 159 }
139 160
140 void G4DNAScavengerMaterial::PrintInfo() 161 void G4DNAScavengerMaterial::PrintInfo()
141 { 162 {
142 auto pConfinedBox = fpChemistryInfo->GetChem 163 auto pConfinedBox = fpChemistryInfo->GetChemistryBoundary();
143 auto iter = fpChemistryInfo->begin(); << 164 auto iter = fpChemistryInfo->begin();
144 G4cout << "********************************* << 165 G4cout << "**************************************************************"
145 for (; iter != fpChemistryInfo->end(); iter+ << 166 << G4endl;
>> 167 for(; iter != fpChemistryInfo->end(); iter++)
>> 168 {
146 auto containedConf = iter->first; 169 auto containedConf = iter->first;
147 // auto concentration = iter->second; 170 // auto concentration = iter->second;
148 auto concentration = fScavengerTable[conta << 171 auto concentration =
>> 172 fScavengerTable[containedConf] / (Avogadro * pConfinedBox->Volume());
149 G4cout << "Scavenger:" << containedConf->G 173 G4cout << "Scavenger:" << containedConf->GetName() << " : "
150 << concentration / 1.0e-6 /*mm3 to 174 << concentration / 1.0e-6 /*mm3 to L*/ << " (M) with : "
151 << fScavengerTable[containedConf] < 175 << fScavengerTable[containedConf] << " (molecules)"
152 << "in: " << pConfinedBox->Volume() << 176 << "in: " << pConfinedBox->Volume() / (um * um * um) << " (um3)"
153 if (fScavengerTable[containedConf] < 1) { << 177 << G4endl;
>> 178 if(fScavengerTable[containedConf] < 1)
>> 179 {
154 G4cout << "!!!!!!!!!!!!! this molecule h 180 G4cout << "!!!!!!!!!!!!! this molecule has less one molecule for "
155 "considered volume" 181 "considered volume"
156 << G4endl; 182 << G4endl;
157 // assert(false); 183 // assert(false);
158 } 184 }
159 if (fVerbose != 0) { << 185 if(fVerbose != 0)
>> 186 {
160 Dump(); 187 Dump();
161 } 188 }
162 } 189 }
163 G4cout << "********************************* << 190 G4cout << "**************************************************************"
>> 191 << G4endl;
164 } 192 }
165 193
166 void G4DNAScavengerMaterial::Reset() 194 void G4DNAScavengerMaterial::Reset()
167 { 195 {
168 if (fpChemistryInfo == nullptr) { << 196 if(fpChemistryInfo == nullptr)
>> 197 {
169 return; 198 return;
170 } 199 }
171 200
172 if (fpChemistryInfo->size() == 0) { << 201 if(fpChemistryInfo->size() == 0)
>> 202 {
173 return; 203 return;
174 } 204 }
175 205
176 fScavengerTable.clear(); 206 fScavengerTable.clear();
177 fCounterMap.clear(); 207 fCounterMap.clear();
178 fpLastSearch.reset(nullptr); 208 fpLastSearch.reset(nullptr);
179 209
180 auto pConfinedBox = fpChemistryInfo->GetChem 210 auto pConfinedBox = fpChemistryInfo->GetChemistryBoundary();
181 auto iter = fpChemistryInfo->begin(); << 211 auto iter = fpChemistryInfo->begin();
182 for (; iter != fpChemistryInfo->end(); iter+ << 212 for(; iter != fpChemistryInfo->end(); iter++)
>> 213 {
183 auto containedConf = iter->first; 214 auto containedConf = iter->first;
184 auto concentration = iter->second; 215 auto concentration = iter->second;
185 fScavengerTable[containedConf] = floor(Avo << 216 fScavengerTable[containedConf] =
>> 217 floor(Avogadro * concentration * pConfinedBox->Volume());
186 fCounterMap[containedConf][1 * picosecond] 218 fCounterMap[containedConf][1 * picosecond] =
187 floor(Avogadro * concentration * pConfin 219 floor(Avogadro * concentration * pConfinedBox->Volume());
188 } 220 }
189 if (fVerbose != 0) { << 221 PrintInfo();
190 PrintInfo(); <<
191 } <<
192 } 222 }
193 223
194 //-------------------------------------------- 224 //------------------------------------------------------------------------------
195 225
196 void G4DNAScavengerMaterial::AddAMoleculeAtTim << 226 void G4DNAScavengerMaterial::AddAMoleculeAtTime(
197 << 227 MolType molecule, G4double time, const G4ThreeVector* /*position*/,
>> 228 int number)
198 { 229 {
199 if (fVerbose != 0) { << 230 if(fVerbose != 0)
200 G4cout << "G4DNAScavengerMaterial::AddAMol << 231 {
201 << " at time : " << G4BestUnit(time << 232 G4cout << "G4DNAScavengerMaterial::AddAMoleculeAtTime : "
>> 233 << molecule->GetName() << " at time : " << G4BestUnit(time, "Time")
>> 234 << G4endl;
202 } 235 }
203 236
204 auto counterMap_i = fCounterMap.find(molecul 237 auto counterMap_i = fCounterMap.find(molecule);
205 238
206 if (counterMap_i == fCounterMap.end()) { << 239 if(counterMap_i == fCounterMap.end())
>> 240 {
207 fCounterMap[molecule][time] = number; 241 fCounterMap[molecule][time] = number;
208 } 242 }
209 else if (counterMap_i->second.empty()) { << 243 else if(counterMap_i->second.empty())
>> 244 {
210 counterMap_i->second[time] = number; 245 counterMap_i->second[time] = number;
211 } 246 }
212 else { << 247 else
>> 248 {
213 auto end = counterMap_i->second.rbegin(); 249 auto end = counterMap_i->second.rbegin();
214 250
215 if (end->first <= time << 251 if(end->first <= time || fabs(end->first - time) <=
216 || fabs(end->first - time) <= G4::Mole << 252 G4::MoleculeCounter::TimePrecision::fPrecision)
217 G4double newValue = end->second + number << 253 {
>> 254 G4double newValue = end->second + number;
218 counterMap_i->second[time] = newValue; 255 counterMap_i->second[time] = newValue;
219 if (newValue != (floor)(fScavengerTable[ << 256 if(newValue != (floor)(fScavengerTable[molecule])) // protection
220 { 257 {
221 G4String errMsg = "You are trying to a << 258 assert(false);
222 G4Exception("AddAMoleculeAtTime", "", <<
223 } 259 }
>> 260 // G4cout<<" AddAMoleculeAtTime : "<<molecule->GetName()<< " :
>> 261 // "<<newValue<<G4endl;
224 } 262 }
>> 263 // else
>> 264 // {
>> 265 //
>> 266 // G4ExceptionDescription errMsg;
>> 267 // errMsg << "Time of species "
>> 268 // << molecule->GetName() << " is "
>> 269 // << G4BestUnit(time, "Time") << " while "
>> 270 // << " global time is "
>> 271 // <<" end->first : "<<end->first
>> 272 // //<<
>> 273 // G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(),
>> 274 // "Time")
>> 275 // << G4endl;
>> 276 // G4Exception("G4DNAScavengerMaterial::AddAMoleculeAtTime",
>> 277 // "TIME_DONT_MATCH",
>> 278 // FatalException, errMsg);
>> 279 // }
225 } 280 }
226 } 281 }
227 282
228 //-------------------------------------------- 283 //------------------------------------------------------------------------------
229 284
230 void G4DNAScavengerMaterial::RemoveAMoleculeAt << 285 void G4DNAScavengerMaterial::RemoveAMoleculeAtTime(
231 << 286 MolType pMolecule, G4double time, const G4ThreeVector* /*position*/,
>> 287 int number)
232 { 288 {
233 NbMoleculeInTime& nbMolPerTime = fCounterMap << 289 NbMoleculeAgainstTime& nbMolPerTime = fCounterMap[pMolecule];
234 290
235 if (fVerbose != 0) { << 291 if(fVerbose != 0)
>> 292 {
236 auto it_ = nbMolPerTime.rbegin(); 293 auto it_ = nbMolPerTime.rbegin();
237 G4cout << "G4DNAScavengerMaterial::RemoveA << 294 G4cout << "G4DNAScavengerMaterial::RemoveAMoleculeAtTime : "
238 << " at time : " << G4BestUnit(time << 295 << pMolecule->GetName() << " at time : " << G4BestUnit(time, "Time")
239 296
240 << " form : " << it_->second << G4e 297 << " form : " << it_->second << G4endl;
241 } 298 }
242 299
243 if (nbMolPerTime.empty()) { << 300 if(nbMolPerTime.empty())
>> 301 {
244 Dump(); 302 Dump();
245 G4String errMsg = "You are trying to remov 303 G4String errMsg = "You are trying to remove molecule " +
246 pMolecule->GetName() + 304 pMolecule->GetName() +
247 " from the counter while 305 " from the counter while this kind of molecules has not "
248 "been registered yet"; 306 "been registered yet";
249 G4Exception("G4DNAScavengerMaterial::Remov << 307 G4Exception("G4DNAScavengerMaterial::RemoveAMoleculeAtTime", "",
>> 308 FatalErrorInArgument, errMsg);
250 309
251 return; 310 return;
252 } 311 }
>> 312 else
>> 313 {
>> 314 auto it = nbMolPerTime.rbegin();
253 315
254 auto it = nbMolPerTime.rbegin(); << 316 if(it == nbMolPerTime.rend())
255 << 317 {
256 if (it == nbMolPerTime.rend()) { << 318 it--;
257 it--; << 319
258 << 320 G4String errMsg = "There was no " + pMolecule->GetName() +
259 G4String errMsg = "There was no " + pMolec << 321 " recorded at the time or even before the time asked";
260 + " recorded at the time << 322 G4Exception("G4DNAScavengerMaterial::RemoveAMoleculeAtTime", "",
261 G4Exception("G4DNAScavengerMaterial::Remov << 323 FatalErrorInArgument, errMsg);
262 } << 324 }
263 325
264 G4double finalN = it->second - number; << 326 G4double finalN = it->second - number;
265 if (finalN < 0) { << 327 if(finalN < 0)
266 Dump(); << 328 {
>> 329 Dump();
267 330
268 G4cout << "fScavengerTable : " << pMolecul << 331 G4cout << "fScavengerTable : " << pMolecule->GetName() << " : "
269 << G4endl; << 332 << (fScavengerTable[pMolecule]) << G4endl;
270 333
271 G4ExceptionDescription errMsg; << 334 G4ExceptionDescription errMsg;
272 errMsg << "After removal of " << number << << 335 errMsg << "After removal of " << number << " species of "
273 << " " << it->second << " " << pMol << 336 << " " << it->second << " " << pMolecule->GetName()
274 << G4BestUnit(time, "Time") << " is << 337 << " the final number at time " << G4BestUnit(time, "Time")
275 G4cout << " Global time is " << G4BestUnit << 338 << " is less than zero and so not valid." << G4endl;
276 << ". Previous selected time is " < << 339 G4cout << " Global time is "
277 G4Exception("G4DNAScavengerMaterial::Remov << 340 << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
278 } << 341 << ". Previous selected time is " << G4BestUnit(it->first, "Time")
279 nbMolPerTime[time] = finalN; << 342 << G4endl;
>> 343 G4Exception("G4DNAScavengerMaterial::RemoveAMoleculeAtTime", "N_INF_0",
>> 344 FatalException, errMsg);
>> 345 }
>> 346 nbMolPerTime[time] = finalN;
280 347
281 if (finalN != (floor)(fScavengerTable[pMolec << 348 if(finalN != (floor)(fScavengerTable[pMolecule])) // protection
282 { << 349 {
283 assert(false); << 350 assert(false);
>> 351 }
284 } 352 }
285 } 353 }
286 354
287 void G4DNAScavengerMaterial::Dump() 355 void G4DNAScavengerMaterial::Dump()
288 { 356 {
289 auto pConfinedBox = fpChemistryInfo->GetChem 357 auto pConfinedBox = fpChemistryInfo->GetChemistryBoundary();
290 auto V = pConfinedBox->Volume(); << 358 auto V = pConfinedBox->Volume();
291 for (const auto& it : fCounterMap) { << 359 for(auto it : fCounterMap)
>> 360 {
292 auto pReactant = it.first; 361 auto pReactant = it.first;
293 362
294 G4cout << " --- > For " << pReactant->GetN 363 G4cout << " --- > For " << pReactant->GetName() << G4endl;
295 364
296 for (const auto& it2 : it.second) { << 365 for(auto it2 : it.second)
>> 366 {
297 G4cout << " " << G4BestUnit(it2.first, " 367 G4cout << " " << G4BestUnit(it2.first, "Time") << " "
298 << it2.second / (Avogadro * V * 1 368 << it2.second / (Avogadro * V * 1.0e-6 /*mm3 to L*/) << G4endl;
299 } 369 }
300 } 370 }
301 } 371 }
302 372
303 int64_t G4DNAScavengerMaterial::GetNMoleculesA << 373 int G4DNAScavengerMaterial::GetNMoleculesAtTime(MolType molecule, G4double time)
304 { 374 {
305 if (!fCounterAgainstTime) { << 375 if(!fCounterAgainstTime)
>> 376 {
306 G4cout << "fCounterAgainstTime == false" < 377 G4cout << "fCounterAgainstTime == false" << G4endl;
307 assert(false); 378 assert(false);
308 } 379 }
309 380
310 G4bool sameTypeOfMolecule = SearchTimeMap(mo 381 G4bool sameTypeOfMolecule = SearchTimeMap(molecule);
311 auto output = SearchUpperBoundTime(time, sam << 382 return SearchUpperBoundTime(time, sameTypeOfMolecule);
312 if (output < 0) { <<
313 G4ExceptionDescription errMsg; <<
314 errMsg << "N molecules not valid < 0 : " < <<
315 G4Exception("G4DNAScavengerMaterial::GetNM <<
316 } <<
317 return output; <<
318 } 383 }
319 384
320 G4bool G4DNAScavengerMaterial::SearchTimeMap(M 385 G4bool G4DNAScavengerMaterial::SearchTimeMap(MolType molecule)
321 { 386 {
322 if (fpLastSearch == nullptr) { << 387 if(fpLastSearch == nullptr)
>> 388 {
323 fpLastSearch = std::make_unique<Search>(); 389 fpLastSearch = std::make_unique<Search>();
324 } 390 }
325 else { << 391 else
326 if (fpLastSearch->fLowerBoundSet && fpLast << 392 {
>> 393 if(fpLastSearch->fLowerBoundSet &&
>> 394 fpLastSearch->fLastMoleculeSearched->first == molecule)
>> 395 {
327 return true; 396 return true;
328 } 397 }
329 } 398 }
330 399
331 auto mol_it = fCounterMap.find(molecule); << 400 auto mol_it = fCounterMap.find(molecule);
332 fpLastSearch->fLastMoleculeSearched = mol_it 401 fpLastSearch->fLastMoleculeSearched = mol_it;
333 402
334 if (mol_it != fCounterMap.end()) { << 403 if(mol_it != fCounterMap.end())
335 fpLastSearch->fLowerBoundTime = fpLastSear << 404 {
>> 405 fpLastSearch->fLowerBoundTime =
>> 406 fpLastSearch->fLastMoleculeSearched->second.end();
336 fpLastSearch->fLowerBoundSet = true; 407 fpLastSearch->fLowerBoundSet = true;
337 } 408 }
338 else { << 409 else
>> 410 {
339 fpLastSearch->fLowerBoundSet = false; 411 fpLastSearch->fLowerBoundSet = false;
340 } 412 }
341 413
342 return false; 414 return false;
343 } 415 }
344 416
345 //-------------------------------------------- 417 //------------------------------------------------------------------------------
346 418
347 int64_t G4DNAScavengerMaterial::SearchUpperBou << 419 int G4DNAScavengerMaterial::SearchUpperBoundTime(G4double time,
>> 420 G4bool sameTypeOfMolecule)
348 { 421 {
349 auto mol_it = fpLastSearch->fLastMoleculeSea 422 auto mol_it = fpLastSearch->fLastMoleculeSearched;
350 if (mol_it == fCounterMap.end()) { << 423 if(mol_it == fCounterMap.end())
>> 424 {
351 return 0; 425 return 0;
352 } 426 }
353 427
354 NbMoleculeInTime& timeMap = mol_it->second; << 428 NbMoleculeAgainstTime& timeMap = mol_it->second;
355 if (timeMap.empty()) { << 429 if(timeMap.empty())
>> 430 {
356 return 0; 431 return 0;
357 } 432 }
358 433
359 if (sameTypeOfMolecule) { << 434 if(sameTypeOfMolecule)
360 if (fpLastSearch->fLowerBoundSet && fpLast << 435 {
361 if (fpLastSearch->fLowerBoundTime->first << 436 if(fpLastSearch->fLowerBoundSet &&
>> 437 fpLastSearch->fLowerBoundTime != timeMap.end())
>> 438 {
>> 439 if(fpLastSearch->fLowerBoundTime->first < time)
>> 440 {
362 auto upperToLast = fpLastSearch->fLowe 441 auto upperToLast = fpLastSearch->fLowerBoundTime;
363 upperToLast++; 442 upperToLast++;
364 443
365 if (upperToLast == timeMap.end()) { << 444 if(upperToLast == timeMap.end())
>> 445 {
366 return fpLastSearch->fLowerBoundTime 446 return fpLastSearch->fLowerBoundTime->second;
367 } 447 }
368 448
369 if (upperToLast->first > time) { << 449 if(upperToLast->first > time)
>> 450 {
370 return fpLastSearch->fLowerBoundTime 451 return fpLastSearch->fLowerBoundTime->second;
371 } 452 }
372 } 453 }
373 } 454 }
374 } 455 }
375 456
376 auto up_time_it = timeMap.upper_bound(time); 457 auto up_time_it = timeMap.upper_bound(time);
377 458
378 if (up_time_it == timeMap.end()) { << 459 if(up_time_it == timeMap.end())
>> 460 {
379 auto last_time = timeMap.rbegin(); 461 auto last_time = timeMap.rbegin();
380 return last_time->second; 462 return last_time->second;
381 } 463 }
382 if (up_time_it == timeMap.begin()) { << 464 if(up_time_it == timeMap.begin())
>> 465 {
383 return 0; 466 return 0;
384 } 467 }
385 468
386 up_time_it--; 469 up_time_it--;
387 470
388 fpLastSearch->fLowerBoundTime = up_time_it; 471 fpLastSearch->fLowerBoundTime = up_time_it;
389 fpLastSearch->fLowerBoundSet = true; << 472 fpLastSearch->fLowerBoundSet = true;
390 473
391 return fpLastSearch->fLowerBoundTime->second 474 return fpLastSearch->fLowerBoundTime->second;
392 } <<
393 <<
394 void G4DNAScavengerMaterial::WaterEquilibrium( <<
395 { <<
396 auto convertFactor = Avogadro * fpChemistryI <<
397 G4double kw = 1.01e-14; <<
398 fScavengerTable[fHOm] = (kw / ((G4double)fSc <<
399 G4cout << "pH : " << GetpH() << G4endl; <<
400 return; <<
401 } <<
402 <<
403 void G4DNAScavengerMaterial::SetpH(const G4int <<
404 { <<
405 auto volume = fpChemistryInfo->GetChemistryB <<
406 fScavengerTable[fH3Op] = floor(Avogadro * st <<
407 fScavengerTable[fHOm] = floor(Avogadro * std <<
408 } <<
409 <<
410 G4double G4DNAScavengerMaterial::GetpH() <<
411 { <<
412 G4double volumeInLiter = fpChemistryInfo->Ge <<
413 G4double Cion = (G4double)fScavengerTable[fH <<
414 G4double pH = std::log10(Cion); <<
415 // G4cout<<"OH- : "<<fScavengerTable[fHOm]<< <<
416 // "<<-pH<<G4endl; <<
417 if (fScavengerTable[fH3Op] < 0) // protect <<
418 { <<
419 G4Exception("G4DNAScavengerMaterial::GetpH <<
420 "H3O+ < 0"); <<
421 fScavengerTable[fH3Op] = 0; <<
422 } <<
423 if (fScavengerTable[fHOm] < 0) // protect m <<
424 { <<
425 G4Exception("G4DNAScavengerMaterial::GetpH <<
426 "HO- < 0"); <<
427 fScavengerTable[fHOm] = 0; <<
428 } <<
429 return -pH; <<
430 } 475 }