Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNAIRT Class Reference

#include <G4DNAIRT.hh>

+ Inheritance diagram for G4DNAIRT:

Public Member Functions

 G4DNAIRT ()
 
 G4DNAIRT (G4VDNAReactionModel *)
 
 ~G4DNAIRT () override
 
 G4DNAIRT (const G4DNAIRT &other)=delete
 
G4DNAIRToperator= (const G4DNAIRT &other)=delete
 
G4bool TestReactibility (const G4Track &, const G4Track &, G4double, G4bool) override
 
std::vector< std::unique_ptr< G4ITReactionChange > > FindReaction (G4ITReactionSet *, const G4double, const G4double, const G4bool) override
 
std::unique_ptr< G4ITReactionChangeMakeReaction (const G4Track &, const G4Track &) override
 
void SetReactionModel (G4VDNAReactionModel *)
 
void Initialize () override
 
void SpaceBinning ()
 
void IRTSampling ()
 
void Sampling (G4Track *)
 
G4double GetIndependentReactionTime (const G4MolecularConfiguration *, const G4MolecularConfiguration *, G4double)
 
G4int FindBin (G4int, G4double, G4double, G4double)
 
G4double SamplePDC (G4double, G4double)
 
- Public Member Functions inherited from G4VITReactionProcess
 G4VITReactionProcess ()=default
 
virtual ~G4VITReactionProcess ()=default
 
 G4VITReactionProcess (const G4VITReactionProcess &other)=delete
 
G4VITReactionProcessoperator= (const G4VITReactionProcess &other)=delete
 
virtual void Initialize ()
 
virtual G4bool IsApplicable (const G4ITType &, const G4ITType &) const
 
virtual G4bool TestReactibility (const G4Track &, const G4Track &, double, bool)=0
 
virtual std::vector< std::unique_ptr< G4ITReactionChange > > FindReaction (G4ITReactionSet *, const double, const double, const bool)=0
 
virtual std::unique_ptr< G4ITReactionChangeMakeReaction (const G4Track &, const G4Track &)=0
 
virtual void SetReactionTable (const G4ITReactionTable *)
 

Protected Attributes

const G4DNAMolecularReactionTable *& fMolReactionTable
 
G4VDNAReactionModelfpReactionModel
 
- Protected Attributes inherited from G4VITReactionProcess
const G4ITReactionTablefpReactionTable = nullptr
 
G4String fName
 

Detailed Description

Definition at line 64 of file G4DNAIRT.hh.

Constructor & Destructor Documentation

◆ G4DNAIRT() [1/3]

G4DNAIRT::G4DNAIRT ( )

Definition at line 50 of file G4DNAIRT.cc.

50 :
52fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)),
53fpReactionModel(nullptr),
54fTrackHolder(G4ITTrackHolder::Instance()),
55fReactionSet(nullptr)
56{
58 timeMax = G4Scheduler::Instance()->GetEndTime();
59
60 fXMin = 1e9*nm;
61 fYMin = 1e9*nm;
62 fZMin = 1e9*nm;
63
64 fXMax = 0e0*nm;
65 fYMax = 0e0*nm;
66 fZMax = 0e0*nm;
67
68 fNx = 0;
69 fNy = 0;
70 fNz = 0;
71
72 xiniIndex = 0, yiniIndex = 0, ziniIndex = 0;
73 xendIndex = 0, yendIndex = 0, zendIndex = 0;
74
75 fRCutOff =
76 1.45 * nm + 2 * std::sqrt(8*9.46e9*nm*nm/s * timeMax); // 95% confidence level
77
78 erfc = new G4ErrorFunction();
79}
const G4DNAMolecularReactionTable *& fMolReactionTable
Definition: G4DNAIRT.hh:93
G4VDNAReactionModel * fpReactionModel
Definition: G4DNAIRT.hh:94
static G4ITTrackHolder * Instance()
static G4Scheduler * Instance()
Definition: G4Scheduler.cc:101
G4double GetEndTime() const
Definition: G4Scheduler.hh:349
G4double GetStartTime() const
Definition: G4Scheduler.hh:344
const G4ITReactionTable * fpReactionTable
G4VITReactionProcess()=default

◆ G4DNAIRT() [2/3]

G4DNAIRT::G4DNAIRT ( G4VDNAReactionModel pReactionModel)
explicit

Definition at line 82 of file G4DNAIRT.cc.

83 : G4DNAIRT()
84{
85 fpReactionModel = pReactionModel;
86}
G4DNAIRT()
Definition: G4DNAIRT.cc:50

◆ ~G4DNAIRT()

G4DNAIRT::~G4DNAIRT ( )
override

Definition at line 88 of file G4DNAIRT.cc.

89{
90 delete erfc;
91}

◆ G4DNAIRT() [3/3]

G4DNAIRT::G4DNAIRT ( const G4DNAIRT other)
delete

Member Function Documentation

◆ FindBin()

G4int G4DNAIRT::FindBin ( G4int  n,
G4double  xmin,
G4double  xmax,
G4double  value 
)

(xmax < value) ) //value >= xmax )

Definition at line 344 of file G4DNAIRT.cc.

344 {
345
346 G4int bin = -1;
347 if ( value <= xmin )
348 bin = 0; //1;
349 else if ( value >= xmax) //!(xmax < value) ) //value >= xmax )
350 bin = n-1; //n;
351 else
352 bin = G4int( n * ( value - xmin )/( xmax - xmin ) ); //bin = 1 + G4int( n * ( value - xmin )/( xmax - xmin ) );
353
354 if ( bin < 0 ) bin = 0;
355 if ( bin >= n ) bin = n-1;
356
357 return bin;
358}
int G4int
Definition: G4Types.hh:85

Referenced by IRTSampling(), MakeReaction(), and Sampling().

◆ FindReaction()

std::vector< std::unique_ptr< G4ITReactionChange > > G4DNAIRT::FindReaction ( G4ITReactionSet pReactionSet,
const  G4double,
const G4double  fGlobalTime,
const  G4bool 
)
overridevirtual

Implements G4VITReactionProcess.

Definition at line 510 of file G4DNAIRT.cc.

515{
516 std::vector<std::unique_ptr<G4ITReactionChange>> fReactionInfo;
517 fReactionInfo.clear();
518
519 if (pReactionSet == nullptr)
520 {
521 return fReactionInfo;
522 }
523
524 auto fReactionsetInTime = pReactionSet->GetReactionsPerTime();
525 assert(fReactionsetInTime.begin() != fReactionsetInTime.end());
526
527 auto it_begin = fReactionsetInTime.begin();
528 while(it_begin != fReactionsetInTime.end())
529 {
530 G4double irt = it_begin->get()->GetTime();
531
532 if(fGlobalTime < irt) break;
533
534 pReactionSet->SelectThisReaction(*it_begin);
535
536 G4Track* pTrackA = it_begin->get()->GetReactants().first;
537 G4Track* pTrackB = it_begin->get()->GetReactants().second;
538 auto pReactionChange = MakeReaction(*pTrackA, *pTrackB);
539
540 if(pReactionChange){
541 fReactionInfo.push_back(std::move(pReactionChange));
542 }
543
544 fReactionsetInTime = pReactionSet->GetReactionsPerTime();
545 it_begin = fReactionsetInTime.begin();
546 }
547
548 return fReactionInfo;
549}
double G4double
Definition: G4Types.hh:83
std::unique_ptr< G4ITReactionChange > MakeReaction(const G4Track &, const G4Track &) override
Definition: G4DNAIRT.cc:393
void SelectThisReaction(G4ITReactionPtr reaction)
G4ITReactionPerTime & GetReactionsPerTime()

◆ GetIndependentReactionTime()

G4double G4DNAIRT::GetIndependentReactionTime ( const G4MolecularConfiguration molA,
const G4MolecularConfiguration molB,
G4double  distance 
)

Definition at line 285 of file G4DNAIRT.cc.

285 {
286 const auto pMoleculeA = molA;
287 const auto pMoleculeB = molB;
288 auto fReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
289 G4int reactionType = fReactionData->GetReactionType();
290 G4double r0 = distance;
291 if(r0 == 0) r0 += 1e-3*nm;
292 G4double irt = -1 * ps;
295 if(D == 0) D += 1e-20*(m2/s);
296 G4double rc = fReactionData->GetOnsagerRadius();
297
298 if ( reactionType == 0){
299 G4double sigma = fReactionData->GetEffectiveReactionRadius();
300
301 if(sigma > r0) return 0; // contact reaction
302 if( rc != 0) r0 = -rc / (1-std::exp(rc/r0));
303
304 G4double Winf = sigma/r0;
306
307 if ( W > 0 && W < Winf ) irt = (0.25/D) * std::pow( (r0-sigma)/erfc->erfcInv(r0*W/sigma), 2 );
308
309 return irt;
310 }
311 else if ( reactionType == 1 ){
312 G4double sigma = fReactionData->GetReactionRadius();
313 G4double kact = fReactionData->GetActivationRateConstant();
314 G4double kdif = fReactionData->GetDiffusionRateConstant();
315 G4double kobs = fReactionData->GetObservedReactionRateConstant();
316
317 G4double a, b, Winf;
318
319 if ( rc == 0 ) {
320 a = 1/sigma * kact / kobs;
321 b = (r0 - sigma) / 2;
322 } else {
323 G4double v = kact/Avogadro/(4*CLHEP::pi*pow(sigma,2) * exp(-rc / sigma));
324 G4double alpha = v+rc*D/(pow(sigma,2)*(1-exp(-rc/sigma)));
325 a = 4*pow(sigma,2)*alpha/(D*pow(rc,2))*pow(sinh(rc/(2*sigma)),2);
326 b = rc/4*(cosh(rc/(2*r0))/sinh(rc/(2*r0))-cosh(rc/(2*sigma))/sinh(rc/(2*sigma)));
327 r0 = -rc/(1-std::exp(rc/r0));
328 sigma = fReactionData->GetEffectiveReactionRadius();
329 }
330
331 if(sigma > r0){
332 if(fReactionData->GetProbability() > G4UniformRand()) return 0;
333 else return irt;
334 }
335 Winf = sigma / r0 * kobs / kdif;
336
337 if(Winf > G4UniformRand()) irt = SamplePDC(a,b)/D;
338 return irt;
339 }
340
341 return -1 * ps;
342}
G4double D(G4double temp)
#define G4UniformRand()
Definition: Randomize.hh:52
G4double SamplePDC(G4double, G4double)
Definition: G4DNAIRT.cc:360
Data * GetReactionData(Reactant *, Reactant *) const
static G4double erfcInv(G4double x)
#define W
Definition: crc32.c:84

Referenced by Sampling().

◆ Initialize()

void G4DNAIRT::Initialize ( )
overridevirtual

First initialization (done once for all at the begin of the run) eg. check if the reaction table is given ...

Reimplemented from G4VITReactionProcess.

Definition at line 93 of file G4DNAIRT.cc.

93 {
94
95 fTrackHolder = G4ITTrackHolder::Instance();
96
97 fReactionSet = G4ITReactionSet::Instance();
98 fReactionSet->CleanAllReaction();
99 fReactionSet->SortByTime();
100
101 spaceBinned.clear();
102
103 timeMin = G4Scheduler::Instance()->GetStartTime();
104 timeMax = G4Scheduler::Instance()->GetEndTime();
105
106 xiniIndex = 0;
107 yiniIndex = 0;
108 ziniIndex = 0;
109 xendIndex = 0;
110 yendIndex = 0;
111 zendIndex = 0;
112
113 fXMin = 1e9*nm;
114 fYMin = 1e9*nm;
115 fZMin = 1e9*nm;
116
117 fXMax = 0e0*nm;
118 fYMax = 0e0*nm;
119 fZMax = 0e0*nm;
120
121 fNx = 0;
122 fNy = 0;
123 fNz = 0;
124
125 SpaceBinning(); // 1. binning the space
126 IRTSampling(); // 2. Sampling of the IRT
127
128 //hoang : if the first IRTSampling won't give any reactions, end the simu.
129 if(fReactionSet->Empty())
130 {
131 for (auto pTrack : *fTrackHolder->GetMainList())
132 {
133 pTrack->SetGlobalTime(G4Scheduler::Instance()->GetEndTime());
134 }
135 }
136}
void IRTSampling()
Definition: G4DNAIRT.cc:161
void SpaceBinning()
Definition: G4DNAIRT.cc:138
static G4ITReactionSet * Instance()
void CleanAllReaction()
G4TrackList * GetMainList(Key)

◆ IRTSampling()

void G4DNAIRT::IRTSampling ( )

Definition at line 161 of file G4DNAIRT.cc.

161 {
162
163 auto it_begin = fTrackHolder->GetMainList()->begin();
164 while(it_begin != fTrackHolder->GetMainList()->end()){
165 G4int I = FindBin(fNx, fXMin, fXMax, it_begin->GetPosition().x());
166 G4int J = FindBin(fNy, fYMin, fYMax, it_begin->GetPosition().y());
167 G4int K = FindBin(fNz, fZMin, fZMax, it_begin->GetPosition().z());
168
169 spaceBinned[I][J][K].push_back(*it_begin);
170
171 Sampling(*it_begin);
172 ++it_begin;
173 }
174}
G4int FindBin(G4int, G4double, G4double, G4double)
Definition: G4DNAIRT.cc:344
void Sampling(G4Track *)
Definition: G4DNAIRT.cc:176
iterator begin()
iterator end()

Referenced by Initialize().

◆ MakeReaction()

std::unique_ptr< G4ITReactionChange > G4DNAIRT::MakeReaction ( const G4Track trackA,
const G4Track trackB 
)
overridevirtual

Implements G4VITReactionProcess.

Definition at line 393 of file G4DNAIRT.cc.

395{
396
397 std::unique_ptr<G4ITReactionChange> pChanges(new G4ITReactionChange());
398 pChanges->Initialize(trackA, trackB);
399
400 const auto pMoleculeA = GetMolecule(trackA)->GetMolecularConfiguration();
401 const auto pMoleculeB = GetMolecule(trackB)->GetMolecularConfiguration();
402 const auto pReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
403
405 G4double effectiveReactionRadius = pReactionData->GetEffectiveReactionRadius();
406
407 const G4double D1 = pMoleculeA->GetDiffusionCoefficient();
408 const G4double D2 = pMoleculeB->GetDiffusionCoefficient();
409
410 G4ThreeVector r1 = trackA.GetPosition();
411 G4ThreeVector r2 = trackB.GetPosition();
412
413 if(r1 == r2) r2 += G4ThreeVector(0,0,1e-3*nm);
414
415 G4ThreeVector S1 = r1 - r2;
416
417 G4double r0 = S1.mag();
418
419 S1.setMag(effectiveReactionRadius);
420
421 G4double dt = globalTime - trackA.GetGlobalTime();
422
423 if(dt != 0 && (D1 + D2) != 0 && r0 != 0){
424 G4double s12 = 2.0 * D1 * dt;
425 G4double s22 = 2.0 * D2 * dt;
426 if(s12 == 0) r2 = r1;
427 else if(s22 == 0) r1 = r2;
428 else{
429 G4double alpha = effectiveReactionRadius * r0 / (2*(D1 + D2)*dt);
430 G4ThreeVector S2 = (r1 + (s12 / s22)*r2) + G4ThreeVector(G4RandGauss::shoot(0, s12 + s22 * s22 / s12),
431 G4RandGauss::shoot(0, s12 + s22 * s22 / s12),
432 G4RandGauss::shoot(0, s12 + s22 * s22 / s12));
433
434 if(alpha == 0){
435 return pChanges;
436 }
437 S1.setPhi(rad * G4UniformRand() * 2.0 * CLHEP::pi);
438 S1.setTheta(rad * std::acos(1.0 + 1./alpha * std::log(1.0 - G4UniformRand() * (1 - std::exp(-2.0 * alpha)))));
439
440 r1 = (D1 * S1 + D2 * S2) / (D1 + D2);
441 r2 = D2 * (S2 - S1) / (D1 + D2);
442 }
443 }
444
445 auto pTrackA = const_cast<G4Track*>(pChanges->GetTrackA());
446 auto pTrackB = const_cast<G4Track*>(pChanges->GetTrackB());
447
448 pTrackA->SetPosition(r1);
449 pTrackB->SetPosition(r2);
450
451 pTrackA->SetGlobalTime(globalTime);
452 pTrackB->SetGlobalTime(globalTime);
453
454 pTrackA->SetTrackStatus(fStopButAlive);
455 pTrackB->SetTrackStatus(fStopButAlive);
456
457 const G4int nbProducts = pReactionData->GetNbProducts();
458
459 if(nbProducts){
460
461 const G4double sqrD1 = D1 == 0. ? 0. : std::sqrt(D1);
462 const G4double sqrD2 = D2 == 0. ? 0. : std::sqrt(D2);
463 if((sqrD1 + sqrD2) == 0){
464 return pChanges;
465 }
466 const G4double inv_numerator = 1./(sqrD1 + sqrD2);
467 const G4ThreeVector reactionSite = sqrD2 * inv_numerator * trackA.GetPosition()
468 + sqrD1 * inv_numerator * trackB.GetPosition();
469
470 std::vector<G4ThreeVector> position;
471
472 if(nbProducts == 1){
473 position.push_back(reactionSite);
474 }else if(nbProducts == 2){
475 position.push_back(trackA.GetPosition());
476 position.push_back(trackB.GetPosition());
477 }else if (nbProducts == 3){
478 position.push_back(reactionSite);
479 position.push_back(trackA.GetPosition());
480 position.push_back(trackB.GetPosition());
481 }
482
483 for(G4int u = 0; u < nbProducts; u++){
484
485 auto product = new G4Molecule(pReactionData->GetProduct(u));
486 auto productTrack = product->BuildTrack(globalTime,
487 position[u]);
488
489 productTrack->SetTrackStatus(fAlive);
490 fTrackHolder->Push(productTrack);
491
492 pChanges->AddSecondary(productTrack);
493
494 G4int I = FindBin(fNx, fXMin, fXMax, position[u].x());
495 G4int J = FindBin(fNy, fYMin, fYMax, position[u].y());
496 G4int K = FindBin(fNz, fZMin, fZMax, position[u].z());
497
498 spaceBinned[I][J][K].push_back(productTrack);
499
500 Sampling(productTrack);
501 }
502 }
503
504 fTrackHolder->MergeSecondariesWithMainList();
505 pChanges->KillParents(true);
506 return pChanges;
507}
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:74
CLHEP::Hep3Vector G4ThreeVector
@ fAlive
@ fStopButAlive
void setTheta(double)
double mag() const
void setMag(double)
Definition: ThreeVector.cc:20
void setPhi(double)
virtual void Push(G4Track *)
void MergeSecondariesWithMainList()
const G4MolecularConfiguration * GetMolecularConfiguration() const
Definition: G4Molecule.cc:530
G4double GetGlobalTime() const
Definition: G4Scheduler.hh:364
void SetPosition(const G4ThreeVector &aValue)
const G4ThreeVector & GetPosition() const
G4double GetGlobalTime() const

Referenced by FindReaction().

◆ operator=()

G4DNAIRT & G4DNAIRT::operator= ( const G4DNAIRT other)
delete

◆ SamplePDC()

G4double G4DNAIRT::SamplePDC ( G4double  a,
G4double  b 
)

Definition at line 360 of file G4DNAIRT.cc.

360 {
361
362 G4double p = 2.0 * std::sqrt(2.0*b/a);
363 G4double q = 2.0 / std::sqrt(2.0*b/a);
364 G4double M = max(1.0/(a*a),3.0*b/a);
365
366 G4double X, U, lambdax;
367
368 G4int ntrials = 0;
369 while(1) {
370
371 // Generate X
372 U = G4UniformRand();
373 if ( U < p/(p + q * M) ) X = pow(U * (p + q * M) / 2, 2);
374 else X = pow(2/((1-U)*(p+q*M)/M),2);
375
376 U = G4UniformRand();
377
378 lambdax = std::exp(-b*b/X) * ( 1.0 - a * std::sqrt(CLHEP::pi * X) * erfc->erfcx(b/std::sqrt(X) + a*std::sqrt(X)));
379
380 if ((X <= 2.0*b/a && U <= lambdax) ||
381 (X >= 2.0*b/a && U*M/X <= lambdax)) break;
382
383 ntrials++;
384
385 if ( ntrials > 10000 ){
386 G4cout<<"Totally rejected"<<'\n';
387 return -1.0;
388 }
389 }
390 return X;
391}
#define M(row, col)
G4GLOB_DLL std::ostream G4cout
static G4double erfcx(G4double x)
T max(const T t1, const T t2)
brief Return the largest of the two arguments

Referenced by GetIndependentReactionTime().

◆ Sampling()

void G4DNAIRT::Sampling ( G4Track track)

Definition at line 176 of file G4DNAIRT.cc.

176 {
177 G4Molecule* molA = G4Molecule::GetMolecule(track);
178 const G4MolecularConfiguration* molConfA = molA->GetMolecularConfiguration();
179 if(molConfA->GetDiffusionCoefficient() == 0) return;
180
181 const vector<const G4MolecularConfiguration*>* reactivesVector =
183
184 if(reactivesVector == nullptr) return;
185
187 G4double minTime = timeMax;
188
189 xiniIndex = FindBin(fNx, fXMin, fXMax, track->GetPosition().x()-fRCutOff);
190 xendIndex = FindBin(fNx, fXMin, fXMax, track->GetPosition().x()+fRCutOff);
191 yiniIndex = FindBin(fNy, fYMin, fYMax, track->GetPosition().y()-fRCutOff);
192 yendIndex = FindBin(fNy, fYMin, fYMax, track->GetPosition().y()+fRCutOff);
193 ziniIndex = FindBin(fNz, fZMin, fZMax, track->GetPosition().z()-fRCutOff);
194 zendIndex = FindBin(fNz, fZMin, fZMax, track->GetPosition().z()+fRCutOff);
195
196 for ( int ii = xiniIndex; ii <= xendIndex; ii++ ) {
197 for ( int jj = yiniIndex; jj <= yendIndex; jj++ ) {
198 for ( int kk = ziniIndex; kk <= zendIndex; kk++ ) {
199
200 std::vector<G4Track*> spaceBin = spaceBinned[ii][jj][kk];
201 for ( int n = 0; n < (int)spaceBinned[ii][jj][kk].size(); n++ ) {
202 if(!spaceBin[n] || track == spaceBin[n]) continue;
203 if(spaceBin[n]->GetTrackStatus() == fStopButAlive) continue;
204
205 G4Molecule* molB = G4Molecule::GetMolecule(spaceBin[n]);
206 if(!molB) continue;
207
208 const G4MolecularConfiguration* molConfB = molB->GetMolecularConfiguration();
209 if(molConfB->GetDiffusionCoefficient() == 0) continue;
210
211 auto it = std::find(reactivesVector->begin(), reactivesVector->end(), molConfB);
212 if(it == reactivesVector->end()) continue;
213
214 G4ThreeVector orgPosB = spaceBin[n]->GetPosition();
215 G4double dt = track->GetGlobalTime() - spaceBin[n]->GetGlobalTime();
216 G4ThreeVector newPosB = orgPosB;
217
218 if(dt > 0){
219 G4double sigma, x, y, z;
220 G4double diffusionCoefficient = G4Molecule::GetMolecule(spaceBin[n])->GetDiffusionCoefficient();
221
222 sigma = std::sqrt(2.0 * diffusionCoefficient * dt);
223
224 x = G4RandGauss::shoot(0., 1.0)*sigma;
225 y = G4RandGauss::shoot(0., 1.0)*sigma;
226 z = G4RandGauss::shoot(0., 1.0)*sigma;
227
228 newPosB = orgPosB + G4ThreeVector(x,y,z);
229 }else if(dt < 0) continue;
230
231 G4double r0 = (newPosB - track->GetPosition()).mag();
233 molConfB,
234 r0);
235 if(irt>=0 && irt<timeMax - globalTime)
236 {
237 irt += globalTime;
238 if(irt < minTime) minTime = irt;
239#ifdef DEBUG
240 G4cout<<irt<<'\t'<<molConfA->GetName()<<" "<<track->GetTrackID()<<'\t'<<molConfB->GetName()<<" "<<spaceBin[n]->GetTrackID()<<'\n';
241#endif
242 fReactionSet->AddReaction(irt,track,spaceBin[n]);
243 }
244 }
245 spaceBin.clear();
246 }
247 }
248 }
249
250// Scavenging & first order reactions
251
252 auto fReactionDatas = fMolReactionTable->GetReactionData(molConfA);
253 G4double index = -1;
254 //change the scavenging filter of the IRT beyond 1 us proposed by Naoki and Jose
255 if(timeMax > 1*us)
256 {
257 minTime = timeMax;
258 }
259 //
260
261 for(size_t u=0; u<fReactionDatas->size();u++){
262 if((*fReactionDatas)[u]->GetReactant2()->GetDiffusionCoefficient() == 0){
263 G4double kObs = (*fReactionDatas)[u]->GetObservedReactionRateConstant();
264 if(kObs == 0) continue;
265 G4double time = -(std::log(1.0 - G4UniformRand())/kObs) + globalTime;
266 if( time < minTime && time >= globalTime && time < timeMax){
267 minTime = time;
268 index = (G4int)u;
269 }
270 }
271 }
272
273 if(index != -1){
274#ifdef DEBUG
275 G4cout<<"scavenged: "<<minTime<<'\t'<<molConfA->GetName()<<it_begin->GetTrackID()<<'\n';
276#endif
277 G4Molecule* fakeMol = new G4Molecule((*fReactionDatas)[index]->GetReactant2());
278 G4Track* fakeTrack = fakeMol->BuildTrack(globalTime,track->GetPosition());
279 fTrackHolder->Push(fakeTrack);
280 fReactionSet->AddReaction(minTime, track, fakeTrack);
281 }
282}
double z() const
double x() const
double y() const
G4double GetIndependentReactionTime(const G4MolecularConfiguration *, const G4MolecularConfiguration *, G4double)
Definition: G4DNAIRT.cc:285
const ReactantList * CanReactWith(Reactant *) const
void AddReaction(G4double time, G4Track *trackA, G4Track *trackB)
const G4String & GetName() const
static G4Molecule * GetMolecule(const G4Track *)
Definition: G4Molecule.cc:88
G4Track * BuildTrack(G4double globalTime, const G4ThreeVector &Position)
Definition: G4Molecule.cc:371
G4double GetDiffusionCoefficient() const
Definition: G4Molecule.cc:514
G4int GetTrackID() const

Referenced by IRTSampling(), and MakeReaction().

◆ SetReactionModel()

void G4DNAIRT::SetReactionModel ( G4VDNAReactionModel model)

Definition at line 559 of file G4DNAIRT.cc.

560{
561 fpReactionModel = model;
562}

◆ SpaceBinning()

void G4DNAIRT::SpaceBinning ( )

Definition at line 138 of file G4DNAIRT.cc.

138 {
139 auto it_begin = fTrackHolder->GetMainList()->begin();
140 while(it_begin != fTrackHolder->GetMainList()->end()){
141
142 G4ThreeVector position = it_begin->GetPosition();
143
144 if ( fXMin > position.x() ) fXMin = position.x();
145 if ( fYMin > position.y() ) fYMin = position.y();
146 if ( fZMin > position.z() ) fZMin = position.z();
147
148 if ( fXMax < position.x() ) fXMax = position.x();
149 if ( fYMax < position.y() ) fYMax = position.y();
150 if ( fZMax < position.z() ) fZMax = position.z();
151
152 ++it_begin;
153 }
154
155 fNx = G4int((fXMax-fXMin)/fRCutOff) == 0 ? 1 : G4int((fXMax-fXMin)/fRCutOff);
156 fNy = G4int((fYMax-fYMin)/fRCutOff) == 0 ? 1 : G4int((fYMax-fYMin)/fRCutOff);
157 fNz = G4int((fZMax-fZMin)/fRCutOff) == 0 ? 1 : G4int((fZMax-fZMin)/fRCutOff);
158
159}

Referenced by Initialize().

◆ TestReactibility()

G4bool G4DNAIRT::TestReactibility ( const G4Track ,
const G4Track ,
G4double  ,
G4bool   
)
overridevirtual

Implements G4VITReactionProcess.

Definition at line 551 of file G4DNAIRT.cc.

555{
556 return true;
557}

Member Data Documentation

◆ fMolReactionTable

const G4DNAMolecularReactionTable*& G4DNAIRT::fMolReactionTable
protected

Definition at line 93 of file G4DNAIRT.hh.

Referenced by GetIndependentReactionTime(), MakeReaction(), and Sampling().

◆ fpReactionModel

G4VDNAReactionModel* G4DNAIRT::fpReactionModel
protected

Definition at line 94 of file G4DNAIRT.hh.

Referenced by G4DNAIRT(), and SetReactionModel().


The documentation for this class was generated from the following files: