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