Warning, file /geant4/examples/advanced/dna/moleculardna/src/IRTDamageReactionModel.cc was not indexed
or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).
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 fpTrack = &track;
0188 fpDNAPhyVolume = std::get<const G4VPhysicalVolume*>(vp);
0189 MakeReaction(track);
0190 RecordDNADamage();
0191 G4Scheduler::Instance()->SetInteractionStep(true);
0192 return true;
0193 }
0194
0195
0196
0197 G4double IRTDamageReactionModel::CalculateReactionTime(const G4Track& track, DNANode& vp)
0198 {
0199 fpTrack = nullptr;
0200 fminTimeStep = DBL_MAX;
0201 fReactionTime = DBL_MAX;
0202 fpDNAPhyVolume = nullptr;
0203 if (fpDNAGeometry == nullptr) {
0204 G4ExceptionDescription exceptionDescription;
0205 exceptionDescription << "no fpDNAGeometry" << G4endl;
0206 G4Exception(
0207 "IRTDamageReactionModel"
0208 "::CalculateReactionTime()",
0209 "MolecularIRTDamageReactionModel007", FatalException, exceptionDescription);
0210 }
0211
0212 auto pMoleculeA = GetMolecule(track);
0213 auto pMolConfA = pMoleculeA->GetMolecularConfiguration();
0214 const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA);
0215
0216 if (pReactantList == nullptr) {
0217 G4ExceptionDescription exceptionDescription;
0218 exceptionDescription << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0219 G4cout << "!!! WARNING" << G4endl;
0220 G4cout << "IRTDamageReactionModel::CalculateReactionTime will return infinity "
0221 "for the reaction because the molecule "
0222 << pMoleculeA->GetName() << " does not have any reactants given in the reaction table."
0223 << G4endl;
0224 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0225 G4Exception(
0226 "IRTDamageReactionModel"
0227 "::CalculateReactionTime()",
0228 "MolecularIRTDamageReactionModel003", FatalException, exceptionDescription);
0229 return -1;
0230 }
0231
0232 size_t nbReactives = pReactantList->size();
0233
0234 if (nbReactives == 0) {
0235 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0236 G4cout << "!!! WARNING" << G4endl;
0237 G4cout << "IRTDamageReactionModel::CalculateReactionTime will return infinity "
0238 "for the reaction because the molecule "
0239 << pMoleculeA->GetName() << " does not have any reactants given in the reaction table."
0240 << "This message can also result from a wrong implementation of the "
0241 "reaction table."
0242 << G4endl;
0243 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0244 return -1;
0245 }
0246 const G4VTouchable* touchable = track.GetTouchable();
0247 if (touchable == nullptr) {
0248 return -1;
0249 }
0250
0251 const G4LogicalVolume* logicalVolume = touchable->GetVolume()->GetLogicalVolume();
0252 const G4ThreeVector& globalPos_track = track.GetPosition();
0253 const G4ThreeVector& localPos =
0254 touchable->GetHistory()->GetTopTransform().TransformPoint(globalPos_track);
0255
0256 G4double D = pMoleculeA->GetDiffusionCoefficient();
0257 std::vector<G4VPhysicalVolume*> result_pv;
0258 result_pv.clear();
0259 fpDNAGeometry->FindNearbyMolecules(logicalVolume, localPos, result_pv,
0260 G4IRTUtils::GetRCutOff(100 * ns));
0261
0262 if (result_pv.empty()) {
0263 return -1;
0264 }
0265 for (const auto& physicalVolume : result_pv) {
0266 const G4Material* material = physicalVolume->GetLogicalVolume()->GetMaterial();
0267
0268 G4MolecularConfiguration* dna_molConf =
0269 G4DNAMolecularMaterial::Instance()->GetMolecularConfiguration(material);
0270 auto it = std::find(pReactantList->begin(), pReactantList->end(), dna_molConf);
0271 if (it == pReactantList->end()) {
0272 continue;
0273 }
0274
0275 G4ThreeVector localPos_DNA = physicalVolume->GetTranslation();
0276 G4ThreeVector globalPos_DNA =
0277 touchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(localPos_DNA);
0278 G4double distance = (localPos_DNA - localPos).mag();
0279
0280 G4double distance2 = distance * distance;
0281 auto reactionData =
0282 G4DNAMolecularReactionTable::Instance()->GetReactionData(pMolConfA, dna_molConf);
0283
0284 G4double kobs = reactionData->GetObservedReactionRateConstant();
0285 G4double Reff = kobs / (4 * CLHEP::pi * D * CLHEP::Avogadro);
0286
0287 if (distance2 < Reff * Reff) {
0288 fminTimeStep = 0.;
0289 vp = physicalVolume;
0290 }
0291 else {
0292 G4double tempMinET = GetTimeToEncounter(pMolConfA, dna_molConf, distance);
0293 if (tempMinET > G4Scheduler::Instance()->GetEndTime() || tempMinET < 0) {
0294 continue;
0295 }
0296 if (tempMinET >= fminTimeStep) {
0297 continue;
0298 }
0299 fminTimeStep = tempMinET;
0300 vp = physicalVolume;
0301 }
0302 }
0303 return fminTimeStep;
0304 }
0305
