Geant4 11.3.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNAWaterDissociationDisplacer.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// Author: Mathieu Karamitros
28//
29// WARNING : This class is released as a prototype.
30// It might strongly evolve or even disapear in the next releases.
31//
32// History:
33// -----------
34// 10 Oct 2011 M.Karamitros created
35//
36// -------------------------------------------------------------------
37
38#include "G4DNAWaterDissociationDisplacer.hh"
39#include "G4PhysicalConstants.hh"
40#include "G4SystemOfUnits.hh"
41#include "G4H2O.hh"
42#include "G4H2.hh"
43#include "G4Hydrogen.hh"
44#include "G4Oxygen.hh"
45#include "G4OH.hh"
46#include "G4H3O.hh"
47#include "G4Electron_aq.hh"
48#include "G4H2O2.hh"
49#include "Randomize.hh"
50#include "G4Molecule.hh"
51#include "G4MolecularConfiguration.hh"
52#include "G4RandomDirection.hh"
53#include "G4Exp.hh"
54#include "G4UnitsTable.hh"
55#include "G4EmParameters.hh"
56#include "G4DNAOneStepThermalizationModel.hh"
57
58using namespace std;
59
60//------------------------------------------------------------------------------
61
62G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
63 Ionisation_DissociationDecay)
64
65G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
66 A1B1_DissociationDecay)
67
68G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
69 B1A1_DissociationDecay)
70
71G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
72 B1A1_DissociationDecay2)
73
74G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
75 AutoIonisation)
76
77G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
78 DissociativeAttachment)
79
80G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
81 DoubleIonisation_DissociationDecay1)
82
83G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
84 DoubleIonisation_DissociationDecay2)
85
86G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
87 DoubleIonisation_DissociationDecay3)
88
89G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
90 TripleIonisation_DissociationDecay)
91
92G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
93 QuadrupleIonisation_DissociationDecay)
94
95/*
96//------------------------------------------------------------------------------
97#ifdef _WATER_DISPLACER_USE_KREIPL_
98// This function is used to generate the following density distribution:
99// f(r) := 4*r * exp(-2*r)
100// reference:
101// Kreipl, M. S., Friedland, W. & Paretzke, H. G.
102// Time- and space-resolved Monte Carlo study of water radiolysis for photon,
103// electron and ion irradiation.
104// Radiat. Environ. Biophys. 48, 11–20 (2009).
105G4double G4DNAWaterDissociationDisplacer::ElectronProbaDistribution(G4double r)
106{
107 return (2.*r+1.)*G4Exp(-2.*r);
108}
109#endif
110
111//------------------------------------------------------------------------------
112#ifdef _WATER_DISPLACER_USE_TERRISOL_
113// This function is used to generate the following density distribution:
114// f(r) := 4*x^2/(sqrt(%pi)*(b)^3)*exp(-x^2/(b)^2)
115// with b=27.22 for 7 eV
116// reference
117// Terrissol M, Beaudre A (1990) Simulation of space and time evolution
118// of radiolytic species induced by electrons in water.
119// Radiat Prot Dosimetry 31:171–175
120
121G4double G4DNAWaterDissociationDisplacer::ElectronProbaDistribution(G4double r)
122{
123#define b 27.22 // nanometer
124 static constexpr double sqrt_pi = 1.77245; // sqrt(CLHEP::pi);
125 static constexpr double b_to3 = 20168.1; // pow(b,3.);
126 static constexpr double b_to2 = 740.928; // pow(b,2.);
127 static constexpr double inverse_b_to2 = 1. / b_to2;
128
129 static constexpr double main_factor = 4. / (sqrt_pi * b_to3);
130 static constexpr double factorA = sqrt_pi * b_to3 / 4.;
131 static constexpr double factorB = b_to2 / 2.;
132
133 return (main_factor *
134 (factorA * erf(r / b)
135 - factorB * r * G4Exp(-pow(r, 2.) * inverse_b_to2)));
136}
137
138#endif
139//------------------------------------------------------------------------------
140*/
141G4DNAWaterDissociationDisplacer::G4DNAWaterDissociationDisplacer()
142 :
143
144 ke(1.7*eV)
145/*#ifdef _WATER_DISPLACER_USE_KREIPL_
146 fFastElectronDistrib(0., 5., 0.001)
147#elif defined _WATER_DISPLACER_USE_TERRISOL_
148 fFastElectronDistrib(0., 100., 0.001)
149#endif*/
150{/*
151 fProba1DFunction =
152 std::bind((G4double(*)(G4double))
153 &G4DNAWaterDissociationDisplacer::ElectronProbaDistribution,
154 std::placeholders::_1);
155
156 size_t nBins = 500;
157 fElectronThermalization.reserve(nBins);
158 double eps = 1. / ((int) nBins);
159 double proba = eps;
160
161 fElectronThermalization.push_back(0.);
162
163 for (size_t i = 1; i < nBins; ++i)
164 {
165 double r = fFastElectronDistrib.Revert(proba, fProba1DFunction);
166 fElectronThermalization.push_back(r * nanometer);
167 proba += eps;
168// G4cout << G4BestUnit(r*nanometer, "Length") << G4endl;
169 }*/
170 dnaSubType = G4EmParameters::Instance()->DNAeSolvationSubType();
171// SetVerbose(1);
172}
173
174//------------------------------------------------------------------------------
175
179
180//------------------------------------------------------------------------------
181
182G4ThreeVector G4DNAWaterDissociationDisplacer::GetMotherMoleculeDisplacement(
183 const G4MolecularDissociationChannel* theDecayChannel) const
184{
185 G4int decayType = theDecayChannel->GetDisplacementType();
186 G4double RMSMotherMoleculeDisplacement(0.0);
187
188 switch (decayType) {
189 case Ionisation_DissociationDecay:
190 RMSMotherMoleculeDisplacement = 2.0 * nanometer;
191 break;
192 case A1B1_DissociationDecay:
193 RMSMotherMoleculeDisplacement = 0.0 * nanometer;
194 break;
195 case B1A1_DissociationDecay:
196 RMSMotherMoleculeDisplacement = 0.0 * nanometer;
197 break;
198 case B1A1_DissociationDecay2:
199 RMSMotherMoleculeDisplacement = 0.0 * nanometer;
200 break;
201 case AutoIonisation:
202 RMSMotherMoleculeDisplacement = 2.0 * nanometer;
203 break;
204 case DissociativeAttachment:
205 RMSMotherMoleculeDisplacement = 0.0 * nanometer;
206 break;
207 case DoubleIonisation_DissociationDecay1:
208 RMSMotherMoleculeDisplacement = 2.0 * nanometer;
209 break;
210 case DoubleIonisation_DissociationDecay2:
211 RMSMotherMoleculeDisplacement = 2.0 * nanometer;
212 break;
213 case DoubleIonisation_DissociationDecay3:
214 RMSMotherMoleculeDisplacement = 2.0 * nanometer;
215 break;
216 case TripleIonisation_DissociationDecay:
217 RMSMotherMoleculeDisplacement = 2.0 * nanometer;
218 break;
219 case QuadrupleIonisation_DissociationDecay:
220 RMSMotherMoleculeDisplacement = 2.0 * nanometer;
221 break;
222 }
223
224 if (RMSMotherMoleculeDisplacement == 0) {
225 return G4ThreeVector(0, 0, 0);
226 }
227
228 auto RandDirection =
229 radialDistributionOfProducts(RMSMotherMoleculeDisplacement);
230
231 return RandDirection;
232}
233
234//------------------------------------------------------------------------------
235
236vector<G4ThreeVector> G4DNAWaterDissociationDisplacer::GetProductsDisplacement(
237 const G4MolecularDissociationChannel* pDecayChannel) const
238{
239 G4int nbProducts = pDecayChannel->GetNbProducts();
240 vector<G4ThreeVector> theProductDisplacementVector(nbProducts);
241
242 typedef map<const G4MoleculeDefinition*, G4double> RMSmap;
243 RMSmap theRMSmap;
244
245 G4int decayType = pDecayChannel->GetDisplacementType();
246
247 switch (decayType) {
248 case Ionisation_DissociationDecay:
249 {
250 if (fVerbose != 0) {
251 G4cout << "Ionisation_DissociationDecay" << G4endl;
252 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
253 }
254 G4double RdmValue = G4UniformRand();
255
256 if (RdmValue < 0.5) {
257 // H3O
258 theRMSmap[G4H3O::Definition()] = 0.0 * nanometer;
259 // OH
260 theRMSmap[G4OH::Definition()] = 0.8 * nanometer;
261 } else {
262 // H3O
263 theRMSmap[G4H3O::Definition()] = 0.8 * nanometer;
264 // OH
265 theRMSmap[G4OH::Definition()] = 0.0 * nanometer;
266 }
267
268 for (G4int i = 0; i < nbProducts; i++) {
269 auto pProduct = pDecayChannel->GetProduct(i);
270 G4double theRMSDisplacement = theRMSmap[pProduct->GetDefinition()];
271
272 if (theRMSDisplacement == 0.0) {
273 theProductDisplacementVector[i] = G4ThreeVector();
274 } else {
275 auto RandDirection =
276 radialDistributionOfProducts(theRMSDisplacement);
277 theProductDisplacementVector[i] = RandDirection;
278 }
279 }
280 break;
281 }
282 case A1B1_DissociationDecay:
283 {
284 if (fVerbose != 0) {
285 G4cout << "A1B1_DissociationDecay" << G4endl;
286 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
287 }
288
289 constexpr G4double theRMSDisplacement = 2.4 * nanometer;
290 auto RandDirection = radialDistributionOfProducts(theRMSDisplacement);
291
292 for (G4int i = 0; i < nbProducts; i++) {
293 auto pProduct = pDecayChannel->GetProduct(i);
294
295 if (pProduct->GetDefinition() == G4OH::Definition()) {
296 theProductDisplacementVector[i] = -1.0 / 18.0 * RandDirection;
297 } else if (pProduct->GetDefinition() == G4Hydrogen::Definition()) {
298 theProductDisplacementVector[i] = +17.0 / 18.0 * RandDirection;
299 }
300 }
301 break;
302 }
303 case B1A1_DissociationDecay:
304 {
305 if (fVerbose != 0) {
306 G4cout << "B1A1_DissociationDecay" << G4endl;
307 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
308 }
309
310 constexpr G4double theRMSDisplacement = 0.8 * nanometer;
311 auto RandDirection = radialDistributionOfProducts(theRMSDisplacement);
312
313 G4int NbOfOH = 0;
314 for (G4int i = 0; i < nbProducts; ++i) {
315 auto pProduct = pDecayChannel->GetProduct(i);
316 if (pProduct->GetDefinition() == G4H2::Definition()) {
317 // In the paper of Kreipl (2009)
318 // theProductDisplacementVector[i] = -2.0 / 18.0 * RandDirection;
319
320 // Based on momentum conservation
321 theProductDisplacementVector[i] = -16.0 / 18.0 * RandDirection;
322 } else if (pProduct->GetDefinition() == G4OH::Definition()) {
323 // In the paper of Kreipl (2009)
324 // G4ThreeVector OxygenDisplacement = +16.0 / 18.0 * RandDirection;
325
326 // Based on momentum conservation
327 G4ThreeVector OxygenDisplacement = +2.0 / 18.0 * RandDirection;
328 constexpr G4double OHRMSDisplacement = 1.1 * nanometer;
329
330 auto OHDisplacement =
331 radialDistributionOfProducts(OHRMSDisplacement);
332
333 if (NbOfOH == 0) {
334 OHDisplacement = 0.5 * OHDisplacement;
335 } else {
336 OHDisplacement = -0.5 * OHDisplacement;
337 }
338
339 theProductDisplacementVector[i] =
340 OHDisplacement + OxygenDisplacement;
341
342 ++NbOfOH;
343 }
344 }
345 break;
346 }
347 case B1A1_DissociationDecay2:
348 {
349 if (fVerbose != 0){
350 G4cout << "B1A1_DissociationDecay2" << G4endl;
351 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
352 }
353
354 G4int NbOfH = 0;
355 for (G4int i = 0; i < nbProducts; ++i) {
356 auto pProduct = pDecayChannel->GetProduct(i);
357 if (pProduct->GetDefinition() == G4Oxygen::Definition()) {
358 // O(3p)
359 theProductDisplacementVector[i] = G4ThreeVector(0, 0, 0);
360 } else if (pProduct->GetDefinition() == G4Hydrogen::Definition()) {
361 // H
362 constexpr G4double HRMSDisplacement = 1.6 * nanometer;
363
364 auto HDisplacement =
365 radialDistributionOfProducts(HRMSDisplacement);
366
367 if (NbOfH == 0) {
368 HDisplacement = 0.5 * HDisplacement;
369 } else {
370 HDisplacement = -0.5 * HDisplacement;
371 }
372 theProductDisplacementVector[i] = HDisplacement;
373
374 ++NbOfH;
375 }
376 }
377 break;
378 }
379 case AutoIonisation:
380 {
381 if (fVerbose != 0) {
382 G4cout << "AutoIonisation" << G4endl;
383 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
384 }
385
386 G4double RdmValue = G4UniformRand();
387
388 if (RdmValue < 0.5) {
389 // H3O
390 theRMSmap[G4H3O::Definition()] = 0.0 * nanometer;
391 // OH
392 theRMSmap[G4OH::Definition()] = 0.8 * nanometer;
393 } else {
394 // H3O
395 theRMSmap[G4H3O::Definition()] = 0.8 * nanometer;
396 // OH
397 theRMSmap[G4OH::Definition()] = 0.0 * nanometer;
398 }
399
400 for (G4int i = 0; i < nbProducts; i++) {
401 auto pProduct = pDecayChannel->GetProduct(i);
402 auto theRMSDisplacement = theRMSmap[pProduct->GetDefinition()];
403
404 if (theRMSDisplacement == 0) {
405 theProductDisplacementVector[i] = G4ThreeVector();
406 } else {
407 auto RandDirection =
408 radialDistributionOfProducts(theRMSDisplacement);
409 theProductDisplacementVector[i] = RandDirection;
410 }
411 if (pProduct->GetDefinition() == G4Electron_aq::Definition()) {
412 theProductDisplacementVector[i] = radialDistributionOfElectron();
413 }
414 }
415 break;
416 }
417 case DissociativeAttachment:
418 {
419 if (fVerbose != 0) {
420 G4cout << "DissociativeAttachment" << G4endl;
421 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
422 }
423 constexpr G4double theRMSDisplacement = 0.8 * nanometer;
424 auto RandDirection = radialDistributionOfProducts(theRMSDisplacement);
425
426 G4int NbOfOH = 0;
427 for (G4int i = 0; i < nbProducts; ++i) {
428 auto pProduct = pDecayChannel->GetProduct(i);
429 if (pProduct->GetDefinition() == G4H2::Definition()) {
430 // In the paper of Kreipl (2009)
431 // theProductDisplacementVector[i] = -2.0 / 18.0 * RandDirection;
432
433 // Based on momentum conservation
434 theProductDisplacementVector[i] = -16.0 / 18.0 * RandDirection;
435 } else if (pProduct->GetDefinition() == G4OH::Definition()) {
436 // In the paper of Kreipl (2009)
437 // G4ThreeVector OxygenDisplacement = +16.0 / 18.0 * RandDirection;
438
439 // Based on momentum conservation
440 G4ThreeVector OxygenDisplacement = +2.0 / 18.0 * RandDirection;
441 constexpr G4double OHRMSDisplacement = 1.1 * nanometer;
442
443 auto OHDisplacement =
444 radialDistributionOfProducts(OHRMSDisplacement);
445
446 if (NbOfOH == 0) {
447 OHDisplacement = 0.5 * OHDisplacement;
448 } else {
449 OHDisplacement = -0.5 * OHDisplacement;
450 }
451 theProductDisplacementVector[i] = OHDisplacement +
452 OxygenDisplacement;
453 ++NbOfOH;
454 }
455 }
456 break;
457 }
458 case DoubleIonisation_DissociationDecay1:
459 {
460 if (fVerbose != 0) {
461 G4cout << "DoubleIonisation_DissociationDecay1" << G4endl;
462 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
463 }
464
465 // Ref.) B. Gervais, et al., DOI: 10.1016/j.radphyschem.2005.09.01
466 // Decay Channel #1: H2O^2+ -> 2H+ + O(3P) -> 2H3O+ + O(3P)
467
468 theRMSmap[G4H3O::Definition()] = 1.2 * nanometer;
469 theRMSmap[G4Oxygen::Definition()] = 0.0 * nanometer;
470
471 for (G4int i = 0, num_H3O = 0; i < nbProducts; i++) {
472
473 const auto prod_def = pDecayChannel->GetProduct(i)->GetDefinition();
474
475 if (prod_def == G4H3O::Definition()) {
476 num_H3O++;
477 if (num_H3O == 2) {
478 constexpr G4double H3Op_rms = 0.3 * nanometer;
479 theRMSmap[G4H3O::Definition()] = H3Op_rms;
480 }
481 }
482
483 theProductDisplacementVector[i]
484 = radialDistributionOfProducts(theRMSmap[prod_def]);
485
486 }
487
488 break;
489 }
490 case DoubleIonisation_DissociationDecay2:
491 {
492 if (fVerbose != 0) {
493 G4cout << "DoubleIonisation_DissociationDecay2" << G4endl;
494 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
495 }
496
497 // Ref.) B. Gervais, et al., DOI: 10.1016/j.radphyschem.2005.09.01
498 // Decay Channel #2: H2O^2+ -> H+ + H* + O+ -> 2H3O+ + H* + *OH + O(3P)
499
500 theRMSmap[G4H3O::Definition()] = 1.2 * nanometer;
501 theRMSmap[G4Oxygen::Definition()] = 0.0 * nanometer;
502 theRMSmap[G4Hydrogen::Definition()] = 0.8 * nanometer;
503 theRMSmap[G4OH::Definition()] = 0.3 * nanometer;
504
505 const auto OH_disp = radialDistributionOfProducts(
506 theRMSmap[G4OH::Definition()]);
507
508 for (G4int i = 0, num_H3O = 0; i < nbProducts; i++) {
509
510 const auto prod_def = pDecayChannel->GetProduct(i)->GetDefinition();
511
512 if (prod_def == G4H3O::Definition()) {
513 num_H3O++;
514 if (num_H3O == 2) {
515 constexpr G4double OH_rms = 0.3 * nanometer;
516 theProductDisplacementVector[i]
517 = radialDistributionOfProducts(OH_rms);
518 theProductDisplacementVector[i] += OH_disp;
519 continue;
520 }
521 } else if (prod_def == G4OH::Definition()) {
522 theProductDisplacementVector[i] = OH_disp;
523 continue;
524 }
525
526 theProductDisplacementVector[i]
527 = radialDistributionOfProducts(theRMSmap[prod_def]);
528
529 }
530
531 break;
532 }
533 case DoubleIonisation_DissociationDecay3:
534 {
535 if (fVerbose != 0) {
536 G4cout << "DoubleIonisation_DissociationDecay3" << G4endl;
537 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
538 }
539
540 // Ref.) B. Gervais, et al., DOI: 10.1016/j.radphyschem.2005.09.01
541 // Decay Channel #3: H2O^2+ -> H+ + OH+ -> 2H3O+ + O(3P)
542
543 theRMSmap[G4H3O::Definition()] = 1.2 * nanometer;
544 theRMSmap[G4Oxygen::Definition()] = 0.0 * nanometer;
545
546 for (G4int i = 0; i < nbProducts; i++) {
547
548 const auto prod_def = pDecayChannel->GetProduct(i)->GetDefinition();
549
550 theProductDisplacementVector[i]
551 = radialDistributionOfProducts(theRMSmap[prod_def]);
552
553 }
554
555 break;
556 }
557 case TripleIonisation_DissociationDecay:
558 {
559 if (fVerbose != 0) {
560 G4cout << "TripleIonisation_DissociationDecay" << G4endl;
561 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
562 }
563
564 // Ref.) B. Gervais, et al., DOI: 10.1016/j.radphyschem.2005.09.01
565 // Decay Channel: H2O^3+ -> 3H3O+ + *OH + O(3P)
566 //
567 // ** Detaied decay chain **
568 // H2O^3+ -> O+ + H3O+ + H3O+
569 // |
570 // --> H2O+ + O(3P)
571 // |
572 // --> *OH + H3O+
573
574 theRMSmap[G4H3O::Definition()] = 1.2 * nanometer;
575 theRMSmap[G4Oxygen::Definition()] = 0.0 * nanometer;
576 theRMSmap[G4OH::Definition()] = 0.3 * nanometer;
577
578 const auto OH_disp = radialDistributionOfProducts(
579 theRMSmap[G4OH::Definition()]);
580
581 for (G4int i = 0, num_H3O = 0; i < nbProducts; i++) {
582
583 const auto prod_def = pDecayChannel->GetProduct(i)->GetDefinition();
584
585 if (prod_def == G4H3O::Definition()) {
586 num_H3O++;
587 if (num_H3O == 3) {
588 constexpr G4double H3Op_rms = 0.3 * nanometer;
589 theProductDisplacementVector[i]
590 = radialDistributionOfProducts(H3Op_rms);
591 theProductDisplacementVector[i] += OH_disp;
592 continue;
593 }
594 } else if (prod_def == G4OH::Definition()) {
595 theProductDisplacementVector[i] = OH_disp;
596 continue;
597 }
598
599 theProductDisplacementVector[i]
600 = radialDistributionOfProducts(theRMSmap[prod_def]);
601
602 }
603
604 break;
605 }
606 case QuadrupleIonisation_DissociationDecay:
607 {
608 if (fVerbose != 0) {
609 G4cout << "QuadrupleIonisation_DissociationDecay" << G4endl;
610 G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
611 }
612
613 // Ref.) B. Gervais, et al., DOI: 10.1016/j.radphyschem.2005.09.01
614 // Decay Channel: H2O^4+ -> 4H3O+ + 2*OH + O(3P)
615 //
616 // ** Detaied decay chain **
617 // H2O^4+ -> O^2+ + H3O+ + H3O+
618 // |
619 // --> H2O+ + H2O+ + O(3P)
620 // | |
621 // | --> *OH + H3O+
622 // |
623 // --> *OH + H3O+
624
625 theRMSmap[G4H3O::Definition()] = 1.2 * nanometer;
626 theRMSmap[G4Oxygen::Definition()] = 0.0 * nanometer;
627 theRMSmap[G4OH::Definition()] = 0.3 * nanometer;
628
629 const auto OH_disp1 = radialDistributionOfProducts(
630 theRMSmap[G4OH::Definition()]);
631 const auto OH_disp2 = radialDistributionOfProducts(
632 theRMSmap[G4OH::Definition()]);
633
634 for (G4int i = 0, num_H3O = 0, num_OH = 0; i < nbProducts; i++) {
635
636 const auto prod_def = pDecayChannel->GetProduct(i)->GetDefinition();
637
638 if (prod_def == G4H3O::Definition()) {
639
640 num_H3O++;
641 constexpr G4double H3Op_rms = 0.3 * nanometer;
642
643 if (num_H3O == 3) {
644 theProductDisplacementVector[i]
645 = radialDistributionOfProducts(H3Op_rms);
646 theProductDisplacementVector[i] += OH_disp1;
647 continue;
648 } else if (num_H3O == 4) {
649 theProductDisplacementVector[i]
650 = radialDistributionOfProducts(H3Op_rms);
651 theProductDisplacementVector[i] += OH_disp2;
652 continue;
653 }
654
655 } else if (prod_def == G4OH::Definition()) {
656
657 num_OH++;
658 if (num_OH == 1) {
659 theProductDisplacementVector[i] = OH_disp1;
660 } else {
661 theProductDisplacementVector[i] = OH_disp2;
662 }
663
664 continue;
665 }
666
667 theProductDisplacementVector[i]
668 = radialDistributionOfProducts(theRMSmap[prod_def]);
669
670 }
671
672 break;
673 }
674
675 }
676 return theProductDisplacementVector;
677}
678
679//------------------------------------------------------------------------------
680
681G4ThreeVector
682G4DNAWaterDissociationDisplacer::
683radialDistributionOfProducts(G4double Rrms) const
684{
685 static const double inverse_sqrt_3 = 1.0 / sqrt(3.0);
686 double sigma = Rrms * inverse_sqrt_3;
687 double x = G4RandGauss::shoot(0.0, sigma);
688 double y = G4RandGauss::shoot(0.0, sigma);
689 double z = G4RandGauss::shoot(0.0, sigma);
690 return G4ThreeVector(x, y, z);
691}
692
693//------------------------------------------------------------------------------
694
695G4ThreeVector
696G4DNAWaterDissociationDisplacer::radialDistributionOfElectron() const
697{/*
698 G4double rand_value = G4UniformRand();
699 size_t nBins = fElectronThermalization.size();
700 size_t bin = size_t(floor(rand_value * nBins));
701 size_t bin_p1 = min(bin + 1, nBins - 1);
702
703 return (fElectronThermalization[bin] * (1. - rand_value)
704 + fElectronThermalization[bin_p1] * rand_value) *
705 G4RandomDirection();*/
706
707 G4ThreeVector pdf = G4ThreeVector(0, 0, 0);
708
709 if (dnaSubType == fRitchie1994eSolvation) {
710 DNA::Penetration::Ritchie1994::GetPenetration(ke, pdf);
711 } else if (dnaSubType == fTerrisol1990eSolvation) {
712 DNA::Penetration::Terrisol1990::GetPenetration(ke, pdf);
713 } else if (dnaSubType == fMeesungnoensolid2002eSolvation) {
714 DNA::Penetration::Meesungnoen2002_amorphous::GetPenetration(ke, pdf);
715 } else if (dnaSubType == fKreipl2009eSolvation) {
716 DNA::Penetration::Kreipl2009::GetPenetration(ke, pdf);
717 } else {
718 DNA::Penetration::Meesungnoen2002::GetPenetration(ke, pdf);
719 }
720 return pdf;
721}
G4CT_COUNT_IMPL
#define G4CT_COUNT_IMPL(enumName, flagName)
Definition G4CTCounter.hh:137
fKreipl2009eSolvation
@ fKreipl2009eSolvation
Definition G4DNAModelSubType.hh:48
fMeesungnoensolid2002eSolvation
@ fMeesungnoensolid2002eSolvation
Definition G4DNAModelSubType.hh:49
fRitchie1994eSolvation
@ fRitchie1994eSolvation
Definition G4DNAModelSubType.hh:45
fTerrisol1990eSolvation
@ fTerrisol1990eSolvation
Definition G4DNAModelSubType.hh:46
G4Electron_aq.hh
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
G4DNAWaterDissociationDisplacer::GetMotherMoleculeDisplacement
G4ThreeVector GetMotherMoleculeDisplacement(const G4MolecularDissociationChannel *) const override
Definition G4DNAWaterDissociationDisplacer.cc:182
G4DNAWaterDissociationDisplacer::GetProductsDisplacement
std::vector< G4ThreeVector > GetProductsDisplacement(const G4MolecularDissociationChannel *) const override
Definition G4DNAWaterDissociationDisplacer.cc:236
G4DNAWaterDissociationDisplacer::radialDistributionOfProducts
G4ThreeVector radialDistributionOfProducts(G4double r_rms) const
Definition G4DNAWaterDissociationDisplacer.cc:683
G4Electron_aq::Definition
static G4Electron_aq * Definition()
Definition G4Electron_aq.cc:45
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition G4EmParameters.cc:70
G4EmParameters::DNAeSolvationSubType
G4DNAModelSubType DNAeSolvationSubType() const
Definition G4EmParameters.cc:1142
G4H3O::Definition
static G4H3O * Definition()
Definition G4H3O.cc:46
G4Hydrogen::Definition
static G4Hydrogen * Definition()
Definition G4Hydrogen.cc:45
G4MolecularConfiguration::GetDefinition
const G4MoleculeDefinition * GetDefinition() const
Definition G4MolecularConfiguration.hh:475
G4OH::Definition
static G4OH * Definition()
Definition G4OH.cc:45
G4Oxygen::Definition
static G4Oxygen * Definition()
Definition G4Oxygen.cc:44
DNA::Penetration::Kreipl2009::GetPenetration
static void GetPenetration(G4double energy, G4ThreeVector &displacement)
Definition G4DNAOneStepThermalizationModel.cc:157
DNA::Penetration::Meesungnoen2002_amorphous::GetPenetration
static void GetPenetration(G4double energy, G4ThreeVector &displacement)
Definition G4DNAOneStepThermalizationModel.cc:150
DNA::Penetration::Meesungnoen2002::GetPenetration
static void GetPenetration(G4double energy, G4ThreeVector &displacement)
Definition G4DNAOneStepThermalizationModel.cc:144
DNA::Penetration::Ritchie1994::GetPenetration
static void GetPenetration(G4double energy, G4ThreeVector &displacement)
Definition G4DNAOneStepThermalizationModel.cc:176
DNA::Penetration::Terrisol1990::GetPenetration
static void GetPenetration(G4double energy, G4ThreeVector &displacement)
Definition G4DNAOneStepThermalizationModel.cc:238