| 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 #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........oooOO0OOooo........oooOO0OOooo......
49
50 IRTDamageReactionModel::IRTDamageReactionModel(const G4String& name)
51 : G4VDNAHitModel(name), fMolecularReactionTable(G4DNAMolecularReactionTable::Instance())
52 {
53 auto det = dynamic_cast<const DetectorConstruction*>(
54 G4RunManager::GetRunManager()->GetUserDetectorConstruction());
55 fpDNAGeometry = det->GetDNAGeometry();
56 }
57
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59
60 G4double IRTDamageReactionModel::GetTimeToEncounter(const G4MolecularConfiguration* molA,
61 const G4MolecularConfiguration* molB,
62 const G4double& distance)
63 {
64 G4double irt = -1;
65 const auto pMoleculeA = molA;
66 const auto pMoleculeB = molB;
67 auto reactionData = fMolecularReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
68 G4double D = molA->GetDiffusionCoefficient() + molB->GetDiffusionCoefficient();
69 auto kobs = reactionData->GetObservedReactionRateConstant();
70 if (D == 0) {
71 G4ExceptionDescription exceptionDescription;
72 exceptionDescription << "D = 0"
73 << " for : " << molA->GetName() << " and " << molB->GetName();
74 G4Exception(
75 "IRTDamageReactionModel"
76 "::GetTimeToEncounter()",
77 "MolecularIRTDamageReactionModel002", FatalException, exceptionDescription);
78 return -1;
79 }
80 G4double Reff = kobs / (4 * CLHEP::pi * D * CLHEP::Avogadro);
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 exceptionDescription;
93 exceptionDescription << "distance = " << distance << " is uncorrected with "
94 << " Reff = " << Reff << " for : " << molA->GetName() << " and "
95 << molB->GetName();
96 G4Exception(
97 "IRTDamageReactionModel"
98 "::GetTimeToEncounter()",
99 "MolecularIRTDamageReactionModel001", FatalException, exceptionDescription);
100 }
101
102 G4double U = G4UniformRand();
103
104 if (Winf != 0 && U < Winf) {
105 irt = (1.0 / (4 * D)) * std::pow((distance - Reff) / G4ErrorFunction::erfcInv(U / Winf), 2);
106 }
107 return irt;
108 }
109
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111
112 void IRTDamageReactionModel::RecordDNADamage() const // let's suppose that this is not so much
113 {
114 if (fpDNAPhyVolume == nullptr || fpTrack == nullptr) {
115 G4ExceptionDescription exceptionDescription;
116 exceptionDescription << "fpDNAPhyVolume == nullptr or fpTrack == nullptr";
117 G4Exception(
118 "IRTDamageReactionModel"
119 "RecordDNA",
120 "NO_VOLUME001", FatalException, exceptionDescription);
121 }
122 const G4VTouchable* touchable = fpTrack->GetTouchable();
123 if (touchable == nullptr) {
124 return;
125 }
126 auto pPhyVolum = const_cast<G4VPhysicalVolume*>(fpDNAPhyVolume);
127 const G4MolecularConfiguration* radical = GetMolecule(fpTrack)->GetMolecularConfiguration();
128
129 // particle position
130 const G4ThreeVector& pos_particle = fpTrack->GetPosition();
131 G4AffineTransform transform = touchable->GetHistory()->GetTopTransform();
132 G4ThreeVector localpos_particle = transform.TransformPoint(pos_particle);
133
134 // DNA position
135 G4ThreeVector localPos_DNA = pPhyVolum->GetTranslation();
136 G4ThreeVector globalPos_DNA =
137 touchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(localPos_DNA);
138
139 const int64_t idx = fpDNAGeometry->GetGlobalUniqueID(pPhyVolum, touchable);
140
141 int64_t bp = fpDNAGeometry->GetBasePairFromUniqueID(idx);
142 G4int chainID = fpDNAGeometry->GetChainIDFromUniqueID(idx);
143 G4int strandID = fpDNAGeometry->GetStrandIDFromUniqueID(idx);
144 molecule mol = fpDNAGeometry->GetMoleculeFromUniqueID(idx);
145 G4int placementIdx = fpDNAGeometry->GetPlacementIndexFromUniqueID(idx);
146
147 G4String chromosome =
148 fpDNAGeometry->GetChromosomeMapper()->GetCurrentChromosomeKey(globalPos_DNA);
149
150 auto dnaHit = new DNAHit(mol, placementIdx, chainID, strandID, bp, globalPos_DNA, localPos_DNA,
151 chromosome, radical);
152
153 dynamic_cast<EventAction*>(G4EventManager::GetEventManager()->GetUserEventAction())
154 ->AddDNAHit(dnaHit);
155 }
156
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
158
159 void IRTDamageReactionModel::MakeReaction(const G4Track& track)
160 {
161 // G4Track *pTrackA = const_cast<G4Track *>(&track);
162 if (track.GetTrackStatus() == fStopAndKill) {
163 G4ExceptionDescription exceptionDescription;
164 exceptionDescription << "this track is killed";
165 G4Exception(
166 "IRTDamageReactionModel"
167 "MakeReaction",
168 "NO_TRACK02", FatalException, exceptionDescription);
169 }
170 if (G4Scheduler::Instance()->GetVerbose() != 0) {
171 G4cout << "At time : " << std::setw(7)
172 << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
173 << " Reaction : " << GetIT(track)->GetName() << " (" << track.GetTrackID() << ") + "
174 << fpDNAPhyVolume->GetName() << " -> ";
175 G4cout << " Damaged " + fpDNAPhyVolume->GetName();
176 G4cout << G4endl;
177 }
178 }
179
180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
181
182 /// DoReaction means : kill species and record DNA damage
183 G4bool IRTDamageReactionModel::DoReaction(const G4Track& track, const G4double& reactionTime,
184 const DNANode& vp)
185 {
186 fReactionTime = reactionTime;
187
188 if (fReactionTime == G4Scheduler::Instance()->GetLimitingTimeStep()) {
189 return false;
190 }
191
192 fpTrack = &track;
193 fpDNAPhyVolume = std::get<const G4VPhysicalVolume*>(vp);
194 MakeReaction(track);
195 RecordDNADamage();
196 return true;
197 }
198
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
200
201 G4double IRTDamageReactionModel::CalculateReactionTime(const G4Track& track, DNANode& vp)
202 {
203 fpTrack = nullptr;
204 fminTimeStep = DBL_MAX;
205 fReactionTime = DBL_MAX;
206 fpDNAPhyVolume = nullptr;
207 if (fpDNAGeometry == nullptr) {
208 G4ExceptionDescription exceptionDescription;
209 exceptionDescription << "no fpDNAGeometry" << G4endl;
210 G4Exception(
211 "IRTDamageReactionModel"
212 "::CalculateReactionTime()",
213 "MolecularIRTDamageReactionModel007", FatalException, exceptionDescription);
214 }
215
216 auto pMoleculeA = GetMolecule(track);
217 auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
218 const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
219
220 if (pReactantList == nullptr) {
221 G4ExceptionDescription exceptionDescription;
222 exceptionDescription << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
223 G4cout << "!!! WARNING" << G4endl;
224 G4cout << "IRTDamageReactionModel::CalculateReactionTime will return infinity "
225 "for the reaction because the molecule "
226 << pMoleculeA->GetName() << " does not have any reactants given in the reaction table."
227 << G4endl;
228 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
229 G4Exception(
230 "IRTDamageReactionModel"
231 "::CalculateReactionTime()",
232 "MolecularIRTDamageReactionModel003", FatalException, exceptionDescription);
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::CalculateReactionTime will return infinity "
242 "for the reaction because the molecule "
243 << pMoleculeA->GetName() << " does not have any reactants given in the reaction table."
244 << "This message can also result from a wrong implementation of the "
245 "reaction table."
246 << G4endl;
247 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
248 return -1;
249 }
250 const G4VTouchable* touchable = track.GetTouchable();
251 if (touchable == nullptr) {
252 return -1;
253 }
254
255 const G4LogicalVolume* logicalVolume = touchable->GetVolume()->GetLogicalVolume();
256 const G4ThreeVector& globalPos_track = track.GetPosition();
257 const G4ThreeVector& localPos =
258 touchable->GetHistory()->GetTopTransform().TransformPoint(globalPos_track);
259
260 G4double D = pMoleculeA->GetDiffusionCoefficient();
261 std::vector<G4VPhysicalVolume*> result_pv;
262 result_pv.clear();
263 fpDNAGeometry->FindNearbyMolecules(logicalVolume, localPos, result_pv,
264 G4IRTUtils::GetRCutOff(100 * ns));
265
266 if (result_pv.empty()) {
267 return -1;
268 }
269 for (const auto& physicalVolume : result_pv) {
270 const G4Material* material = physicalVolume->GetLogicalVolume()->GetMaterial();
271
272 G4MolecularConfiguration* dna_molConf =
273 G4DNAMolecularMaterial::Instance()->GetMolecularConfiguration(material);
274 auto it = std::find(pReactantList->begin(), pReactantList->end(), dna_molConf);
275 if (it == pReactantList->end()) {
276 continue;
277 }
278
279 G4ThreeVector localPos_DNA = physicalVolume->GetTranslation();
280 G4ThreeVector globalPos_DNA =
281 touchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(localPos_DNA);
282 G4double distance = (localPos_DNA - localPos).mag();
283
284 G4double distance2 = distance * distance;
285 auto reactionData =
286 G4DNAMolecularReactionTable::Instance()->GetReactionData(pMolConfA, dna_molConf);
287
288 G4double kobs = reactionData->GetObservedReactionRateConstant();
289 G4double Reff = kobs / (4 * CLHEP::pi * D * CLHEP::Avogadro);
290
291 if (distance2 < Reff * Reff) {
292 fminTimeStep = 0.;
293 vp = physicalVolume;
294 }
295 else {
296 G4double tempMinET = GetTimeToEncounter(pMolConfA, dna_molConf, distance);
297 if (tempMinET > G4Scheduler::Instance()->GetEndTime() || tempMinET < 0) {
298 continue;
299 }
300 if (tempMinET >= fminTimeStep) {
301 continue;
302 }
303 fminTimeStep = tempMinET;
304 vp = physicalVolume;
305 }
306 }
307 if (fminTimeStep > G4Scheduler::Instance()->GetLimitingTimeStep()
308 && fminTimeStep < G4Scheduler::Instance()->GetEndTime())
309 {
310 fminTimeStep = G4Scheduler::Instance()->GetLimitingTimeStep();
311 }
312 return fminTimeStep;
313 }
314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
315