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