File indexing completed on 2025-02-23 09:22:01
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 #include "IRTDamageReactionModel.hh"
0027
0028 #include "DNAGeometry.hh"
0029 #include "DNAHit.hh"
0030 #include "DetectorConstruction.hh"
0031 #include "EventAction.hh"
0032
0033 #include "G4DNAMolecularMaterial.hh"
0034 #include "G4DNAMolecularReactionTable.hh"
0035 #include "G4ErrorFunction.hh"
0036 #include "G4EventManager.hh"
0037 #include "G4IRTUtils.hh"
0038 #include "G4Molecule.hh"
0039 #include "G4PhysicalConstants.hh"
0040 #include "G4RunManager.hh"
0041 #include "G4Scheduler.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4UnitsTable.hh"
0044 #include "Randomize.hh"
0045
0046 #include <vector>
0047
0048
0049
0050 IRTDamageReactionModel::IRTDamageReactionModel(const G4String& name)
0051 : G4VDNAHitModel(name), fMolecularReactionTable(G4DNAMolecularReactionTable::Instance())
0052 {
0053 auto det = dynamic_cast<const DetectorConstruction*>(
0054 G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0055 fpDNAGeometry = det->GetDNAGeometry();
0056 }
0057
0058
0059
0060 G4double IRTDamageReactionModel::GetTimeToEncounter(const G4MolecularConfiguration* molA,
0061 const G4MolecularConfiguration* molB,
0062 const G4double& distance)
0063 {
0064 G4double irt = -1;
0065 const auto pMoleculeA = molA;
0066 const auto pMoleculeB = molB;
0067 auto reactionData = fMolecularReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
0068 G4double D = molA->GetDiffusionCoefficient() + molB->GetDiffusionCoefficient();
0069 auto kobs = reactionData->GetObservedReactionRateConstant();
0070 if (D == 0) {
0071 G4ExceptionDescription exceptionDescription;
0072 exceptionDescription << "D = 0"
0073 << " for : " << molA->GetName() << " and " << molB->GetName();
0074 G4Exception(
0075 "IRTDamageReactionModel"
0076 "::GetTimeToEncounter()",
0077 "MolecularIRTDamageReactionModel002", FatalException, exceptionDescription);
0078 return -1;
0079 }
0080 G4double Reff = kobs / (4 * CLHEP::pi * D * CLHEP::Avogadro);
0081
0082 if (distance < Reff) {
0083 return 0;
0084 }
0085
0086 G4double Winf = 0;
0087
0088 if (distance != 0) {
0089 Winf = Reff / distance;
0090 }
0091 else {
0092 G4ExceptionDescription exceptionDescription;
0093 exceptionDescription << "distance = " << distance << " is uncorrected with "
0094 << " Reff = " << Reff << " for : " << molA->GetName() << " and "
0095 << molB->GetName();
0096 G4Exception(
0097 "IRTDamageReactionModel"
0098 "::GetTimeToEncounter()",
0099 "MolecularIRTDamageReactionModel001", FatalException, exceptionDescription);
0100 }
0101
0102 G4double U = G4UniformRand();
0103
0104 if (Winf != 0 && U < Winf) {
0105 irt = (1.0 / (4 * D)) * std::pow((distance - Reff) / G4ErrorFunction::erfcInv(U / Winf), 2);
0106 }
0107 return irt;
0108 }
0109
0110
0111
0112 void IRTDamageReactionModel::RecordDNADamage() const
0113 {
0114 if (fpDNAPhyVolume == nullptr || fpTrack == nullptr) {
0115 G4ExceptionDescription exceptionDescription;
0116 exceptionDescription << "fpDNAPhyVolume == nullptr or fpTrack == nullptr";
0117 G4Exception(
0118 "IRTDamageReactionModel"
0119 "RecordDNA",
0120 "NO_VOLUME001", FatalException, exceptionDescription);
0121 }
0122 const G4VTouchable* touchable = fpTrack->GetTouchable();
0123 if (touchable == nullptr) {
0124 return;
0125 }
0126 auto pPhyVolum = const_cast<G4VPhysicalVolume*>(fpDNAPhyVolume);
0127 const G4MolecularConfiguration* radical = GetMolecule(fpTrack)->GetMolecularConfiguration();
0128
0129
0130 const G4ThreeVector& pos_particle = fpTrack->GetPosition();
0131 G4AffineTransform transform = touchable->GetHistory()->GetTopTransform();
0132 G4ThreeVector localpos_particle = transform.TransformPoint(pos_particle);
0133
0134
0135 G4ThreeVector localPos_DNA = pPhyVolum->GetTranslation();
0136 G4ThreeVector globalPos_DNA =
0137 touchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(localPos_DNA);
0138
0139 const int64_t idx = fpDNAGeometry->GetGlobalUniqueID(pPhyVolum, touchable);
0140
0141 int64_t bp = fpDNAGeometry->GetBasePairFromUniqueID(idx);
0142 G4int chainID = fpDNAGeometry->GetChainIDFromUniqueID(idx);
0143 G4int strandID = fpDNAGeometry->GetStrandIDFromUniqueID(idx);
0144 molecule mol = fpDNAGeometry->GetMoleculeFromUniqueID(idx);
0145 G4int placementIdx = fpDNAGeometry->GetPlacementIndexFromUniqueID(idx);
0146
0147 G4String chromosome =
0148 fpDNAGeometry->GetChromosomeMapper()->GetCurrentChromosomeKey(globalPos_DNA);
0149
0150 auto dnaHit = new DNAHit(mol, placementIdx, chainID, strandID, bp, globalPos_DNA, localPos_DNA,
0151 chromosome, radical);
0152
0153 dynamic_cast<EventAction*>(G4EventManager::GetEventManager()->GetUserEventAction())
0154 ->AddDNAHit(dnaHit);
0155 }
0156
0157
0158
0159 void IRTDamageReactionModel::MakeReaction(const G4Track& track)
0160 {
0161
0162 if (track.GetTrackStatus() == fStopAndKill) {
0163 G4ExceptionDescription exceptionDescription;
0164 exceptionDescription << "this track is killed";
0165 G4Exception(
0166 "IRTDamageReactionModel"
0167 "MakeReaction",
0168 "NO_TRACK02", FatalException, exceptionDescription);
0169 }
0170 if (G4Scheduler::Instance()->GetVerbose() != 0) {
0171 G4cout << "At time : " << std::setw(7)
0172 << G4BestUnit(G4Scheduler::Instance()->GetGlobalTime(), "Time")
0173 << " Reaction : " << GetIT(track)->GetName() << " (" << track.GetTrackID() << ") + "
0174 << fpDNAPhyVolume->GetName() << " -> ";
0175 G4cout << " Damaged " + fpDNAPhyVolume->GetName();
0176 G4cout << G4endl;
0177 }
0178 }
0179
0180
0181
0182
0183 G4bool IRTDamageReactionModel::DoReaction(const G4Track& track, const G4double& reactionTime,
0184 const DNANode& vp)
0185 {
0186 fReactionTime = reactionTime;
0187
0188 if (fReactionTime == G4Scheduler::Instance()->GetLimitingTimeStep()) {
0189 return false;
0190 }
0191
0192 fpTrack = &track;
0193 fpDNAPhyVolume = std::get<const G4VPhysicalVolume*>(vp);
0194 MakeReaction(track);
0195 RecordDNADamage();
0196 return true;
0197 }
0198
0199
0200
0201 G4double IRTDamageReactionModel::CalculateReactionTime(const G4Track& track, DNANode& vp)
0202 {
0203 fpTrack = nullptr;
0204 fminTimeStep = DBL_MAX;
0205 fReactionTime = DBL_MAX;
0206 fpDNAPhyVolume = nullptr;
0207 if (fpDNAGeometry == nullptr) {
0208 G4ExceptionDescription exceptionDescription;
0209 exceptionDescription << "no fpDNAGeometry" << G4endl;
0210 G4Exception(
0211 "IRTDamageReactionModel"
0212 "::CalculateReactionTime()",
0213 "MolecularIRTDamageReactionModel007", FatalException, exceptionDescription);
0214 }
0215
0216 auto pMoleculeA = GetMolecule(track);
0217 auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
0218 const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
0219
0220 if (pReactantList == nullptr) {
0221 G4ExceptionDescription exceptionDescription;
0222 exceptionDescription << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0223 G4cout << "!!! WARNING" << G4endl;
0224 G4cout << "IRTDamageReactionModel::CalculateReactionTime will return infinity "
0225 "for the reaction because the molecule "
0226 << pMoleculeA->GetName() << " does not have any reactants given in the reaction table."
0227 << G4endl;
0228 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0229 G4Exception(
0230 "IRTDamageReactionModel"
0231 "::CalculateReactionTime()",
0232 "MolecularIRTDamageReactionModel003", FatalException, exceptionDescription);
0233 return -1;
0234 }
0235
0236 size_t nbReactives = pReactantList->size();
0237
0238 if (nbReactives == 0) {
0239 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0240 G4cout << "!!! WARNING" << G4endl;
0241 G4cout << "IRTDamageReactionModel::CalculateReactionTime will return infinity "
0242 "for the reaction because the molecule "
0243 << pMoleculeA->GetName() << " does not have any reactants given in the reaction table."
0244 << "This message can also result from a wrong implementation of the "
0245 "reaction table."
0246 << G4endl;
0247 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0248 return -1;
0249 }
0250 const G4VTouchable* touchable = track.GetTouchable();
0251 if (touchable == nullptr) {
0252 return -1;
0253 }
0254
0255 const G4LogicalVolume* logicalVolume = touchable->GetVolume()->GetLogicalVolume();
0256 const G4ThreeVector& globalPos_track = track.GetPosition();
0257 const G4ThreeVector& localPos =
0258 touchable->GetHistory()->GetTopTransform().TransformPoint(globalPos_track);
0259
0260 G4double D = pMoleculeA->GetDiffusionCoefficient();
0261 std::vector<G4VPhysicalVolume*> result_pv;
0262 result_pv.clear();
0263 fpDNAGeometry->FindNearbyMolecules(logicalVolume, localPos, result_pv,
0264 G4IRTUtils::GetRCutOff(100 * ns));
0265
0266 if (result_pv.empty()) {
0267 return -1;
0268 }
0269 for (const auto& physicalVolume : result_pv) {
0270 const G4Material* material = physicalVolume->GetLogicalVolume()->GetMaterial();
0271
0272 G4MolecularConfiguration* dna_molConf =
0273 G4DNAMolecularMaterial::Instance()->GetMolecularConfiguration(material);
0274 auto it = std::find(pReactantList->begin(), pReactantList->end(), dna_molConf);
0275 if (it == pReactantList->end()) {
0276 continue;
0277 }
0278
0279 G4ThreeVector localPos_DNA = physicalVolume->GetTranslation();
0280 G4ThreeVector globalPos_DNA =
0281 touchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(localPos_DNA);
0282 G4double distance = (localPos_DNA - localPos).mag();
0283
0284 G4double distance2 = distance * distance;
0285 auto reactionData =
0286 G4DNAMolecularReactionTable::Instance()->GetReactionData(pMolConfA, dna_molConf);
0287
0288 G4double kobs = reactionData->GetObservedReactionRateConstant();
0289 G4double Reff = kobs / (4 * CLHEP::pi * D * CLHEP::Avogadro);
0290
0291 if (distance2 < Reff * Reff) {
0292 fminTimeStep = 0.;
0293 vp = physicalVolume;
0294 }
0295 else {
0296 G4double tempMinET = GetTimeToEncounter(pMolConfA, dna_molConf, distance);
0297 if (tempMinET > G4Scheduler::Instance()->GetEndTime() || tempMinET < 0) {
0298 continue;
0299 }
0300 if (tempMinET >= fminTimeStep) {
0301 continue;
0302 }
0303 fminTimeStep = tempMinET;
0304 vp = physicalVolume;
0305 }
0306 }
0307 if (fminTimeStep > G4Scheduler::Instance()->GetLimitingTimeStep()
0308 && fminTimeStep < G4Scheduler::Instance()->GetEndTime())
0309 {
0310 fminTimeStep = G4Scheduler::Instance()->GetLimitingTimeStep();
0311 }
0312 return fminTimeStep;
0313 }
0314
