| 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 // $Id: G4DNASmoluchowskiReactionModel.cc 65695 2012-11-27 11:39:12Z gunter $
26 // 27 //
27 #include "G4DNASmoluchowskiReactionModel.hh" 28 #include "G4DNASmoluchowskiReactionModel.hh"
28 #include "Randomize.hh" 29 #include "Randomize.hh"
29 #include "G4Track.hh" 30 #include "G4Track.hh"
30 #include "G4DNAMolecularReactionTable.hh" 31 #include "G4DNAMolecularReactionTable.hh"
31 #include "G4UnitsTable.hh" <<
32 #include "G4Molecule.hh" <<
33 #include "G4Exp.hh" <<
34 32
35 G4DNASmoluchowskiReactionModel::G4DNASmoluchow << 33 G4DNASmoluchowskiReactionModel::G4DNASmoluchowskiReactionModel() : G4VDNAReactionModel()
>> 34 {
>> 35 fReactionData = 0 ;
>> 36 }
>> 37
>> 38 G4DNASmoluchowskiReactionModel::G4DNASmoluchowskiReactionModel(const G4DNASmoluchowskiReactionModel& __right) :
>> 39 G4VDNAReactionModel(__right)
>> 40 {
>> 41 fReactionData = 0 ;
>> 42 }
>> 43
>> 44 G4DNASmoluchowskiReactionModel& G4DNASmoluchowskiReactionModel::operator=(const G4DNASmoluchowskiReactionModel& right)
>> 45 {
>> 46 if(this == &right) return *this;
>> 47 fReactionData = 0;
>> 48 return *this;
>> 49 }
36 50
37 G4DNASmoluchowskiReactionModel::~G4DNASmolucho << 51 G4DNASmoluchowskiReactionModel::~G4DNASmoluchowskiReactionModel()
>> 52 {
>> 53 fReactionData = 0 ;
>> 54 }
38 55
39 void G4DNASmoluchowskiReactionModel::Initialis << 56 void G4DNASmoluchowskiReactionModel::Initialise(const G4Molecule* __molecule, const G4Track&)
40 <<
41 { 57 {
42 fpReactionData = fpReactionTable->GetReact << 58 fReactionData = fReactionTable->GetReactionData(__molecule);
43 } 59 }
44 60
45 void G4DNASmoluchowskiReactionModel::Initialis << 61 void G4DNASmoluchowskiReactionModel::InitialiseToPrint(const G4Molecule* __molecule)
46 { 62 {
47 fpReactionData = fpReactionTable->GetReact << 63 fReactionData = fReactionTable->GetReactionData(__molecule);
48 } 64 }
49 65
50 G4double G4DNASmoluchowskiReactionModel::GetRe << 66 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4Molecule* mol1,
51 << 67 const G4Molecule* mol2)
52 { 68 {
53 G4double __output = fpReactionTable->GetRe << 69 G4double __output = fReactionTable -> GetReactionData(mol1,mol2)->GetReducedReactionRadius();
54 return __output; << 70 return __output ;
55 } 71 }
56 72
57 G4double G4DNASmoluchowskiReactionModel::GetRe << 73 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4int __i)
58 { 74 {
59 G4double __output = (*fpReactionData)[__i] << 75 G4double __output = (*fReactionData)[__i] -> GetReducedReactionRadius();
60 return __output; << 76 return __output ;
61 } 77 }
62 78
63 G4bool G4DNASmoluchowskiReactionModel::FindRea 79 G4bool G4DNASmoluchowskiReactionModel::FindReaction(const G4Track& __trackA,
64 << 80 const G4Track& __trackB,
65 << 81 const G4double __R,
66 << 82 G4double& __r,
67 << 83 const G4bool __alongStepReaction)
68 { << 84 {
69 const G4double R2 = __reactionRadius * __r << 85 G4double __postStepSeparation = 0;
70 G4double postStepSeparation = 0.; << 86 bool __do_break = false ;
71 bool do_break = false; << 87 G4double __R2 = __R*__R ;
72 int k = 0; << 88 int k = 0 ;
73 89
74 for (; k < 3; ++k) << 90 for(; k < 3 ; k++)
75 { 91 {
76 postStepSeparation += std::pow( << 92 __postStepSeparation += std::pow(__trackA.GetPosition()[k] - __trackB.GetPosition()[k],2);
77 __trackA.GetPosition()[k] - __trac <<
78 93
79 if (postStepSeparation > R2) << 94 if(__postStepSeparation > __R2)
80 { 95 {
81 do_break = true; << 96 __do_break = true ;
82 break; << 97 break ;
83 } 98 }
84 } 99 }
85 100
86 if (!do_break) << 101 if(__do_break == false)
87 { 102 {
88 // The loop was not break << 103 // The loop was not break
89 // => r^2 < R^2 << 104 // => __r^2 < __R^2
90 __separationDistance = std::sqrt(postS << 105 __r = std::sqrt(__postStepSeparation);
91 return true; 106 return true;
92 } 107 }
93 if (__alongStepReaction) << 108 else if(__alongStepReaction == true)
94 { 109 {
95 //Along step check and the loop has br << 110 //G4cout << "alongStepReaction==true" << G4endl;
>> 111 //Along step cheack and
>> 112 // the loop has break
96 113
97 // Continue loop 114 // Continue loop
98 for (; k < 3; ++k) << 115 for(; k < 3 ; k++)
99 { 116 {
100 postStepSeparation += std::pow( << 117 __postStepSeparation += std::pow(__trackA.GetPosition()[k] - __trackB.GetPosition()[k],2);
101 __trackA.GetPosition()[k] - __ <<
102 } 118 }
103 // Use Green approach : the Brownian b 119 // Use Green approach : the Brownian bridge
104 __separationDistance = (postStepSepara << 120 __r = (__postStepSeparation = std::sqrt(__postStepSeparation) );
>> 121
>> 122 G4Molecule* __moleculeA = GetMolecule(__trackA);
>> 123 G4Molecule* __moleculeB = GetMolecule(__trackB);
105 124
106 auto pMoleculeA = GetMolecule(__trackA << 125 G4double __D = __moleculeA->GetDiffusionCoefficient() + __moleculeB->GetDiffusionCoefficient();
107 auto pMoleculeB = GetMolecule(__trackB <<
108 126
109 G4double D = pMoleculeA->GetDiffusionC << 127 G4ThreeVector __preStepPositionA = __trackA.GetStep()->GetPreStepPoint() ->GetPosition();
110 + pMoleculeB->GetDiffusio << 128 G4ThreeVector __preStepPositionB = __trackB.GetStep()->GetPreStepPoint() ->GetPosition();
111 129
112 const auto& preStepPositionA = __track << 130 if(__preStepPositionA == __trackA.GetPosition())
113 const auto& preStepPositionB = __track << 131 {
>> 132 G4ExceptionDescription exceptionDescription ;
>> 133 exceptionDescription << "The molecule : " << __moleculeA->GetName();
>> 134 exceptionDescription << " did not move since the previous step ";
>> 135 G4Exception("G4DNASmoluchowskiReactionModel::FindReaction","G4DNASmoluchowskiReactionModel",
>> 136 FatalErrorInArgument,exceptionDescription);
>> 137 }
114 138
115 G4double preStepSeparation = (preStepP << 139 G4double __preStepSeparation = (__preStepPositionA - __preStepPositionB).mag();
116 140
117 //=================================== << 141 G4double __probabiltyOfEncounter = std::exp(-(__preStepSeparation - __R)*(__postStepSeparation - __R)
118 // Brownian bridge << 142 / (__D* (__trackB.GetStep()->GetDeltaTime())));
119 G4double __probabiltyOfEncounter = G4E <<
120 <<
121 <<
122 ); <<
123 G4double __selectedPOE = G4UniformRand 143 G4double __selectedPOE = G4UniformRand();
124 144
125 if (__selectedPOE <= __probabiltyOfEnc << 145 if(__selectedPOE<=__probabiltyOfEncounter) return true;
126 { <<
127 return true; <<
128 } <<
129 //=================================== <<
130 } 146 }
131 147
132 return false; << 148 return false ;
133 } 149 }
134 150