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