| 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
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.hh"
38 #include "G4Scheduler.hh"
39 #include "G4UnitsTable.hh"
40
41 G4DNAMakeReaction::G4DNAMakeReaction()
42 :
43 fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
44 , fpReactionModel(nullptr)
45 , fpTimeStepper(nullptr)
46 , fTimeStep(0)
47 {
48 }
49
50 G4DNAMakeReaction::G4DNAMakeReaction(G4VDNAReactionModel* pReactionModel)
51 : G4DNAMakeReaction()
52 {
53 fpReactionModel = pReactionModel;
54 }
55
56 void G4DNAMakeReaction::SetTimeStepComputer(G4VITTimeStepComputer* pStepper)
57 {
58 fpTimeStepper = pStepper;
59 }
60
61 G4bool G4DNAMakeReaction::TestReactibility(const G4Track& /*trackA*/,
62 const G4Track& /*trackB*/,
63 G4double currentStepTime,
64 G4bool /*userStepTimeLimit*/) /*const*/
65 {
66 fTimeStep = currentStepTime;
67 return true;
68 }
69
70 std::unique_ptr<G4ITReactionChange>
71 G4DNAMakeReaction::MakeReaction(const G4Track &trackA,
72 const G4Track &trackB)
73 {
74 auto & tA = const_cast<G4Track&>(trackA);
75 auto & tB = const_cast<G4Track&>(trackB);
76 UpdatePositionForReaction( tA , tB );//TODO: should change it
77
78 std::unique_ptr<G4ITReactionChange> pChanges(new G4ITReactionChange());
79 pChanges->Initialize(trackA, trackB);
80
81 const auto pMoleculeA = GetMolecule(trackA)->GetMolecularConfiguration();
82 const auto pMoleculeB = GetMolecule(trackB)->GetMolecularConfiguration();
83
84 const auto pReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
85 const G4int nbProducts = pReactionData->GetNbProducts();
86 if (nbProducts != 0)
87 {
88 const G4double D1 = pMoleculeA->GetDiffusionCoefficient();
89 const G4double D2 = pMoleculeB->GetDiffusionCoefficient();
90 const G4double sqrD1 = D1 == 0. ? 0. : std::sqrt(D1);
91 const G4double sqrD2 = D2 == 0. ? 0. : std::sqrt(D2);
92 const G4double inv_numerator = 1./(sqrD1 + sqrD2);
93 const G4ThreeVector reactionSite = sqrD2 * inv_numerator * tA.GetPosition()
94 + sqrD1 * inv_numerator * tB.GetPosition();
95
96 G4double u = G4UniformRand();
97 auto randP = (1-u) * tA.GetPosition() + u * tB.GetPosition();
98
99 for (G4int j = 0; j < nbProducts; ++j)
100 {
101 auto pProduct = new G4Molecule(pReactionData->GetProduct(j));
102 auto pProductTrack = pProduct->BuildTrack(trackA.GetGlobalTime(), (reactionSite + randP)/2);
103 pProductTrack->SetTrackStatus(fAlive);
104 G4ITTrackHolder::Instance()->Push(pProductTrack);
105 pChanges->AddSecondary(pProductTrack);
106 }
107 }
108 pChanges->KillParents(true);
109 return pChanges;
110 }
111
112 void G4DNAMakeReaction::SetReactionModel(G4VDNAReactionModel* pReactionModel)
113 {
114 fpReactionModel = pReactionModel;
115 }
116
117 void G4DNAMakeReaction::UpdatePositionForReaction(G4Track& trackA,
118 G4Track& trackB)
119 {
120 const auto pMoleculeA = GetMolecule(trackA)->GetMolecularConfiguration();
121 const auto pMoleculeB = GetMolecule(trackB)->GetMolecularConfiguration();
122 G4double D1 = pMoleculeA->GetDiffusionCoefficient();
123 G4double D2 = pMoleculeB->GetDiffusionCoefficient();
124
125 G4double reactionRadius = fpReactionModel->GetReactionRadius( pMoleculeA, pMoleculeB );
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 exceptionDescription;
146 exceptionDescription << "Two particles are overlap: "
147 <<GetMolecule(trackA)->GetName()
148 <<" and "<<GetMolecule(trackB)->GetName()
149 <<" at "<<trackA.GetPosition();
150 G4Exception("G4DNAMakeReaction::PrepareForReaction()",
151 "G4DNAMakeReaction003",
152 FatalErrorInArgument,exceptionDescription);
153 }
154 S1.setMag(reactionRadius);
155
156 const G4double dt = fTimeStep;//irt - actualize molecule time
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 ) / s22;
163 G4double alpha = reactionRadius * distance / (2 * (D1 + D2) * dt );
164
165 G4ThreeVector S2 = (p1 + ( s12 / s22 ) * p2) +
166 G4ThreeVector(G4RandGauss::shoot(0.0, sigma),
167 G4RandGauss::shoot(0.0, sigma),
168 G4RandGauss::shoot(0.0, sigma));
169
170 S1.setPhi(rad * G4UniformRand() * 2.0 * CLHEP::pi);
171
172 S1.setTheta(rad * std::acos( 1.0 + (1. / alpha) *
173 std::log(1.0 - G4UniformRand() *
174 (1.-std::exp(-2.0 * alpha)))));
175
176 const G4ThreeVector R1 = (D1 * S1 + D2 * S2) / (D1 + D2);
177 const G4ThreeVector R2 = D2 * (S2 - S1) / (D1 + D2);
178
179 trackA.SetPosition(R1);
180 trackB.SetPosition(R2);
181 }
182 }
183
184 std::vector<std::unique_ptr<G4ITReactionChange>>
185 G4DNAMakeReaction::FindReaction(G4ITReactionSet* pReactionSet,
186 const G4double currentStepTime,
187 const G4double /*globalTime*/,
188 const G4bool /*reachedUserStepTimeLimit*/)
189 {
190 std::vector<std::unique_ptr<G4ITReactionChange>> ReactionInfo;
191 ReactionInfo.clear();
192 auto stepper = dynamic_cast<G4DNAIndependentReactionTimeStepper*>(fpTimeStepper);
193 if(stepper == nullptr){
194 return ReactionInfo;
195 }else
196 {
197 do{
198 auto pReactionChange = stepper->
199 FindReaction(pReactionSet,currentStepTime);
200 if (pReactionChange != nullptr)
201 {
202 ReactionInfo.push_back(std::move(pReactionChange));
203 }
204 }while (!pReactionSet->GetReactionsPerTime().empty());
205 }
206 return ReactionInfo;
207 }
208