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