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