| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26
27
28 #include "G4DNAMakeReaction.hh"
29 #include "G4DNAMolecularReactionTable.hh"
30 #include "G4VDNAReactionModel.hh"
31 #include "G4Molecule.hh"
32 #include "G4MoleculeFinder.hh"
33 #include "G4ITReactionChange.hh"
34 #include "Randomize.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4ITReaction.hh"
37 #include "G4DNAIndependentReactionTimeStepper.
38 #include "G4Scheduler.hh"
39 #include "G4UnitsTable.hh"
40
41 G4DNAMakeReaction::G4DNAMakeReaction()
42 :
43 fMolReactionTable(reference_cast<const G4
44 , fpReactionModel(nullptr)
45 , fpTimeStepper(nullptr)
46 , fTimeStep(0)
47 {
48 }
49
50 G4DNAMakeReaction::G4DNAMakeReaction(G4VDNARea
51 : G4DNAMakeReaction()
52 {
53 fpReactionModel = pReactionModel;
54 }
55
56 void G4DNAMakeReaction::SetTimeStepComputer(G4
57 {
58 fpTimeStepper = pStepper;
59 }
60
61 G4bool G4DNAMakeReaction::TestReactibility(con
62 con
63 G4d
64 G4b
65 {
66 fTimeStep = currentStepTime;
67 return true;
68 }
69
70 std::unique_ptr<G4ITReactionChange>
71 G4DNAMakeReaction::MakeReaction(const G4Track
72 const G4Track
73 {
74 auto & tA = const_cast<G4Track&>(trackA);
75 auto & tB = const_cast<G4Track&>(trackB);
76 UpdatePositionForReaction( tA , tB );//TOD
77
78 std::unique_ptr<G4ITReactionChange> pChang
79 pChanges->Initialize(trackA, trackB);
80
81 const auto pMoleculeA = GetMolecule(trackA
82 const auto pMoleculeB = GetMolecule(trackB
83
84 const auto pReactionData = fMolReactionTab
85 const G4int nbProducts = pReactionData->Ge
86 if (nbProducts != 0)
87 {
88 const G4double D1 = pMoleculeA->GetDif
89 const G4double D2 = pMoleculeB->GetDif
90 const G4double sqrD1 = D1 == 0. ? 0. :
91 const G4double sqrD2 = D2 == 0. ? 0. :
92 const G4double inv_numerator = 1./(sqr
93 const G4ThreeVector reactionSite = sqr
94 + sqr
95
96 G4double u = G4UniformRand();
97 auto randP = (1-u) * tA.GetPosition()
98
99 for (G4int j = 0; j < nbProducts; ++j)
100 {
101 auto pProduct = new G4Molecule(pRe
102 auto pProductTrack = pProduct->Bui
103 pProductTrack->SetTrackStatus(fAli
104 G4ITTrackHolder::Instance()->Push(
105 pChanges->AddSecondary(pProductTra
106 }
107 }
108 pChanges->KillParents(true);
109 return pChanges;
110 }
111
112 void G4DNAMakeReaction::SetReactionModel(G4VDN
113 {
114 fpReactionModel = pReactionModel;
115 }
116
117 void G4DNAMakeReaction::UpdatePositionForReact
118 G4T
119 {
120 const auto pMoleculeA = GetMolecule(trackA
121 const auto pMoleculeB = GetMolecule(trackB
122 G4double D1 = pMoleculeA->GetDiffusionCoef
123 G4double D2 = pMoleculeB->GetDiffusionCoef
124
125 G4double reactionRadius = fpReactionModel-
126 G4ThreeVector p1 = trackA.GetPosition();
127 G4ThreeVector p2 = trackB.GetPosition();
128
129 G4ThreeVector S1 = p1 - p2;
130 G4double distance = S1.mag();
131
132 if(D1 == 0)
133 {
134 trackB.SetPosition(p1);
135 return;
136 }
137 if(D2 == 0)
138 {
139 trackA.SetPosition(p2);
140 return;
141 }
142
143 if(distance == 0)
144 {
145 G4ExceptionDescription exceptionDescri
146 exceptionDescription << "Two particles
147 <<GetMolecule(tra
148 <<" and "<<GetMol
149 <<" at "<<trackA.
150 G4Exception("G4DNAMakeReaction::Prepar
151 "G4DNAMakeReaction003",
152 FatalErrorInArgument,excep
153 }
154 S1.setMag(reactionRadius);
155
156 const G4double dt = fTimeStep;//irt - actu
157
158 if(dt > 0)// irt > 0
159 {
160 G4double s12 = 2.0 * D1 * dt;
161 G4double s22 = 2.0 * D2 * dt;
162 G4double sigma = s12 + ( s12 * s12 ) /
163 G4double alpha = reactionRadius * dist
164
165 G4ThreeVector S2 = (p1 + ( s12 / s22 )
166 G4ThreeVector(G4Ran
167 G4Ran
168 G4Ran
169
170 S1.setPhi(rad * G4UniformRand() * 2.0
171
172 S1.setTheta(rad * std::acos( 1.0 + (1.
173 std
174
175
176 const G4ThreeVector R1 = (D1 * S1 + D2
177 const G4ThreeVector R2 = D2 * (S2 - S1
178
179 trackA.SetPosition(R1);
180 trackB.SetPosition(R2);
181 }
182 }
183
184 std::vector<std::unique_ptr<G4ITReactionChange
185 G4DNAMakeReaction::FindReaction(G4ITReactionSe
186 const G4double
187 const G4double
188 const G4bool /
189 {
190 std::vector<std::unique_ptr<G4ITReactionCh
191 ReactionInfo.clear();
192 auto stepper = dynamic_cast<G4DNAIndepende
193 if(stepper == nullptr){
194 return ReactionInfo;
195 }else
196 {
197 do{
198 auto pReactionChange = stepper->
199 FindReactio
200 if (pReactionChange != nullptr)
201 {
202 ReactionInfo.push_back(std::move
203 }
204 }while (!pReactionSet->GetReactionsPer
205 }
206 return ReactionInfo;
207 }
208