Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 // $Id: G4DNAWaterDissociationDisplacer.cc 101362 2016-11-15 15:26:23Z gcosmo $
27 //
28 // Author: Mathieu Karamitros
29 //
30 // WARNING : This class is released as a prototype.
31 // It might strongly evolve or even disapear in the next releases.
32 //
33 // History:
34 // -----------
35 // 10 Oct 2011 M.Karamitros created
36 //
37 // -------------------------------------------------------------------
38 
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4H2O.hh"
43 #include "G4H2.hh"
44 #include "G4Hydrogen.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"
52 #include "G4RandomDirection.hh"
53 #include "G4Exp.hh"
54 #include "G4UnitsTable.hh"
55 
56 using namespace std;
57 
58 //------------------------------------------------------------------------------
59 
61  Ionisation_DissociationDecay)
63  A1B1_DissociationDecay)
64 G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
65  B1A1_DissociationDecay)
66 G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
67  AutoIonisation)
68 G4CT_COUNT_IMPL(G4DNAWaterDissociationDisplacer,
69  DissociativeAttachment)
70 
71 //------------------------------------------------------------------------------
72 #ifdef _WATER_DISPLACER_USE_KREIPL_
73 // This function is used to generate the following density distribution:
74 // f(r) := 4*r * exp(-2*r)
75 // reference:
76 // Kreipl, M. S., Friedland, W. & Paretzke, H. G.
77 // Time- and space-resolved Monte Carlo study of water radiolysis for photon,
78 // electron and ion irradiation.
79 // Radiat. Environ. Biophys. 48, 11–20 (2009).
81 {
82  return (2.*r+1.)*G4Exp(-2.*r);
83 }
84 #endif
85 
86 //------------------------------------------------------------------------------
87 #ifdef _WATER_DISPLACER_USE_TERRISOL_
88 // This function is used to generate the following density distribution:
89 // f(r) := 4*x^2/(sqrt(%pi)*(b)^3)*exp(-x^2/(b)^2)
90 // with b=27.22 for 7 eV
91 // reference
92 // Terrissol M, Beaudre A (1990) Simulation of space and time evolution
93 // of radiolytic species induced by electrons in water.
94 // Radiat Prot Dosimetry 31:171–175
95 
97 {
98 #define b 27.22 //*nanometer
99  static constexpr double sqrt_pi=1.77245; // sqrt(CLHEP::pi);
100  static constexpr double b_to3 = 20168.1; // pow(b,3.);
101  static constexpr double b_to2 = 740.928; // pow(b,2.);
102  static constexpr double inverse_b_to2 = 1./b_to2;
103 
104  static constexpr double main_factor=4./(sqrt_pi*b_to3);
105  static constexpr double factorA=sqrt_pi*b_to3/4.;
106  static constexpr double factorB=b_to2/2.;
107 
108  return (main_factor*
109  (factorA*erf(r/b)
110  - factorB*r*G4Exp(-pow(r,2.)*inverse_b_to2)));
111 }
112 #endif
113 //------------------------------------------------------------------------------
114 
117 #ifdef _WATER_DISPLACER_USE_KREIPL_
118  fFastElectronDistrib(0., 5., 0.001)
120  fFastElectronDistrib(0., 100., 0.001)
121 #endif
122 {
123  fProba1DFunction =
124  std::bind((G4double(*)(G4double))
126  std::placeholders::_1);
127 
128  size_t nBins = 500;
129  fElectronThermalization.reserve(nBins);
130  double eps = 1./((int)nBins);
131  double proba = eps;
132 
133  fElectronThermalization.push_back(0.);
134 
135  for(size_t i = 1 ; i < nBins ; ++i){
136  double r = fFastElectronDistrib.Revert(proba, fProba1DFunction);
137  fElectronThermalization.push_back(r*nanometer);
138  proba+=eps;
139 // G4cout << G4BestUnit(r*nanometer, "Length") << G4endl;
140  }
141 // SetVerbose(1);
142 }
143 
144 //------------------------------------------------------------------------------
145 
147 {;}
148 
149 //------------------------------------------------------------------------------
150 
154  theDecayChannel) const
155 {
156  G4int decayType = theDecayChannel->GetDisplacementType();
157  G4double RMSMotherMoleculeDisplacement(0.);
158 
159  switch(decayType){
160  case Ionisation_DissociationDecay:
161  RMSMotherMoleculeDisplacement = 2.0 * nanometer;
162  break;
163  case A1B1_DissociationDecay:
164  RMSMotherMoleculeDisplacement = 0. * nanometer;
165  break;
166  case B1A1_DissociationDecay:
167  RMSMotherMoleculeDisplacement = 0. * nanometer;
168  break;
169  case AutoIonisation:
170  RMSMotherMoleculeDisplacement = 2.0 * nanometer;
171  break;
172  case DissociativeAttachment:
173  RMSMotherMoleculeDisplacement = 0. * nanometer;
174  break;
175  }
176 
177  if(RMSMotherMoleculeDisplacement == 0){
178  return G4ThreeVector(0, 0, 0);
179  }
180  auto RandDirection =
181  radialDistributionOfProducts(RMSMotherMoleculeDisplacement);
182 
183  return RandDirection;
184 }
185 
186 //------------------------------------------------------------------------------
187 
188 vector<G4ThreeVector>
191  theDecayChannel) const
192 {
193  G4int nbProducts = theDecayChannel->GetNbProducts();
194  vector<G4ThreeVector> theProductDisplacementVector(nbProducts);
195 
196  typedef map<const G4MoleculeDefinition*, G4double> RMSmap;
197  RMSmap theRMSmap;
198 
199  G4int decayType = theDecayChannel->GetDisplacementType();
200 
201  switch(decayType){
202  case Ionisation_DissociationDecay:
203  {
204  if (fVerbose){
205  G4cout << "Ionisation_DissociationDecay" << G4endl;
206  G4cout << "Channel's name: " << theDecayChannel->GetName() << G4endl;
207  }
208  G4double RdmValue = G4UniformRand();
209 
210  if(RdmValue< 0.5){
211  // H3O
212  theRMSmap[G4H3O::Definition()] = 0.* nanometer;
213  // OH
214  theRMSmap[G4OH::Definition()] = 0.8* nanometer;
215  }
216  else{
217  // H3O
218  theRMSmap[G4H3O::Definition()] = 0.8* nanometer;
219  // OH
220  theRMSmap[G4OH::Definition()] = 0.* nanometer;
221  }
222 
223  for(int i = 0; i < nbProducts; i++){
224  G4MolecularConfiguration* product = theDecayChannel->GetProduct(i);
225  G4double theRMSDisplacement = theRMSmap[product->GetDefinition()];
226 
227  if(theRMSDisplacement==0.){
228  theProductDisplacementVector[i] = G4ThreeVector();
229  }
230  else{
231  auto RandDirection = radialDistributionOfProducts(theRMSDisplacement);
232  theProductDisplacementVector[i] = RandDirection;
233  }
234  }
235  break;
236  }
237  case A1B1_DissociationDecay:
238  {
239  if(fVerbose){
240  G4cout<<"A1B1_DissociationDecay"<<G4endl;
241  G4cout << "Channel's name: " << theDecayChannel->GetName() << G4endl;
242  }
243  G4double theRMSDisplacement = 2.4 * nanometer;
244  auto RandDirection =
245  radialDistributionOfProducts(theRMSDisplacement);
246 
247  for(G4int i =0; i < nbProducts; i++){
248  G4MolecularConfiguration* product = theDecayChannel->GetProduct(i);
249 
250  if(product->GetDefinition()== G4OH::Definition()){
251  // OH
252  theProductDisplacementVector[i] = -1./18.*RandDirection;
253  }
254  else if(product->GetDefinition() == G4Hydrogen::Definition()){
255  // H
256  theProductDisplacementVector[i] = +17./18.*RandDirection;
257  }
258  }
259  break;
260  }
261  case B1A1_DissociationDecay:
262  {
263  if(fVerbose){
264  G4cout<<"B1A1_DissociationDecay"<<G4endl;
265  G4cout << "Channel's name: " << theDecayChannel->GetName() << G4endl;
266  }
267 
268  G4double theRMSDisplacement = 0.8 * nanometer;
269  auto RandDirection =
270  radialDistributionOfProducts(theRMSDisplacement);
271 
272  G4int NbOfOH = 0;
273  for(G4int i =0; i < nbProducts; ++i)
274  {
275  G4MolecularConfiguration* product = theDecayChannel->GetProduct(i);
276  if(product->GetDefinition() == G4H2::Definition()){
277  // H2
278  theProductDisplacementVector[i] = -2./18.*RandDirection;
279  }
280  else if(product->GetDefinition() == G4OH::Definition()){
281  // OH
282  G4ThreeVector OxygenDisplacement = +16./18.*RandDirection;
283  G4double OHRMSDisplacement = 1.1 * nanometer;
284 
285  auto OHDisplacement =
286  radialDistributionOfProducts(OHRMSDisplacement);
287 
288  if(NbOfOH==0){
289  OHDisplacement = 0.5*OHDisplacement;
290  }
291  else{
292  OHDisplacement = -0.5*OHDisplacement;
293  }
294 
295  theProductDisplacementVector[i] =
296  OHDisplacement + OxygenDisplacement;
297 
298  ++NbOfOH;
299  }
300  }
301  break;
302  }
303  case AutoIonisation:
304  {
305  if(fVerbose){
306  G4cout<<"AutoIonisation"<<G4endl;
307  G4cout << "Channel's name: " << theDecayChannel->GetName() << G4endl;
308  }
309 
310  G4double RdmValue = G4UniformRand();
311 
312  if(RdmValue< 0.5){
313  // H3O
314  theRMSmap[G4H3O::Definition()] = 0.* nanometer;
315  // OH
316  theRMSmap[G4OH::Definition()] = 0.8* nanometer;
317  }
318  else{
319  // H3O
320  theRMSmap[G4H3O::Definition()] = 0.8* nanometer;
321  // OH
322  theRMSmap[G4OH::Definition()] = 0.* nanometer;
323  }
324 
325  for(G4int i =0; i < nbProducts; i++){
326  G4MolecularConfiguration* product = theDecayChannel->GetProduct(i);
327  auto theRMSDisplacement = theRMSmap[product->GetDefinition()];
328 
329  if(theRMSDisplacement==0){
330  theProductDisplacementVector[i] = G4ThreeVector();
331  }
332  else{
333  auto RandDirection =
334  radialDistributionOfProducts(theRMSDisplacement);
335  theProductDisplacementVector[i] = RandDirection;
336  }
337  if(product->GetDefinition() == G4Electron_aq::Definition()){
338  theProductDisplacementVector[i]=radialDistributionOfElectron();
339  }
340  }
341  break;
342  }
343  case DissociativeAttachment:
344  {
345  if(fVerbose){
346  G4cout<<"DissociativeAttachment"<<G4endl;
347  G4cout << "Channel's name: " << theDecayChannel->GetName() << G4endl;
348  }
349  G4double theRMSDisplacement = 0.8 * nanometer;
350  auto RandDirection =
351  radialDistributionOfProducts(theRMSDisplacement);
352 
353  G4int NbOfOH = 0;
354  for(G4int i =0; i < nbProducts; ++i){
355  G4MolecularConfiguration* product = theDecayChannel->GetProduct(i);
356  if(product->GetDefinition() == G4H2::Definition()){
357  theProductDisplacementVector[i] = -2./18.*RandDirection;
358  }
359  else if(product->GetDefinition() == G4OH::Definition()){
360  G4ThreeVector OxygenDisplacement = +16./18.*RandDirection;
361  G4double OHRMSDisplacement = 1.1 * nanometer;
362 
363  auto OHDisplacement =
364  radialDistributionOfProducts(OHRMSDisplacement);
365 
366  if(NbOfOH==0){
367  OHDisplacement = 0.5*OHDisplacement;
368  }
369  else{
370  OHDisplacement = -0.5*OHDisplacement;
371  }
372  theProductDisplacementVector[i] = OHDisplacement +
373  OxygenDisplacement;
374  ++NbOfOH;
375  }
376  }
377  break;
378  }
379  }
380  return theProductDisplacementVector;
381 }
382 
383 //------------------------------------------------------------------------------
384 
388 {
389  static const double inverse_sqrt_3 = 1./sqrt(3.);
390  double sigma = Rrms*inverse_sqrt_3;
391  double x = G4RandGauss::shoot(0.,sigma);
392  double y = G4RandGauss::shoot(0.,sigma);
393  double z = G4RandGauss::shoot(0.,sigma);
394  return G4ThreeVector(x,y,z);
395 }
396 
397 //------------------------------------------------------------------------------
398 
401 {
402  G4double rand_value = G4UniformRand();
403  size_t nBins = fElectronThermalization.size();
404  size_t bin = size_t(floor(rand_value*nBins));
405  size_t bin_p1 = min(bin+1,nBins-1);
406 
407  return (fElectronThermalization[bin]*(1.-rand_value)
408  + fElectronThermalization[bin_p1]*rand_value)*
410 // return fElectronThermalization[bin]*G4RandomDirection();
411 }
ThreeVector shoot(const G4int Ap, const G4int Af)
static G4Electron_aq * Definition()
tuple bin
Definition: plottest35.py:22
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector G4RandomDirection()
static constexpr double nanometer
Definition: G4SIunits.hh:101
static const G4double eps
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
static G4double ElectronProbaDistribution(G4double r)
G4MolecularConfiguration * GetProduct(int) const
static G4H3O * Definition()
Definition: G4H3O.cc:47
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
#define _WATER_DISPLACER_USE_TERRISOL_
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
std::vector< G4ThreeVector > GetProductsDisplacement(const G4MolecularDissociationChannel *) const override
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
const G4MoleculeDefinition * GetDefinition() const
tuple z
Definition: test.py:28
G4ThreeVector GetMotherMoleculeDisplacement(const G4MolecularDissociationChannel *) const override
#define G4endl
Definition: G4ios.hh:61
static G4OH * Definition()
Definition: G4OH.cc:46
double G4double
Definition: G4Types.hh:76
static G4H2 * Definition()
Definition: G4H2.cc:46
G4ThreeVector radialDistributionOfProducts(G4double r_rms) const
static G4Hydrogen * Definition()
Definition: G4Hydrogen.cc:46
#define G4CT_COUNT_IMPL(enumName, flagName)
Definition: G4CTCounter.hh:125