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