| 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 #include "IRTDamageReactionModel.hh"
27
28 #include "DNAGeometry.hh"
29 #include "DNAHit.hh"
30 #include "DetectorConstruction.hh"
31 #include "EventAction.hh"
32
33 #include "G4DNAMolecularMaterial.hh"
34 #include "G4DNAMolecularReactionTable.hh"
35 #include "G4ErrorFunction.hh"
36 #include "G4EventManager.hh"
37 #include "G4IRTUtils.hh"
38 #include "G4Molecule.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4RunManager.hh"
41 #include "G4Scheduler.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4UnitsTable.hh"
44 #include "Randomize.hh"
45
46 #include <vector>
47
48 //....oooOO0OOooo........oooOO0OOooo........oo
49
50 IRTDamageReactionModel::IRTDamageReactionModel
51 : G4VDNAHitModel(name), fMolecularReactionTa
52 {
53 auto det = dynamic_cast<const DetectorConstr
54 G4RunManager::GetRunManager()->GetUserDete
55 fpDNAGeometry = det->GetDNAGeometry();
56 }
57
58 //....oooOO0OOooo........oooOO0OOooo........oo
59
60 G4double IRTDamageReactionModel::GetTimeToEnco
61
62
63 {
64 G4double irt = -1;
65 const auto pMoleculeA = molA;
66 const auto pMoleculeB = molB;
67 auto reactionData = fMolecularReactionTable-
68 G4double D = molA->GetDiffusionCoefficient()
69 auto kobs = reactionData->GetObservedReactio
70 if (D == 0) {
71 G4ExceptionDescription exceptionDescriptio
72 exceptionDescription << "D = 0"
73 << " for : " << molA-
74 G4Exception(
75 "IRTDamageReactionModel"
76 "::GetTimeToEncounter()",
77 "MolecularIRTDamageReactionModel002", Fa
78 return -1;
79 }
80 G4double Reff = kobs / (4 * CLHEP::pi * D *
81
82 if (distance < Reff) {
83 return 0; //
84 }
85
86 G4double Winf = 0;
87
88 if (distance != 0) {
89 Winf = Reff / distance;
90 }
91 else {
92 G4ExceptionDescription exceptionDescriptio
93 exceptionDescription << "distance = " << d
94 << " Reff = " << Reff
95 << molB->GetName();
96 G4Exception(
97 "IRTDamageReactionModel"
98 "::GetTimeToEncounter()",
99 "MolecularIRTDamageReactionModel001", Fa
100 }
101
102 G4double U = G4UniformRand();
103
104 if (Winf != 0 && U < Winf) {
105 irt = (1.0 / (4 * D)) * std::pow((distance
106 }
107 return irt;
108 }
109
110 //....oooOO0OOooo........oooOO0OOooo........oo
111
112 void IRTDamageReactionModel::RecordDNADamage()
113 {
114 if (fpDNAPhyVolume == nullptr || fpTrack ==
115 G4ExceptionDescription exceptionDescriptio
116 exceptionDescription << "fpDNAPhyVolume ==
117 G4Exception(
118 "IRTDamageReactionModel"
119 "RecordDNA",
120 "NO_VOLUME001", FatalException, exceptio
121 }
122 const G4VTouchable* touchable = fpTrack->Get
123 if (touchable == nullptr) {
124 return;
125 }
126 auto pPhyVolum = const_cast<G4VPhysicalVolum
127 const G4MolecularConfiguration* radical = Ge
128
129 // particle position
130 const G4ThreeVector& pos_particle = fpTrack-
131 G4AffineTransform transform = touchable->Get
132 G4ThreeVector localpos_particle = transform.
133
134 // DNA position
135 G4ThreeVector localPos_DNA = pPhyVolum->GetT
136 G4ThreeVector globalPos_DNA =
137 touchable->GetHistory()->GetTopTransform()
138
139 const int64_t idx = fpDNAGeometry->GetGlobal
140
141 int64_t bp = fpDNAGeometry->GetBasePairFromU
142 G4int chainID = fpDNAGeometry->GetChainIDFro
143 G4int strandID = fpDNAGeometry->GetStrandIDF
144 molecule mol = fpDNAGeometry->GetMoleculeFro
145 G4int placementIdx = fpDNAGeometry->GetPlace
146
147 G4String chromosome =
148 fpDNAGeometry->GetChromosomeMapper()->GetC
149
150 auto dnaHit = new DNAHit(mol, placementIdx,
151 chromosome, radical
152
153 dynamic_cast<EventAction*>(G4EventManager::G
154 ->AddDNAHit(dnaHit);
155 }
156
157 //....oooOO0OOooo........oooOO0OOooo........oo
158
159 void IRTDamageReactionModel::MakeReaction(cons
160 {
161 // G4Track *pTrackA = const_cast<G4Track *>(
162 if (track.GetTrackStatus() == fStopAndKill)
163 G4ExceptionDescription exceptionDescriptio
164 exceptionDescription << "this track is kil
165 G4Exception(
166 "IRTDamageReactionModel"
167 "MakeReaction",
168 "NO_TRACK02", FatalException, exceptionD
169 }
170 if (G4Scheduler::Instance()->GetVerbose() !=
171 G4cout << "At time : " << std::setw(7)
172 << G4BestUnit(G4Scheduler::Instance
173 << " Reaction : " << GetIT(track)->
174 << fpDNAPhyVolume->GetName() << " -
175 G4cout << " Damaged " + fpDNAPhyVolume->Ge
176 G4cout << G4endl;
177 }
178 }
179
180 //....oooOO0OOooo........oooOO0OOooo........oo
181
182 /// DoReaction means : kill species and record
183 G4bool IRTDamageReactionModel::DoReaction(cons
184 cons
185 {
186 fReactionTime = reactionTime;
187
188 if (fReactionTime == G4Scheduler::Instance()
189 return false;
190 }
191
192 fpTrack = &track;
193 fpDNAPhyVolume = std::get<const G4VPhysicalV
194 MakeReaction(track);
195 RecordDNADamage();
196 return true;
197 }
198
199 //....oooOO0OOooo........oooOO0OOooo........oo
200
201 G4double IRTDamageReactionModel::CalculateReac
202 {
203 fpTrack = nullptr;
204 fminTimeStep = DBL_MAX;
205 fReactionTime = DBL_MAX;
206 fpDNAPhyVolume = nullptr;
207 if (fpDNAGeometry == nullptr) {
208 G4ExceptionDescription exceptionDescriptio
209 exceptionDescription << "no fpDNAGeometry"
210 G4Exception(
211 "IRTDamageReactionModel"
212 "::CalculateReactionTime()",
213 "MolecularIRTDamageReactionModel007", Fa
214 }
215
216 auto pMoleculeA = GetMolecule(track);
217 auto pMolConfA = pMoleculeA->GetMolecularCon
218 const auto pReactantList = fMolecularReactio
219
220 if (pReactantList == nullptr) {
221 G4ExceptionDescription exceptionDescriptio
222 exceptionDescription << "!!!!!!!!!!!!!!!!!
223 G4cout << "!!! WARNING" << G4endl;
224 G4cout << "IRTDamageReactionModel::Calcula
225 "for the reaction because the mo
226 << pMoleculeA->GetName() << " does
227 << G4endl;
228 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl
229 G4Exception(
230 "IRTDamageReactionModel"
231 "::CalculateReactionTime()",
232 "MolecularIRTDamageReactionModel003", Fa
233 return -1;
234 }
235
236 size_t nbReactives = pReactantList->size();
237
238 if (nbReactives == 0) {
239 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl
240 G4cout << "!!! WARNING" << G4endl;
241 G4cout << "IRTDamageReactionModel::Calcula
242 "for the reaction because the mo
243 << pMoleculeA->GetName() << " does
244 << "This message can also result fr
245 "reaction table."
246 << G4endl;
247 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl
248 return -1;
249 }
250 const G4VTouchable* touchable = track.GetTou
251 if (touchable == nullptr) {
252 return -1;
253 }
254
255 const G4LogicalVolume* logicalVolume = touch
256 const G4ThreeVector& globalPos_track = track
257 const G4ThreeVector& localPos =
258 touchable->GetHistory()->GetTopTransform()
259
260 G4double D = pMoleculeA->GetDiffusionCoeffic
261 std::vector<G4VPhysicalVolume*> result_pv;
262 result_pv.clear();
263 fpDNAGeometry->FindNearbyMolecules(logicalVo
264 G4IRTUtil
265
266 if (result_pv.empty()) {
267 return -1;
268 }
269 for (const auto& physicalVolume : result_pv)
270 const G4Material* material = physicalVolum
271
272 G4MolecularConfiguration* dna_molConf =
273 G4DNAMolecularMaterial::Instance()->GetM
274 auto it = std::find(pReactantList->begin()
275 if (it == pReactantList->end()) {
276 continue;
277 }
278
279 G4ThreeVector localPos_DNA = physicalVolum
280 G4ThreeVector globalPos_DNA =
281 touchable->GetHistory()->GetTopTransform
282 G4double distance = (localPos_DNA - localP
283
284 G4double distance2 = distance * distance;
285 auto reactionData =
286 G4DNAMolecularReactionTable::Instance()-
287
288 G4double kobs = reactionData->GetObservedR
289 G4double Reff = kobs / (4 * CLHEP::pi * D
290
291 if (distance2 < Reff * Reff) {
292 fminTimeStep = 0.;
293 vp = physicalVolume;
294 }
295 else {
296 G4double tempMinET = GetTimeToEncounter(
297 if (tempMinET > G4Scheduler::Instance()-
298 continue;
299 }
300 if (tempMinET >= fminTimeStep) {
301 continue;
302 }
303 fminTimeStep = tempMinET;
304 vp = physicalVolume;
305 }
306 }
307 if (fminTimeStep > G4Scheduler::Instance()->
308 && fminTimeStep < G4Scheduler::Instance(
309 {
310 fminTimeStep = G4Scheduler::Instance()->Ge
311 }
312 return fminTimeStep;
313 }
314 //....oooOO0OOooo........oooOO0OOooo........oo
315