| Geant4 Cross Reference |
1 // ******************************************* 1
2 // * License and Disclaimer
3 // *
4 // * The Geant4 software is copyright of th
5 // * the Geant4 Collaboration. It is provided
6 // * conditions of the Geant4 Software License
7 // * LICENSE and available at http://cern.ch/
8 // * include a list of copyright holders.
9 // *
10 // * Neither the authors of this software syst
11 // * institutes,nor the agencies providing fin
12 // * work make any representation or warran
13 // * regarding this software system or assum
14 // * use. Please see the license in the file
15 // * for the full disclaimer and the limitatio
16 // *
17 // * This code implementation is the result
18 // * technical work of the GEANT4 collaboratio
19 // * By using, copying, modifying or distri
20 // * any work based on the software) you ag
21 // * use in resulting scientific publicati
22 // * acceptance of all terms of the Geant4 Sof
23 // *******************************************
24 //
25 //
26 // Author: Hoang TRAN
27
28 #include "G4DiffusionControlledReactionModel.h
29 #include "G4Track.hh"
30 #include "G4DNAMolecularReactionTable.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "G4Exp.hh"
33 #include "G4IRTUtils.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4Electron_aq.hh"
36 #include "Randomize.hh"
37 #include "G4Molecule.hh"
38 #include "G4ErrorFunction.hh"
39
40 G4DiffusionControlledReactionModel::G4Diffusio
41
42 G4DiffusionControlledReactionModel::~G4Diffusi
43
44 void G4DiffusionControlledReactionModel::Initi
45 const G4MolecularConfiguration* pMolecule, c
46 {
47 fpReactionData = fpReactionTable->GetReactio
48 }
49
50 void G4DiffusionControlledReactionModel::Initi
51 const G4MolecularConfiguration* pMolecule)
52 {
53 fpReactionData = fpReactionTable->GetReactio
54 }
55
56 G4double G4DiffusionControlledReactionModel::G
57 const G4MolecularConfiguration* pMol1, const
58 {
59 auto reactionData = fpReactionTable->GetReac
60 if(reactionData == nullptr)
61 {
62 G4ExceptionDescription exceptionDescriptio
63 exceptionDescription << "No reactionData"
64 << " for : " << pMol1
65 << pMol2->GetName();
66 G4Exception("G4DiffusionControlledReaction
67 "::GetReactionRadius()",
68 "G4DiffusionControlledReaction
69 exceptionDescription);
70 return 0.;
71 }
72 return reactionData->GetEffectiveReactionRad
73 }
74
75 G4double G4DiffusionControlledReactionModel::G
76 {
77 auto pMol1 = (*fpReactionData)[i]->GetReacta
78 auto pMol2 = (*fpReactionData)[i]->GetReacta
79 return GetReactionRadius(pMol1, pMol2);
80 }
81
82 G4double G4DiffusionControlledReactionModel::G
83 const G4Track& trackA, const G4Track& trackB
84 {
85 auto pMolConfA = GetMolecule(trackA)->GetMol
86 auto pMolConfB = GetMolecule(trackB)->GetMol
87
88 G4double D =
89 pMolConfA->GetDiffusionCoefficient() + pMo
90
91 if(D == 0)
92 {
93 G4ExceptionDescription exceptionDescriptio
94 exceptionDescription << "The total diffusi
95 << pMolConfA->GetName
96 << pMolConfB->GetName
97 G4Exception("G4DiffusionControlledReaction
98 "::GetTimeToEncounter()",
99 "G4DiffusionControlledReaction
100 exceptionDescription);
101 }
102
103 auto reactionData = G4DNAMolecularReactionTa
104 pMolConfA, pMolConfB);
105 G4double kobs = reactionData->GetObserve
106 G4double distance = (trackA.GetPosition() -
107 G4double SmoluchowskiRadius = reactionData->
108
109 if(distance == 0 || distance < SmoluchowskiR
110 {
111 G4ExceptionDescription exceptionDescriptio
112 exceptionDescription << "distance = " << d
113 << " Reff = " << Smol
114 << " for : " << pMolC
115 << pMolConfB->GetName
116 G4Exception("G4DiffusionControlledReaction
117 "::GetTimeToEncounter()",
118 "G4DiffusionControlledReaction
119 exceptionDescription);
120 }
121 else
122 {
123 G4double Winf = SmoluchowskiRadius / dist
124 G4double U = G4UniformRand();
125 G4double X = 0;
126 G4double irt_1 = -1.0 * ps;
127 if(Winf > 0 && U < Winf)
128 {
129 G4double erfcIn = G4ErrorFunction::erfcI
130 if(erfcIn != 0)
131 {
132 G4double d =
133 (distance - SmoluchowskiRadius) / er
134 irt_1 = (1.0 / (4 * D)) * d * d;
135 }
136 }
137
138 if(reactionData->GetReactionType() == 0)
139 {
140 return irt_1;
141 }
142
143 if(irt_1 < 0)
144 {
145 return irt_1;
146 }
147
148 G4double kdif = 4 * CLHEP::pi * D * Smoluc
149
150 if(pMolConfA == pMolConfB)
151 {
152 kdif /= 2;
153 }
154 G4double kact = G4IRTUtils::GetKact(kobs
155 G4double sumOfk = kact + kdif;
156 if(sumOfk != 0)
157 {
158 G4double rateFactor = kact / sumOfk;
159 if(G4UniformRand() > rateFactor)
160 {
161 return -1.0 * ps;
162 }
163 G4double Y = std::abs(G4RandGauss::shoot
164
165 if(Y > 0)
166 {
167 X = -(G4Log(G4UniformRand())) / Y;
168 }
169 G4double f = X * SmoluchowskiRadius
170 G4double irt_2 = (f * f) / D;
171 return irt_1 + irt_2;
172 }
173 }
174 return -1.0 * ps;
175 }
176