ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAWaterDissociationDisplacer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAWaterDissociationDisplacer.cc
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 "G4OH.hh"
45 #include "G4H3O.hh"
46 #include "G4Electron_aq.hh"
47 #include "G4H2O2.hh"
48 #include "Randomize.hh"
49 #include "G4Molecule.hh"
50 #include "G4MolecularConfiguration.hh"
51 #include "G4RandomDirection.hh"
52 #include "G4Exp.hh"
53 #include "G4UnitsTable.hh"
54 #include "G4EmParameters.hh"
55 #include "G4DNAOneStepThermalizationModel.hh"
56 
57 using namespace std;
58 
59 //------------------------------------------------------------------------------
60 
61 G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
62  Ionisation_DissociationDecay)
63 
64 G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
65  A1B1_DissociationDecay)
66 
67 G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
68  B1A1_DissociationDecay)
69 
70 G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
71  AutoIonisation)
72 
73 G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
74  DissociativeAttachment)
75 /*
76 //------------------------------------------------------------------------------
77 #ifdef _WATER_DISPLACER_USE_KREIPL_
78 // This function is used to generate the following density distribution:
79 // f(r) := 4*r * exp(-2*r)
80 // reference:
81 // Kreipl, M. S., Friedland, W. & Paretzke, H. G.
82 // Time- and space-resolved Monte Carlo study of water radiolysis for photon,
83 // electron and ion irradiation.
84 // Radiat. Environ. Biophys. 48, 11–20 (2009).
85 G4double G4DNAWaterDissociationDisplacer::ElectronProbaDistribution(G4double r)
86 {
87  return (2.*r+1.)*G4Exp(-2.*r);
88 }
89 #endif
90 
91 //------------------------------------------------------------------------------
92 #ifdef _WATER_DISPLACER_USE_TERRISOL_
93 // This function is used to generate the following density distribution:
94 // f(r) := 4*x^2/(sqrt(%pi)*(b)^3)*exp(-x^2/(b)^2)
95 // with b=27.22 for 7 eV
96 // reference
97 // Terrissol M, Beaudre A (1990) Simulation of space and time evolution
98 // of radiolytic species induced by electrons in water.
99 // Radiat Prot Dosimetry 31:171–175
100 
101 G4double G4DNAWaterDissociationDisplacer::ElectronProbaDistribution(G4double r)
102 {
103 #define b 27.22 // nanometer
104  static constexpr double sqrt_pi = 1.77245; // sqrt(CLHEP::pi);
105  static constexpr double b_to3 = 20168.1; // pow(b,3.);
106  static constexpr double b_to2 = 740.928; // pow(b,2.);
107  static constexpr double inverse_b_to2 = 1. / b_to2;
108 
109  static constexpr double main_factor = 4. / (sqrt_pi * b_to3);
110  static constexpr double factorA = sqrt_pi * b_to3 / 4.;
111  static constexpr double factorB = b_to2 / 2.;
112 
113  return (main_factor *
114  (factorA * erf(r / b)
115  - factorB * r * G4Exp(-pow(r, 2.) * inverse_b_to2)));
116 }
117 
118 #endif
119 //------------------------------------------------------------------------------
120 */
121 G4DNAWaterDissociationDisplacer::G4DNAWaterDissociationDisplacer()
122  :
123  G4VMolecularDissociationDisplacer(),
124  ke(1.7*eV)
125 /*#ifdef _WATER_DISPLACER_USE_KREIPL_
126  fFastElectronDistrib(0., 5., 0.001)
127 #elif defined _WATER_DISPLACER_USE_TERRISOL_
128  fFastElectronDistrib(0., 100., 0.001)
129 #endif*/
130 {/*
131  fProba1DFunction =
132  std::bind((G4double(*)(G4double))
133  &G4DNAWaterDissociationDisplacer::ElectronProbaDistribution,
134  std::placeholders::_1);
135 
136  size_t nBins = 500;
137  fElectronThermalization.reserve(nBins);
138  double eps = 1. / ((int) nBins);
139  double proba = eps;
140 
141  fElectronThermalization.push_back(0.);
142 
143  for (size_t i = 1; i < nBins; ++i)
144  {
145  double r = fFastElectronDistrib.Revert(proba, fProba1DFunction);
146  fElectronThermalization.push_back(r * nanometer);
147  proba += eps;
148 // G4cout << G4BestUnit(r*nanometer, "Length") << G4endl;
149  }*/
150  dnaSubType = G4EmParameters::Instance()->DNAeSolvationSubType();
151 // SetVerbose(1);
152 }
153 
154 //------------------------------------------------------------------------------
155 
156 G4DNAWaterDissociationDisplacer::~G4DNAWaterDissociationDisplacer()
157 {
158  ;
159 }
160 
161 //------------------------------------------------------------------------------
162 
163 G4ThreeVector
164 G4DNAWaterDissociationDisplacer::
165 GetMotherMoleculeDisplacement(const G4MolecularDissociationChannel*
166 theDecayChannel) const
167 {
168  G4int decayType = theDecayChannel->GetDisplacementType();
169  G4double RMSMotherMoleculeDisplacement(0.);
170 
171  switch (decayType)
172  {
173  case Ionisation_DissociationDecay:
174  RMSMotherMoleculeDisplacement = 2.0 * nanometer;
175  break;
176  case A1B1_DissociationDecay:
177  RMSMotherMoleculeDisplacement = 0. * nanometer;
178  break;
179  case B1A1_DissociationDecay:
180  RMSMotherMoleculeDisplacement = 0. * nanometer;
181  break;
182  case AutoIonisation:
183  RMSMotherMoleculeDisplacement = 2.0 * nanometer;
184  break;
185  case DissociativeAttachment:
186  RMSMotherMoleculeDisplacement = 0. * nanometer;
187  break;
188  }
189 
190  if (RMSMotherMoleculeDisplacement == 0)
191  {
192  return G4ThreeVector(0, 0, 0);
193  }
194  auto RandDirection =
195  radialDistributionOfProducts(RMSMotherMoleculeDisplacement);
196 
197  return RandDirection;
198 }
199 
200 //------------------------------------------------------------------------------
201 
202 vector<G4ThreeVector>
203 G4DNAWaterDissociationDisplacer::
204 GetProductsDisplacement(const G4MolecularDissociationChannel* pDecayChannel) const
205 {
206  G4int nbProducts = pDecayChannel->GetNbProducts();
207  vector<G4ThreeVector> theProductDisplacementVector(nbProducts);
208 
209  typedef map<const G4MoleculeDefinition*, G4double> RMSmap;
210  RMSmap theRMSmap;
211 
212  G4int decayType = pDecayChannel->GetDisplacementType();
213 
214  switch (decayType)
215  {
216  case Ionisation_DissociationDecay:
217  {
218  if (fVerbose)
219  {
220  G4cout << "Ionisation_DissociationDecay" << G4endl;
221  G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
222  }
223  G4double RdmValue = G4UniformRand();
224 
225  if (RdmValue < 0.5)
226  {
227  // H3O
228  theRMSmap[G4H3O::Definition()] = 0. * nanometer;
229  // OH
230  theRMSmap[G4OH::Definition()] = 0.8 * nanometer;
231  }
232  else
233  {
234  // H3O
235  theRMSmap[G4H3O::Definition()] = 0.8 * nanometer;
236  // OH
237  theRMSmap[G4OH::Definition()] = 0. * nanometer;
238  }
239 
240  for (int i = 0; i < nbProducts; i++)
241  {
242  auto pProduct = pDecayChannel->GetProduct(i);
243  G4double theRMSDisplacement = theRMSmap[pProduct->GetDefinition()];
244 
245  if (theRMSDisplacement == 0.)
246  {
247  theProductDisplacementVector[i] = G4ThreeVector();
248  }
249  else
250  {
251  auto RandDirection = radialDistributionOfProducts(theRMSDisplacement);
252  theProductDisplacementVector[i] = RandDirection;
253  }
254  }
255  break;
256  }
257  case A1B1_DissociationDecay:
258  {
259  if (fVerbose)
260  {
261  G4cout << "A1B1_DissociationDecay" << G4endl;
262  G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
263  }
264  G4double theRMSDisplacement = 2.4 * nanometer;
265  auto RandDirection =
266  radialDistributionOfProducts(theRMSDisplacement);
267 
268  for (G4int i = 0; i < nbProducts; i++)
269  {
270  auto pProduct = pDecayChannel->GetProduct(i);
271 
272  if (pProduct->GetDefinition() == G4OH::Definition())
273  {
274  theProductDisplacementVector[i] = -1. / 18. * RandDirection;
275  }
276  else if (pProduct->GetDefinition() == G4Hydrogen::Definition())
277  {
278  theProductDisplacementVector[i] = +17. / 18. * RandDirection;
279  }
280  }
281  break;
282  }
283  case B1A1_DissociationDecay:
284  {
285  if (fVerbose)
286  {
287  G4cout << "B1A1_DissociationDecay" << G4endl;
288  G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
289  }
290 
291  G4double theRMSDisplacement = 0.8 * nanometer;
292  auto RandDirection =
293  radialDistributionOfProducts(theRMSDisplacement);
294 
295  G4int NbOfOH = 0;
296  for (G4int i = 0; i < nbProducts; ++i)
297  {
298  auto pProduct = pDecayChannel->GetProduct(i);
299  if (pProduct->GetDefinition() == G4H2::Definition())
300  {
301  // H2
302  theProductDisplacementVector[i] = -2. / 18. * RandDirection;
303  }
304  else if (pProduct->GetDefinition() == G4OH::Definition())
305  {
306  // OH
307  G4ThreeVector OxygenDisplacement = +16. / 18. * RandDirection;
308  G4double OHRMSDisplacement = 1.1 * nanometer;
309 
310  auto OHDisplacement =
311  radialDistributionOfProducts(OHRMSDisplacement);
312 
313  if (NbOfOH == 0)
314  {
315  OHDisplacement = 0.5 * OHDisplacement;
316  }
317  else
318  {
319  OHDisplacement = -0.5 * OHDisplacement;
320  }
321 
322  theProductDisplacementVector[i] =
323  OHDisplacement + OxygenDisplacement;
324 
325  ++NbOfOH;
326  }
327  }
328  break;
329  }
330  case AutoIonisation:
331  {
332  if (fVerbose)
333  {
334  G4cout << "AutoIonisation" << G4endl;
335  G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
336  }
337 
338  G4double RdmValue = G4UniformRand();
339 
340  if (RdmValue < 0.5)
341  {
342  // H3O
343  theRMSmap[G4H3O::Definition()] = 0. * nanometer;
344  // OH
345  theRMSmap[G4OH::Definition()] = 0.8 * nanometer;
346  }
347  else
348  {
349  // H3O
350  theRMSmap[G4H3O::Definition()] = 0.8 * nanometer;
351  // OH
352  theRMSmap[G4OH::Definition()] = 0. * nanometer;
353  }
354 
355  for (G4int i = 0; i < nbProducts; i++)
356  {
357  auto pProduct = pDecayChannel->GetProduct(i);
358  auto theRMSDisplacement = theRMSmap[pProduct->GetDefinition()];
359 
360  if (theRMSDisplacement == 0)
361  {
362  theProductDisplacementVector[i] = G4ThreeVector();
363  }
364  else
365  {
366  auto RandDirection =
367  radialDistributionOfProducts(theRMSDisplacement);
368  theProductDisplacementVector[i] = RandDirection;
369  }
370  if (pProduct->GetDefinition() == G4Electron_aq::Definition())
371  {
372  theProductDisplacementVector[i] = radialDistributionOfElectron();
373  }
374  }
375  break;
376  }
377  case DissociativeAttachment:
378  {
379  if (fVerbose)
380  {
381  G4cout << "DissociativeAttachment" << G4endl;
382  G4cout << "Channel's name: " << pDecayChannel->GetName() << G4endl;
383  }
384  G4double theRMSDisplacement = 0.8 * nanometer;
385  auto RandDirection =
386  radialDistributionOfProducts(theRMSDisplacement);
387 
388  G4int NbOfOH = 0;
389  for (G4int i = 0; i < nbProducts; ++i)
390  {
391  auto pProduct = pDecayChannel->GetProduct(i);
392  if (pProduct->GetDefinition() == G4H2::Definition())
393  {
394  theProductDisplacementVector[i] = -2. / 18. * RandDirection;
395  }
396  else if (pProduct->GetDefinition() == G4OH::Definition())
397  {
398  G4ThreeVector OxygenDisplacement = +16. / 18. * RandDirection;
399  G4double OHRMSDisplacement = 1.1 * nanometer;
400 
401  auto OHDisplacement =
402  radialDistributionOfProducts(OHRMSDisplacement);
403 
404  if (NbOfOH == 0)
405  {
406  OHDisplacement = 0.5 * OHDisplacement;
407  }
408  else
409  {
410  OHDisplacement = -0.5 * OHDisplacement;
411  }
412  theProductDisplacementVector[i] = OHDisplacement +
413  OxygenDisplacement;
414  ++NbOfOH;
415  }
416  }
417  break;
418  }
419  }
420  return theProductDisplacementVector;
421 }
422 
423 //------------------------------------------------------------------------------
424 
425 G4ThreeVector
426 G4DNAWaterDissociationDisplacer::
427 radialDistributionOfProducts(G4double Rrms) const
428 {
429  static const double inverse_sqrt_3 = 1. / sqrt(3.);
430  double sigma = Rrms * inverse_sqrt_3;
431  double x = G4RandGauss::shoot(0., sigma);
432  double y = G4RandGauss::shoot(0., sigma);
433  double z = G4RandGauss::shoot(0., sigma);
434  return G4ThreeVector(x, y, z);
435 }
436 
437 //------------------------------------------------------------------------------
438 
439 G4ThreeVector
440 G4DNAWaterDissociationDisplacer::radialDistributionOfElectron() const
441 {/*
442  G4double rand_value = G4UniformRand();
443  size_t nBins = fElectronThermalization.size();
444  size_t bin = size_t(floor(rand_value * nBins));
445  size_t bin_p1 = min(bin + 1, nBins - 1);
446 
447  return (fElectronThermalization[bin] * (1. - rand_value)
448  + fElectronThermalization[bin_p1] * rand_value) *
449  G4RandomDirection();*/
450 
451  G4ThreeVector pdf = G4ThreeVector(0,0,0);
452 
453  if(dnaSubType == fRitchie1994eSolvation) DNA::Penetration::Ritchie1994::GetPenetration(ke,pdf);
454  else if(dnaSubType == fTerrisol1990eSolvation) DNA::Penetration::Terrisol1990::GetPenetration(ke,pdf);
455  else if(dnaSubType == fMeesungnoensolid2002eSolvation) DNA::Penetration::Meesungnoen2002_amorphous::GetPenetration(ke,pdf);
456  else if(dnaSubType == fKreipl2009eSolvation) DNA::Penetration::Kreipl2009::GetPenetration(ke,pdf);
457  else DNA::Penetration::Meesungnoen2002::GetPenetration(ke,pdf);
458 
459  return pdf;
460 }