| Geant4 Cross Reference |
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 //
27 #include "G4DNASmoluchowskiReactionModel.hh"
28 #include "Randomize.hh"
29 #include "G4Track.hh"
30 #include "G4DNAMolecularReactionTable.hh"
31 #include "G4UnitsTable.hh"
32 #include "G4Molecule.hh"
33 #include "G4Exp.hh"
34
35 G4DNASmoluchowskiReactionModel::G4DNASmoluchowskiReactionModel() = default;
36
37 G4DNASmoluchowskiReactionModel::~G4DNASmoluchowskiReactionModel() = default;
38
39 void G4DNASmoluchowskiReactionModel::Initialise(const G4MolecularConfiguration* pMolecule,
40 const G4Track&)
41 {
42 fpReactionData = fpReactionTable->GetReactionData(pMolecule);
43 }
44
45 void G4DNASmoluchowskiReactionModel::InitialiseToPrint(const G4MolecularConfiguration* pMolecule)
46 {
47 fpReactionData = fpReactionTable->GetReactionData(pMolecule);
48 }
49
50 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4MolecularConfiguration* pMol1,
51 const G4MolecularConfiguration* pMol2)
52 {
53 G4double __output = fpReactionTable->GetReactionData(pMol1, pMol2)->GetEffectiveReactionRadius();
54 return __output;
55 }
56
57 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4int& __i)
58 {
59 G4double __output = (*fpReactionData)[__i]->GetEffectiveReactionRadius();
60 return __output;
61 }
62
63 G4bool G4DNASmoluchowskiReactionModel::FindReaction(const G4Track& __trackA,
64 const G4Track& __trackB,
65 const G4double __reactionRadius,
66 G4double& __separationDistance,
67 const G4bool __alongStepReaction)
68 {
69 const G4double R2 = __reactionRadius * __reactionRadius;
70 G4double postStepSeparation = 0.;
71 bool do_break = false;
72 int k = 0;
73
74 for (; k < 3; ++k)
75 {
76 postStepSeparation += std::pow(
77 __trackA.GetPosition()[k] - __trackB.GetPosition()[k], 2);
78
79 if (postStepSeparation > R2)
80 {
81 do_break = true;
82 break;
83 }
84 }
85
86 if (!do_break)
87 {
88 // The loop was not break
89 // => r^2 < R^2
90 __separationDistance = std::sqrt(postStepSeparation);
91 return true;
92 }
93 if (__alongStepReaction)
94 {
95 //Along step check and the loop has break
96
97 // Continue loop
98 for (; k < 3; ++k)
99 {
100 postStepSeparation += std::pow(
101 __trackA.GetPosition()[k] - __trackB.GetPosition()[k], 2);
102 }
103 // Use Green approach : the Brownian bridge
104 __separationDistance = (postStepSeparation = std::sqrt(postStepSeparation));
105
106 auto pMoleculeA = GetMolecule(__trackA);
107 auto pMoleculeB = GetMolecule(__trackB);
108
109 G4double D = pMoleculeA->GetDiffusionCoefficient()
110 + pMoleculeB->GetDiffusionCoefficient();
111
112 const auto& preStepPositionA = __trackA.GetStep()->GetPreStepPoint()->GetPosition();
113 const auto& preStepPositionB = __trackB.GetStep()->GetPreStepPoint()->GetPosition();
114
115 G4double preStepSeparation = (preStepPositionA - preStepPositionB).mag();
116
117 //===================================
118 // Brownian bridge
119 G4double __probabiltyOfEncounter = G4Exp(-(preStepSeparation - __reactionRadius)
120 * (postStepSeparation - __reactionRadius)
121 / (D * (__trackB.GetStep()->GetDeltaTime()))
122 );
123 G4double __selectedPOE = G4UniformRand();
124
125 if (__selectedPOE <= __probabiltyOfEncounter)
126 {
127 return true;
128 }
129 //===================================
130 }
131
132 return false;
133 }
134