| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /* 26 /*
27 * G4DNAIRT.cc 27 * G4DNAIRT.cc
28 * 28 *
29 * Created on: Jul 23, 2019 29 * Created on: Jul 23, 2019
30 * Author: W. G. Shin 30 * Author: W. G. Shin
31 * J. Ramos-Mendez and B. Faddego 31 * J. Ramos-Mendez and B. Faddegon
32 */ 32 */
33 33
34 34
35 #include "G4DNAIRT.hh" 35 #include "G4DNAIRT.hh"
36 #include "G4ErrorFunction.hh" 36 #include "G4ErrorFunction.hh"
37 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
38 #include "G4PhysicalConstants.hh" 38 #include "G4PhysicalConstants.hh"
39 #include "Randomize.hh" 39 #include "Randomize.hh"
40 #include "G4DNAMolecularReactionTable.hh" 40 #include "G4DNAMolecularReactionTable.hh"
41 #include "G4MolecularConfiguration.hh" 41 #include "G4MolecularConfiguration.hh"
42 #include "G4Molecule.hh" 42 #include "G4Molecule.hh"
43 #include "G4ITReactionChange.hh" 43 #include "G4ITReactionChange.hh"
44 #include "G4ITTrackHolder.hh" 44 #include "G4ITTrackHolder.hh"
45 #include "G4ITReaction.hh" 45 #include "G4ITReaction.hh"
46 #include "G4Scheduler.hh" 46 #include "G4Scheduler.hh"
47 47
48 using namespace std; 48 using namespace std;
49 49
50 G4DNAIRT::G4DNAIRT() : 50 G4DNAIRT::G4DNAIRT() :
>> 51 G4VITReactionProcess(),
51 fMolReactionTable(reference_cast<const G4DNAMo 52 fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)),
52 fpReactionModel(nullptr), 53 fpReactionModel(nullptr),
53 fTrackHolder(G4ITTrackHolder::Instance()), 54 fTrackHolder(G4ITTrackHolder::Instance()),
54 fReactionSet(nullptr) << 55 fReactionSet(0)
55 { 56 {
56 timeMin = G4Scheduler::Instance()->GetStartT 57 timeMin = G4Scheduler::Instance()->GetStartTime();
57 timeMax = G4Scheduler::Instance()->GetEndTim 58 timeMax = G4Scheduler::Instance()->GetEndTime();
58 59
59 fXMin = 1e9*nm; 60 fXMin = 1e9*nm;
60 fYMin = 1e9*nm; 61 fYMin = 1e9*nm;
61 fZMin = 1e9*nm; 62 fZMin = 1e9*nm;
62 63
63 fXMax = 0e0*nm; 64 fXMax = 0e0*nm;
64 fYMax = 0e0*nm; 65 fYMax = 0e0*nm;
65 fZMax = 0e0*nm; 66 fZMax = 0e0*nm;
66 67
67 fNx = 0; 68 fNx = 0;
68 fNy = 0; 69 fNy = 0;
69 fNz = 0; 70 fNz = 0;
70 71
71 xiniIndex = 0, yiniIndex = 0, ziniIndex = 0; 72 xiniIndex = 0, yiniIndex = 0, ziniIndex = 0;
72 xendIndex = 0, yendIndex = 0, zendIndex = 0; 73 xendIndex = 0, yendIndex = 0, zendIndex = 0;
73 74
74 fRCutOff = 75 fRCutOff =
75 1.45 * nm + 2 * std::sqrt(8*9.46e9*nm*nm/s 76 1.45 * nm + 2 * std::sqrt(8*9.46e9*nm*nm/s * timeMax); // 95% confidence level
76 77
77 erfc = new G4ErrorFunction(); 78 erfc = new G4ErrorFunction();
78 } 79 }
79 80
80 81
81 G4DNAIRT::G4DNAIRT(G4VDNAReactionModel* pReact 82 G4DNAIRT::G4DNAIRT(G4VDNAReactionModel* pReactionModel)
82 : G4DNAIRT() 83 : G4DNAIRT()
83 { 84 {
84 fpReactionModel = pReactionModel; 85 fpReactionModel = pReactionModel;
85 } 86 }
86 87
87 G4DNAIRT::~G4DNAIRT() 88 G4DNAIRT::~G4DNAIRT()
88 { 89 {
89 delete erfc; 90 delete erfc;
90 } 91 }
91 92
92 void G4DNAIRT::Initialize(){ 93 void G4DNAIRT::Initialize(){
93 94
94 fTrackHolder = G4ITTrackHolder::Instance(); 95 fTrackHolder = G4ITTrackHolder::Instance();
95 96
96 fReactionSet = G4ITReactionSet::Instance(); 97 fReactionSet = G4ITReactionSet::Instance();
97 fReactionSet->CleanAllReaction(); 98 fReactionSet->CleanAllReaction();
98 fReactionSet->SortByTime(); 99 fReactionSet->SortByTime();
99 100
100 spaceBinned.clear(); 101 spaceBinned.clear();
101 102
102 timeMin = G4Scheduler::Instance()->GetStartT 103 timeMin = G4Scheduler::Instance()->GetStartTime();
103 timeMax = G4Scheduler::Instance()->GetEndTim 104 timeMax = G4Scheduler::Instance()->GetEndTime();
104 105
105 xiniIndex = 0; 106 xiniIndex = 0;
106 yiniIndex = 0; 107 yiniIndex = 0;
107 ziniIndex = 0; 108 ziniIndex = 0;
108 xendIndex = 0; 109 xendIndex = 0;
109 yendIndex = 0; 110 yendIndex = 0;
110 zendIndex = 0; 111 zendIndex = 0;
111 112
112 fXMin = 1e9*nm; 113 fXMin = 1e9*nm;
113 fYMin = 1e9*nm; 114 fYMin = 1e9*nm;
114 fZMin = 1e9*nm; 115 fZMin = 1e9*nm;
115 116
116 fXMax = 0e0*nm; 117 fXMax = 0e0*nm;
117 fYMax = 0e0*nm; 118 fYMax = 0e0*nm;
118 fZMax = 0e0*nm; 119 fZMax = 0e0*nm;
119 120
120 fNx = 0; 121 fNx = 0;
121 fNy = 0; 122 fNy = 0;
122 fNz = 0; 123 fNz = 0;
123 124
124 SpaceBinning(); // 1. binning the space 125 SpaceBinning(); // 1. binning the space
125 IRTSampling(); // 2. Sampling of the IRT 126 IRTSampling(); // 2. Sampling of the IRT
126 127
127 //hoang : if the first IRTSampling won't giv <<
128 if(fReactionSet->Empty()) <<
129 { <<
130 for (auto pTrack : *fTrackHolder->GetMainL <<
131 { <<
132 pTrack->SetGlobalTime(G4Scheduler::Insta <<
133 } <<
134 } <<
135 } 128 }
136 129
137 void G4DNAIRT::SpaceBinning(){ 130 void G4DNAIRT::SpaceBinning(){
138 auto it_begin = fTrackHolder->GetMainList()- 131 auto it_begin = fTrackHolder->GetMainList()->begin();
139 while(it_begin != fTrackHolder->GetMainList( 132 while(it_begin != fTrackHolder->GetMainList()->end()){
140 133
141 G4ThreeVector position = it_begin->GetPosi 134 G4ThreeVector position = it_begin->GetPosition();
142 135
143 if ( fXMin > position.x() ) fXMin = positi 136 if ( fXMin > position.x() ) fXMin = position.x();
144 if ( fYMin > position.y() ) fYMin = positi 137 if ( fYMin > position.y() ) fYMin = position.y();
145 if ( fZMin > position.z() ) fZMin = positi 138 if ( fZMin > position.z() ) fZMin = position.z();
146 139
147 if ( fXMax < position.x() ) fXMax = positi 140 if ( fXMax < position.x() ) fXMax = position.x();
148 if ( fYMax < position.y() ) fYMax = positi 141 if ( fYMax < position.y() ) fYMax = position.y();
149 if ( fZMax < position.z() ) fZMax = positi 142 if ( fZMax < position.z() ) fZMax = position.z();
150 143
151 ++it_begin; 144 ++it_begin;
152 } 145 }
153 146
154 fNx = G4int((fXMax-fXMin)/fRCutOff) == 0 ? 1 147 fNx = G4int((fXMax-fXMin)/fRCutOff) == 0 ? 1 : G4int((fXMax-fXMin)/fRCutOff);
155 fNy = G4int((fYMax-fYMin)/fRCutOff) == 0 ? 1 148 fNy = G4int((fYMax-fYMin)/fRCutOff) == 0 ? 1 : G4int((fYMax-fYMin)/fRCutOff);
156 fNz = G4int((fZMax-fZMin)/fRCutOff) == 0 ? 1 149 fNz = G4int((fZMax-fZMin)/fRCutOff) == 0 ? 1 : G4int((fZMax-fZMin)/fRCutOff);
157 150
158 } 151 }
159 152
160 void G4DNAIRT::IRTSampling(){ 153 void G4DNAIRT::IRTSampling(){
161 154
162 auto it_begin = fTrackHolder->GetMainList()- 155 auto it_begin = fTrackHolder->GetMainList()->begin();
163 while(it_begin != fTrackHolder->GetMainList( 156 while(it_begin != fTrackHolder->GetMainList()->end()){
164 G4int I = FindBin(fNx, fXMin, fXMax, it_be 157 G4int I = FindBin(fNx, fXMin, fXMax, it_begin->GetPosition().x());
165 G4int J = FindBin(fNy, fYMin, fYMax, it_be 158 G4int J = FindBin(fNy, fYMin, fYMax, it_begin->GetPosition().y());
166 G4int K = FindBin(fNz, fZMin, fZMax, it_be 159 G4int K = FindBin(fNz, fZMin, fZMax, it_begin->GetPosition().z());
167 160
168 spaceBinned[I][J][K].push_back(*it_begin); 161 spaceBinned[I][J][K].push_back(*it_begin);
169 162
170 Sampling(*it_begin); 163 Sampling(*it_begin);
171 ++it_begin; 164 ++it_begin;
172 } 165 }
173 } 166 }
174 167
175 void G4DNAIRT::Sampling(G4Track* track){ 168 void G4DNAIRT::Sampling(G4Track* track){
176 G4Molecule* molA = G4Molecule::GetMolecule(t 169 G4Molecule* molA = G4Molecule::GetMolecule(track);
177 const G4MolecularConfiguration* molConfA = m 170 const G4MolecularConfiguration* molConfA = molA->GetMolecularConfiguration();
178 if(molConfA->GetDiffusionCoefficient() == 0) 171 if(molConfA->GetDiffusionCoefficient() == 0) return;
179 172
180 const vector<const G4MolecularConfiguration* 173 const vector<const G4MolecularConfiguration*>* reactivesVector =
181 fMolReactionTable->CanReactWith(molConfA); 174 fMolReactionTable->CanReactWith(molConfA);
182 175
183 if(reactivesVector == nullptr) return; 176 if(reactivesVector == nullptr) return;
184 177
185 G4double globalTime = G4Scheduler::Instance( 178 G4double globalTime = G4Scheduler::Instance()->GetGlobalTime();
186 G4double minTime = timeMax; 179 G4double minTime = timeMax;
187 180
188 xiniIndex = FindBin(fNx, fXMin, fXMax, track 181 xiniIndex = FindBin(fNx, fXMin, fXMax, track->GetPosition().x()-fRCutOff);
189 xendIndex = FindBin(fNx, fXMin, fXMax, track 182 xendIndex = FindBin(fNx, fXMin, fXMax, track->GetPosition().x()+fRCutOff);
190 yiniIndex = FindBin(fNy, fYMin, fYMax, track 183 yiniIndex = FindBin(fNy, fYMin, fYMax, track->GetPosition().y()-fRCutOff);
191 yendIndex = FindBin(fNy, fYMin, fYMax, track 184 yendIndex = FindBin(fNy, fYMin, fYMax, track->GetPosition().y()+fRCutOff);
192 ziniIndex = FindBin(fNz, fZMin, fZMax, track 185 ziniIndex = FindBin(fNz, fZMin, fZMax, track->GetPosition().z()-fRCutOff);
193 zendIndex = FindBin(fNz, fZMin, fZMax, track 186 zendIndex = FindBin(fNz, fZMin, fZMax, track->GetPosition().z()+fRCutOff);
194 187
195 for ( int ii = xiniIndex; ii <= xendIndex; i 188 for ( int ii = xiniIndex; ii <= xendIndex; ii++ ) {
196 for ( int jj = yiniIndex; jj <= yendIndex; 189 for ( int jj = yiniIndex; jj <= yendIndex; jj++ ) {
197 for ( int kk = ziniIndex; kk <= zendInde 190 for ( int kk = ziniIndex; kk <= zendIndex; kk++ ) {
198 191
199 std::vector<G4Track*> spaceBin = space 192 std::vector<G4Track*> spaceBin = spaceBinned[ii][jj][kk];
200 for ( int n = 0; n < (int)spaceBinned[ 193 for ( int n = 0; n < (int)spaceBinned[ii][jj][kk].size(); n++ ) {
201 if((spaceBin[n] == nullptr) || track << 194 if(!spaceBin[n] || track == spaceBin[n]) continue;
202 if(spaceBin[n]->GetTrackStatus() == 195 if(spaceBin[n]->GetTrackStatus() == fStopButAlive) continue;
203 196
204 G4Molecule* molB = G4Molecule::GetMo 197 G4Molecule* molB = G4Molecule::GetMolecule(spaceBin[n]);
205 if(molB == nullptr) continue; << 198 if(!molB) continue;
206 199
207 const G4MolecularConfiguration* molC 200 const G4MolecularConfiguration* molConfB = molB->GetMolecularConfiguration();
208 if(molConfB->GetDiffusionCoefficient 201 if(molConfB->GetDiffusionCoefficient() == 0) continue;
209 202
210 auto it = std::find(reactivesVector- 203 auto it = std::find(reactivesVector->begin(), reactivesVector->end(), molConfB);
211 if(it == reactivesVector->end()) con 204 if(it == reactivesVector->end()) continue;
212 205
213 G4ThreeVector orgPosB = spaceBin[n]- 206 G4ThreeVector orgPosB = spaceBin[n]->GetPosition();
214 G4double dt = track->GetGlobalTime() 207 G4double dt = track->GetGlobalTime() - spaceBin[n]->GetGlobalTime();
215 G4ThreeVector newPosB = orgPosB; 208 G4ThreeVector newPosB = orgPosB;
216 209
217 if(dt > 0){ 210 if(dt > 0){
218 G4double sigma, x, y, z; 211 G4double sigma, x, y, z;
219 G4double diffusionCoefficient = G4 212 G4double diffusionCoefficient = G4Molecule::GetMolecule(spaceBin[n])->GetDiffusionCoefficient();
220 213
221 sigma = std::sqrt(2.0 * diffusionC 214 sigma = std::sqrt(2.0 * diffusionCoefficient * dt);
222 215
223 x = G4RandGauss::shoot(0., 1.0)*si 216 x = G4RandGauss::shoot(0., 1.0)*sigma;
224 y = G4RandGauss::shoot(0., 1.0)*si 217 y = G4RandGauss::shoot(0., 1.0)*sigma;
225 z = G4RandGauss::shoot(0., 1.0)*si 218 z = G4RandGauss::shoot(0., 1.0)*sigma;
226 219
227 newPosB = orgPosB + G4ThreeVector( 220 newPosB = orgPosB + G4ThreeVector(x,y,z);
228 }else if(dt < 0) continue; 221 }else if(dt < 0) continue;
229 222
230 G4double r0 = (newPosB - track->GetP 223 G4double r0 = (newPosB - track->GetPosition()).mag();
231 G4double irt = GetIndependentReactio 224 G4double irt = GetIndependentReactionTime(molConfA,
232 225 molConfB,
233 226 r0);
234 if(irt>=0 && irt<timeMax - globalTim 227 if(irt>=0 && irt<timeMax - globalTime)
235 { 228 {
236 irt += globalTime; 229 irt += globalTime;
237 if(irt < minTime) minTime = irt; 230 if(irt < minTime) minTime = irt;
238 #ifdef DEBUG 231 #ifdef DEBUG
239 G4cout<<irt<<'\t'<<molConfA->GetNa 232 G4cout<<irt<<'\t'<<molConfA->GetName()<<" "<<track->GetTrackID()<<'\t'<<molConfB->GetName()<<" "<<spaceBin[n]->GetTrackID()<<'\n';
240 #endif 233 #endif
241 fReactionSet->AddReaction(irt,trac 234 fReactionSet->AddReaction(irt,track,spaceBin[n]);
242 } 235 }
243 } 236 }
244 spaceBin.clear(); 237 spaceBin.clear();
245 } 238 }
246 } 239 }
247 } 240 }
248 241
249 // Scavenging & first order reactions 242 // Scavenging & first order reactions
250 243
251 auto fReactionDatas = fMolReactionTable->Get 244 auto fReactionDatas = fMolReactionTable->GetReactionData(molConfA);
252 G4double index = -1; 245 G4double index = -1;
253 //change the scavenging filter of the IRT be <<
254 if(timeMax > 1*us) <<
255 { <<
256 minTime = timeMax; <<
257 } <<
258 // <<
259 246
260 for(size_t u=0; u<fReactionDatas->size();u++ 247 for(size_t u=0; u<fReactionDatas->size();u++){
261 if((*fReactionDatas)[u]->GetReactant2()->G 248 if((*fReactionDatas)[u]->GetReactant2()->GetDiffusionCoefficient() == 0){
262 G4double kObs = (*fReactionDatas)[u]->Ge 249 G4double kObs = (*fReactionDatas)[u]->GetObservedReactionRateConstant();
263 if(kObs == 0) continue; <<
264 G4double time = -(std::log(1.0 - G4Unifo 250 G4double time = -(std::log(1.0 - G4UniformRand())/kObs) + globalTime;
265 if( time < minTime && time >= globalTime 251 if( time < minTime && time >= globalTime && time < timeMax){
266 minTime = time; 252 minTime = time;
267 index = (G4int)u; << 253 index = (int) u;
268 } 254 }
269 } 255 }
270 } 256 }
271 257
272 if(index != -1){ 258 if(index != -1){
273 #ifdef DEBUG 259 #ifdef DEBUG
274 G4cout<<"scavenged: "<<minTime<<'\t'<<molC 260 G4cout<<"scavenged: "<<minTime<<'\t'<<molConfA->GetName()<<it_begin->GetTrackID()<<'\n';
275 #endif 261 #endif
276 auto fakeMol = new G4Molecule((*fReaction << 262 G4Molecule* fakeMol = new G4Molecule((*fReactionDatas)[index]->GetReactant2());
277 G4Track* fakeTrack = fakeMol->BuildTrack(g 263 G4Track* fakeTrack = fakeMol->BuildTrack(globalTime,track->GetPosition());
278 fTrackHolder->Push(fakeTrack); 264 fTrackHolder->Push(fakeTrack);
279 fReactionSet->AddReaction(minTime, track, 265 fReactionSet->AddReaction(minTime, track, fakeTrack);
280 } 266 }
281 } 267 }
282 268
283 269
284 G4double G4DNAIRT::GetIndependentReactionTime( 270 G4double G4DNAIRT::GetIndependentReactionTime(const G4MolecularConfiguration* molA, const G4MolecularConfiguration* molB, G4double distance) {
285 const auto pMoleculeA = molA; 271 const auto pMoleculeA = molA;
286 const auto pMoleculeB = molB; 272 const auto pMoleculeB = molB;
287 auto fReactionData = fMolReactionTable->GetR 273 auto fReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
288 G4int reactionType = fReactionData->GetReact 274 G4int reactionType = fReactionData->GetReactionType();
289 G4double r0 = distance; 275 G4double r0 = distance;
290 if(r0 == 0) r0 += 1e-3*nm; 276 if(r0 == 0) r0 += 1e-3*nm;
291 G4double irt = -1 * ps; 277 G4double irt = -1 * ps;
292 G4double D = molA->GetDiffusionCoefficient() 278 G4double D = molA->GetDiffusionCoefficient() +
293 molB->GetDiffusionCoefficient() 279 molB->GetDiffusionCoefficient();
294 if(D == 0) D += 1e-20*(m2/s); <<
295 G4double rc = fReactionData->GetOnsagerRadiu 280 G4double rc = fReactionData->GetOnsagerRadius();
296 281
297 if ( reactionType == 0){ 282 if ( reactionType == 0){
298 G4double sigma = fReactionData->GetEffecti 283 G4double sigma = fReactionData->GetEffectiveReactionRadius();
299 284
300 if(sigma > r0) return 0; // contact reacti <<
301 if( rc != 0) r0 = -rc / (1-std::exp(rc/r0) 285 if( rc != 0) r0 = -rc / (1-std::exp(rc/r0));
>> 286 if(sigma > r0) return 0; // contact reaction
302 287
303 G4double Winf = sigma/r0; 288 G4double Winf = sigma/r0;
304 G4double W = G4UniformRand(); 289 G4double W = G4UniformRand();
305 290
306 if ( W > 0 && W < Winf ) irt = (0.25/D) * << 291 if ( W < Winf ) irt = (0.25/D) * std::pow( (r0-sigma)/erfc->erfcInv(r0*W/sigma), 2 );
307 292
308 return irt; 293 return irt;
309 } 294 }
310 if ( reactionType == 1 ){ << 295 else if ( reactionType == 1 ){
311 G4double sigma = fReactionData->GetReactio 296 G4double sigma = fReactionData->GetReactionRadius();
312 G4double kact = fReactionData->GetActivati 297 G4double kact = fReactionData->GetActivationRateConstant();
313 G4double kdif = fReactionData->GetDiffusio 298 G4double kdif = fReactionData->GetDiffusionRateConstant();
314 G4double kobs = fReactionData->GetObserved 299 G4double kobs = fReactionData->GetObservedReactionRateConstant();
315 300
316 G4double a, b, Winf; 301 G4double a, b, Winf;
317 302
318 if ( rc == 0 ) { 303 if ( rc == 0 ) {
319 a = 1/sigma * kact / kobs; 304 a = 1/sigma * kact / kobs;
320 b = (r0 - sigma) / 2; 305 b = (r0 - sigma) / 2;
321 } else { 306 } else {
322 G4double v = kact/Avogadro/(4*CLHEP::pi* 307 G4double v = kact/Avogadro/(4*CLHEP::pi*pow(sigma,2) * exp(-rc / sigma));
323 G4double alpha = v+rc*D/(pow(sigma,2)*(1 308 G4double alpha = v+rc*D/(pow(sigma,2)*(1-exp(-rc/sigma)));
324 a = 4*pow(sigma,2)*alpha/(D*pow(rc,2))*p 309 a = 4*pow(sigma,2)*alpha/(D*pow(rc,2))*pow(sinh(rc/(2*sigma)),2);
325 b = rc/4*(cosh(rc/(2*r0))/sinh(rc/(2*r0) 310 b = rc/4*(cosh(rc/(2*r0))/sinh(rc/(2*r0))-cosh(rc/(2*sigma))/sinh(rc/(2*sigma)));
326 r0 = -rc/(1-std::exp(rc/r0)); 311 r0 = -rc/(1-std::exp(rc/r0));
327 sigma = fReactionData->GetEffectiveReact 312 sigma = fReactionData->GetEffectiveReactionRadius();
328 } 313 }
329 314
330 if(sigma > r0){ 315 if(sigma > r0){
331 if(fReactionData->GetProbability() > G4U 316 if(fReactionData->GetProbability() > G4UniformRand()) return 0;
332 return irt; << 317 else return irt;
333 } 318 }
334 Winf = sigma / r0 * kobs / kdif; 319 Winf = sigma / r0 * kobs / kdif;
335 320
336 if(Winf > G4UniformRand()) irt = SamplePDC 321 if(Winf > G4UniformRand()) irt = SamplePDC(a,b)/D;
337 return irt; << 322 return irt;
338 } << 323 }
339 324
340 return -1 * ps; 325 return -1 * ps;
341 } 326 }
342 327
343 G4int G4DNAIRT::FindBin(G4int n, G4double xmin 328 G4int G4DNAIRT::FindBin(G4int n, G4double xmin, G4double xmax, G4double value) {
344 329
345 G4int bin = -1; 330 G4int bin = -1;
346 if ( value <= xmin ) 331 if ( value <= xmin )
347 bin = 0; //1; 332 bin = 0; //1;
348 else if ( value >= xmax) //!(xmax < value) 333 else if ( value >= xmax) //!(xmax < value) ) //value >= xmax )
349 bin = n-1; //n; 334 bin = n-1; //n;
350 else 335 else
351 bin = G4int( n * ( value - xmin )/( xmax - 336 bin = G4int( n * ( value - xmin )/( xmax - xmin ) ); //bin = 1 + G4int( n * ( value - xmin )/( xmax - xmin ) );
352 337
353 if ( bin < 0 ) bin = 0; 338 if ( bin < 0 ) bin = 0;
354 if ( bin >= n ) bin = n-1; 339 if ( bin >= n ) bin = n-1;
355 340
356 return bin; 341 return bin;
357 } 342 }
358 343
359 G4double G4DNAIRT::SamplePDC(G4double a, G4dou 344 G4double G4DNAIRT::SamplePDC(G4double a, G4double b) {
360 345
361 G4double p = 2.0 * std::sqrt(2.0*b/a); 346 G4double p = 2.0 * std::sqrt(2.0*b/a);
362 G4double q = 2.0 / std::sqrt(2.0*b/a); 347 G4double q = 2.0 / std::sqrt(2.0*b/a);
363 G4double M = max(1.0/(a*a),3.0*b/a); 348 G4double M = max(1.0/(a*a),3.0*b/a);
364 349
365 G4double X, U, lambdax; 350 G4double X, U, lambdax;
366 351
367 G4int ntrials = 0; 352 G4int ntrials = 0;
368 while(true) { << 353 while(1) {
369 354
370 // Generate X 355 // Generate X
371 U = G4UniformRand(); 356 U = G4UniformRand();
372 if ( U < p/(p + q * M) ) X = pow(U * (p + 357 if ( U < p/(p + q * M) ) X = pow(U * (p + q * M) / 2, 2);
373 else X = pow(2/((1-U)*(p+q*M)/M),2); 358 else X = pow(2/((1-U)*(p+q*M)/M),2);
374 359
375 U = G4UniformRand(); 360 U = G4UniformRand();
376 361
377 lambdax = std::exp(-b*b/X) * ( 1.0 - a * s 362 lambdax = std::exp(-b*b/X) * ( 1.0 - a * std::sqrt(CLHEP::pi * X) * erfc->erfcx(b/std::sqrt(X) + a*std::sqrt(X)));
378 363
379 if ((X <= 2.0*b/a && U <= lambdax) || 364 if ((X <= 2.0*b/a && U <= lambdax) ||
380 (X >= 2.0*b/a && U*M/X <= lambdax)) br 365 (X >= 2.0*b/a && U*M/X <= lambdax)) break;
381 366
382 ntrials++; 367 ntrials++;
383 368
384 if ( ntrials > 10000 ){ 369 if ( ntrials > 10000 ){
385 G4cout<<"Totally rejected"<<'\n'; 370 G4cout<<"Totally rejected"<<'\n';
386 return -1.0; 371 return -1.0;
387 } 372 }
388 } 373 }
389 return X; 374 return X;
390 } 375 }
391 376
392 std::unique_ptr<G4ITReactionChange> G4DNAIRT:: 377 std::unique_ptr<G4ITReactionChange> G4DNAIRT::MakeReaction(const G4Track& trackA,
393 378 const G4Track& trackB)
394 { 379 {
395 380
396 std::unique_ptr<G4ITReactionChange> pChanges 381 std::unique_ptr<G4ITReactionChange> pChanges(new G4ITReactionChange());
397 pChanges->Initialize(trackA, trackB); 382 pChanges->Initialize(trackA, trackB);
398 383
399 const auto pMoleculeA = GetMolecule(trackA)- 384 const auto pMoleculeA = GetMolecule(trackA)->GetMolecularConfiguration();
400 const auto pMoleculeB = GetMolecule(trackB)- 385 const auto pMoleculeB = GetMolecule(trackB)->GetMolecularConfiguration();
401 const auto pReactionData = fMolReactionTable 386 const auto pReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
402 387
403 G4double globalTime = G4Scheduler::Instance( 388 G4double globalTime = G4Scheduler::Instance()->GetGlobalTime();
404 G4double effectiveReactionRadius = pReaction 389 G4double effectiveReactionRadius = pReactionData->GetEffectiveReactionRadius();
405 390
406 const G4double D1 = pMoleculeA->GetDiffusion 391 const G4double D1 = pMoleculeA->GetDiffusionCoefficient();
407 const G4double D2 = pMoleculeB->GetDiffusion 392 const G4double D2 = pMoleculeB->GetDiffusionCoefficient();
408 393
409 G4ThreeVector r1 = trackA.GetPosition(); 394 G4ThreeVector r1 = trackA.GetPosition();
410 G4ThreeVector r2 = trackB.GetPosition(); 395 G4ThreeVector r2 = trackB.GetPosition();
411 396
412 if(r1 == r2) r2 += G4ThreeVector(0,0,1e-3*nm 397 if(r1 == r2) r2 += G4ThreeVector(0,0,1e-3*nm);
413 398
414 G4ThreeVector S1 = r1 - r2; 399 G4ThreeVector S1 = r1 - r2;
415 400
416 G4double r0 = S1.mag(); 401 G4double r0 = S1.mag();
417 402
418 S1.setMag(effectiveReactionRadius); 403 S1.setMag(effectiveReactionRadius);
419 404
420 G4double dt = globalTime - trackA.GetGlobalT 405 G4double dt = globalTime - trackA.GetGlobalTime();
421 406
422 if(dt != 0 && (D1 + D2) != 0 && r0 != 0){ << 407 if(dt != 0){
423 G4double s12 = 2.0 * D1 * dt; 408 G4double s12 = 2.0 * D1 * dt;
424 G4double s22 = 2.0 * D2 * dt; 409 G4double s22 = 2.0 * D2 * dt;
425 if(s12 == 0) r2 = r1; 410 if(s12 == 0) r2 = r1;
426 else if(s22 == 0) r1 = r2; 411 else if(s22 == 0) r1 = r2;
427 else{ 412 else{
428 G4double alpha = effectiveReactionRadius 413 G4double alpha = effectiveReactionRadius * r0 / (2*(D1 + D2)*dt);
429 G4ThreeVector S2 = (r1 + (s12 / s22)*r2) 414 G4ThreeVector S2 = (r1 + (s12 / s22)*r2) + G4ThreeVector(G4RandGauss::shoot(0, s12 + s22 * s22 / s12),
430 415 G4RandGauss::shoot(0, s12 + s22 * s22 / s12),
431 416 G4RandGauss::shoot(0, s12 + s22 * s22 / s12));
432 417
433 if(alpha == 0){ <<
434 return pChanges; <<
435 } <<
436 S1.setPhi(rad * G4UniformRand() * 2.0 * 418 S1.setPhi(rad * G4UniformRand() * 2.0 * CLHEP::pi);
437 S1.setTheta(rad * std::acos(1.0 + 1./alp 419 S1.setTheta(rad * std::acos(1.0 + 1./alpha * std::log(1.0 - G4UniformRand() * (1 - std::exp(-2.0 * alpha)))));
438 420
439 r1 = (D1 * S1 + D2 * S2) / (D1 + D2); 421 r1 = (D1 * S1 + D2 * S2) / (D1 + D2);
440 r2 = D2 * (S2 - S1) / (D1 + D2); 422 r2 = D2 * (S2 - S1) / (D1 + D2);
441 } 423 }
442 } 424 }
443 425
444 auto pTrackA = const_cast<G4Track*>(pChanges 426 auto pTrackA = const_cast<G4Track*>(pChanges->GetTrackA());
445 auto pTrackB = const_cast<G4Track*>(pChanges 427 auto pTrackB = const_cast<G4Track*>(pChanges->GetTrackB());
446 428
447 pTrackA->SetPosition(r1); 429 pTrackA->SetPosition(r1);
448 pTrackB->SetPosition(r2); 430 pTrackB->SetPosition(r2);
449 431
450 pTrackA->SetGlobalTime(globalTime); 432 pTrackA->SetGlobalTime(globalTime);
451 pTrackB->SetGlobalTime(globalTime); 433 pTrackB->SetGlobalTime(globalTime);
452 434
453 pTrackA->SetTrackStatus(fStopButAlive); 435 pTrackA->SetTrackStatus(fStopButAlive);
454 pTrackB->SetTrackStatus(fStopButAlive); 436 pTrackB->SetTrackStatus(fStopButAlive);
455 437
456 const G4int nbProducts = pReactionData->GetN 438 const G4int nbProducts = pReactionData->GetNbProducts();
457 439
458 if(nbProducts != 0){ << 440 if(nbProducts){
459 441
460 const G4double sqrD1 = D1 == 0. ? 0. : std 442 const G4double sqrD1 = D1 == 0. ? 0. : std::sqrt(D1);
461 const G4double sqrD2 = D2 == 0. ? 0. : std 443 const G4double sqrD2 = D2 == 0. ? 0. : std::sqrt(D2);
462 if((sqrD1 + sqrD2) == 0){ <<
463 return pChanges; <<
464 } <<
465 const G4double inv_numerator = 1./(sqrD1 + 444 const G4double inv_numerator = 1./(sqrD1 + sqrD2);
466 const G4ThreeVector reactionSite = sqrD2 * 445 const G4ThreeVector reactionSite = sqrD2 * inv_numerator * trackA.GetPosition()
467 + sqrD1 * 446 + sqrD1 * inv_numerator * trackB.GetPosition();
468 447
469 std::vector<G4ThreeVector> position; 448 std::vector<G4ThreeVector> position;
470 449
471 if(nbProducts == 1){ 450 if(nbProducts == 1){
472 position.push_back(reactionSite); 451 position.push_back(reactionSite);
473 }else if(nbProducts == 2){ 452 }else if(nbProducts == 2){
474 position.push_back(trackA.GetPosition()) 453 position.push_back(trackA.GetPosition());
475 position.push_back(trackB.GetPosition()) 454 position.push_back(trackB.GetPosition());
476 }else if (nbProducts == 3){ 455 }else if (nbProducts == 3){
477 position.push_back(reactionSite); 456 position.push_back(reactionSite);
478 position.push_back(trackA.GetPosition()) 457 position.push_back(trackA.GetPosition());
479 position.push_back(trackB.GetPosition()) 458 position.push_back(trackB.GetPosition());
480 } 459 }
481 460
482 for(G4int u = 0; u < nbProducts; u++){ 461 for(G4int u = 0; u < nbProducts; u++){
483 462
484 auto product = new G4Molecule(pReactionD 463 auto product = new G4Molecule(pReactionData->GetProduct(u));
485 auto productTrack = product->BuildTrack( 464 auto productTrack = product->BuildTrack(globalTime,
486 465 position[u]);
487 466
488 productTrack->SetTrackStatus(fAlive); 467 productTrack->SetTrackStatus(fAlive);
489 fTrackHolder->Push(productTrack); 468 fTrackHolder->Push(productTrack);
490 469
491 pChanges->AddSecondary(productTrack); 470 pChanges->AddSecondary(productTrack);
492 471
493 G4int I = FindBin(fNx, fXMin, fXMax, pos 472 G4int I = FindBin(fNx, fXMin, fXMax, position[u].x());
494 G4int J = FindBin(fNy, fYMin, fYMax, pos 473 G4int J = FindBin(fNy, fYMin, fYMax, position[u].y());
495 G4int K = FindBin(fNz, fZMin, fZMax, pos 474 G4int K = FindBin(fNz, fZMin, fZMax, position[u].z());
496 475
497 spaceBinned[I][J][K].push_back(productTr 476 spaceBinned[I][J][K].push_back(productTrack);
498 477
499 Sampling(productTrack); 478 Sampling(productTrack);
500 } 479 }
501 } 480 }
502 481
503 fTrackHolder->MergeSecondariesWithMainList() 482 fTrackHolder->MergeSecondariesWithMainList();
504 pChanges->KillParents(true); 483 pChanges->KillParents(true);
505 return pChanges; 484 return pChanges;
506 } 485 }
507 486
508 487
509 std::vector<std::unique_ptr<G4ITReactionChange 488 std::vector<std::unique_ptr<G4ITReactionChange>> G4DNAIRT::FindReaction(
510 G4ITReactionSet* pReactionSet, 489 G4ITReactionSet* pReactionSet,
511 const G4double /*currentStepTime*/, << 490 const double /*currentStepTime*/,
512 const G4double fGlobalTime, << 491 const double fGlobalTime,
513 const G4bool /*reachedUserStepTimeLimit*/) << 492 const bool /*reachedUserStepTimeLimit*/)
514 { 493 {
515 std::vector<std::unique_ptr<G4ITReactionChan 494 std::vector<std::unique_ptr<G4ITReactionChange>> fReactionInfo;
516 fReactionInfo.clear(); 495 fReactionInfo.clear();
517 496
518 if (pReactionSet == nullptr) 497 if (pReactionSet == nullptr)
519 { 498 {
520 return fReactionInfo; 499 return fReactionInfo;
521 } 500 }
522 501
523 auto fReactionsetInTime = pReactionSet->GetR 502 auto fReactionsetInTime = pReactionSet->GetReactionsPerTime();
524 assert(fReactionsetInTime.begin() != fReacti 503 assert(fReactionsetInTime.begin() != fReactionsetInTime.end());
525 504
526 auto it_begin = fReactionsetInTime.begin(); 505 auto it_begin = fReactionsetInTime.begin();
527 while(it_begin != fReactionsetInTime.end()) 506 while(it_begin != fReactionsetInTime.end())
528 { 507 {
529 G4double irt = it_begin->get()->GetTime(); 508 G4double irt = it_begin->get()->GetTime();
530 509
531 if(fGlobalTime < irt) break; 510 if(fGlobalTime < irt) break;
532 511
533 pReactionSet->SelectThisReaction(*it_begin 512 pReactionSet->SelectThisReaction(*it_begin);
534 513
535 G4Track* pTrackA = it_begin->get()->GetRea 514 G4Track* pTrackA = it_begin->get()->GetReactants().first;
536 G4Track* pTrackB = it_begin->get()->GetRea 515 G4Track* pTrackB = it_begin->get()->GetReactants().second;
537 auto pReactionChange = MakeReaction(*pTrac 516 auto pReactionChange = MakeReaction(*pTrackA, *pTrackB);
538 517
539 if(pReactionChange){ 518 if(pReactionChange){
540 fReactionInfo.push_back(std::move(pReact 519 fReactionInfo.push_back(std::move(pReactionChange));
541 } 520 }
542 521
543 fReactionsetInTime = pReactionSet->GetReac 522 fReactionsetInTime = pReactionSet->GetReactionsPerTime();
544 it_begin = fReactionsetInTime.begin(); 523 it_begin = fReactionsetInTime.begin();
545 } 524 }
546 525
547 return fReactionInfo; 526 return fReactionInfo;
548 } 527 }
549 528
550 G4bool G4DNAIRT::TestReactibility(const G4Trac 529 G4bool G4DNAIRT::TestReactibility(const G4Track& /*trackA*/,
551 const G4Track& /*trackB*/, 530 const G4Track& /*trackB*/,
552 G4double /*currentStepTime*/, << 531 double /*currentStepTime*/,
553 G4bool /*userStepTimeLimit*/) /*const*/ << 532 bool /*userStepTimeLimit*/) /*const*/
554 { 533 {
555 return true; 534 return true;
556 } 535 }
557 536
558 void G4DNAIRT::SetReactionModel(G4VDNAReaction 537 void G4DNAIRT::SetReactionModel(G4VDNAReactionModel* model)
559 { 538 {
560 fpReactionModel = model; 539 fpReactionModel = model;
561 } 540 }
562 541