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 // 20/2/2019 26 // Author: HoangTRAN 27 28 #include "G4DNAIndependentReactionTimeStepper.hh" 29 30 #include "G4ChemicalMoleculeFinder.hh" 31 #include "G4DNAChemistryManager.hh" 32 #include "G4DNAMakeReaction.hh" 33 #include "G4DNAMolecularReactionTable.hh" 34 #include "G4DiffusionControlledReactionModel.hh" 35 #include "G4IRTUtils.hh" 36 #include "G4ITReactionChange.hh" 37 #include "G4Molecule.hh" 38 #include "G4Scheduler.hh" 39 #include "G4UnitsTable.hh" 40 #include "G4VDNAReactionModel.hh" 41 #include "Randomize.hh" 42 43 #include <memory> 44 using namespace std; 45 using namespace CLHEP; 46 47 G4DNAIndependentReactionTimeStepper::Utils::Utils(const G4Track& trackA, const G4Track& trackB) 48 : fTrackA(trackA), fTrackB(trackB) 49 { 50 fpMoleculeA = GetMolecule(trackA); 51 fpMoleculeB = GetMolecule(trackA); 52 } 53 54 G4DNAIndependentReactionTimeStepper::G4DNAIndependentReactionTimeStepper() 55 { 56 fReactionSet->SortByTime(); 57 } 58 59 void G4DNAIndependentReactionTimeStepper::Prepare() 60 { 61 G4VITTimeStepComputer::Prepare(); 62 fSampledPositions.clear(); 63 BuildChemicalMoleculeFinder() 64 } 65 66 void G4DNAIndependentReactionTimeStepper::InitializeForNewTrack() 67 { 68 if (fReactants != nullptr) { 69 fReactants.reset(); 70 } 71 fSampledMinTimeStep = DBL_MAX; 72 fHasAlreadyReachedNullTime = false; 73 } 74 75 G4double G4DNAIndependentReactionTimeStepper::CalculateStep(const G4Track& trackA, 76 const G4double& userMinTimeStep) 77 { 78 auto pMoleculeA = GetMolecule(trackA); 79 InitializeForNewTrack(); 80 fUserMinTimeStep = userMinTimeStep; 81 fCheckedTracks.insert(trackA.GetTrackID()); 82 83 #ifdef G4VERBOSE 84 if (fVerbose != 0) { 85 G4cout << "________________________________________________________________" 86 "_______" 87 << G4endl; 88 G4cout << "G4DNAIndependentReactionTimeStepper::CalculateStep" << G4endl; 89 G4cout << "Check done for molecule : " << pMoleculeA->GetName() << " (" << trackA.GetTrackID() 90 << ") " << G4endl; 91 } 92 #endif 93 94 auto pMolConfA = pMoleculeA->GetMolecularConfiguration(); 95 96 const auto pReactantList = fMolecularReactionTable->CanReactWith(pMolConfA); 97 98 if (pReactantList == nullptr) { 99 if(fVerbose > 1) { 100 G4ExceptionDescription msg; 101 msg << "G4DNAIndependentReactionTimeStepper::CalculateStep will return infinity " 102 "for the reaction because the molecule " 103 << pMoleculeA->GetName() << " does not have any reactants given in the reaction table." 104 << G4endl; 105 G4Exception("G4DNAIndependentReactionTimeStepper::CalculateStep", 106 "G4DNAIndependentReactionTimeStepper03", JustWarning, msg); 107 } 108 return DBL_MAX; 109 } 110 111 auto nbReactives = (G4int)pReactantList->size(); 112 113 if (nbReactives == 0) { 114 if(fVerbose != 0){ 115 G4ExceptionDescription msg; 116 msg << "G4DNAIndependentReactionTimeStepper::CalculateStep will " 117 "return infinity " 118 "for the reaction because the molecule " 119 << pMoleculeA->GetName() << " does not have any reactants given in the reaction table." 120 << "This message can also result from a wrong implementation of " 121 "the reaction table." 122 << G4endl; 123 G4Exception("G4DNAIndependentReactionTimeStepper::CalculateStep", 124 "G4DNAIndependentReactionTimeStepper04", JustWarning, msg); 125 } 126 return DBL_MAX; 127 } 128 fReactants = std::make_shared<vector<G4Track*>>(); 129 fReactionModel->Initialise(pMolConfA, trackA); 130 for (G4int i = 0; i < nbReactives; ++i) { 131 auto pMoleculeB = (*pReactantList)[i]; 132 G4int key = pMoleculeB->GetMoleculeID(); 133 134 // fRCutOff = G4IRTUtils::GetRCutOff(1 * ps); 135 fRCutOff = G4IRTUtils::GetRCutOff(); 136 //______________________________________________________________ 137 // Retrieve reaction range 138 const G4double Reff = fReactionModel->GetReactionRadius(i); 139 std::vector<std::pair<G4TrackList::iterator, G4double>> resultIndices; 140 resultIndices.clear(); 141 G4ChemicalMoleculeFinder::Instance()->FindNearestInRange(trackA, key, fRCutOff, resultIndices); 142 143 if (resultIndices.empty()) { 144 continue; 145 } 146 for (auto& it : resultIndices) { 147 G4Track* pTrackB = *(std::get<0>(it)); 148 149 if (pTrackB == &trackA) { 150 continue; 151 } 152 if (pTrackB == nullptr) { 153 G4ExceptionDescription exceptionDescription; 154 exceptionDescription << "No trackB no valid"; 155 G4Exception( 156 "G4DNAIndependentReactionTimeStepper" 157 "::CalculateStep()", 158 "G4DNAIndependentReactionTimeStepper007", FatalException, exceptionDescription); 159 } 160 else { 161 if (fCheckedTracks.find(pTrackB->GetTrackID()) != fCheckedTracks.end()) { 162 continue; 163 } 164 165 Utils utils(trackA, *pTrackB); 166 167 auto pMolB = GetMolecule(pTrackB); 168 auto pMolConfB = pMolB->GetMolecularConfiguration(); 169 G4double distance = (trackA.GetPosition() - pTrackB->GetPosition()).mag(); 170 if (distance * distance < Reff * Reff) { 171 auto reactionData = fMolecularReactionTable->GetReactionData(pMolConfA, pMolConfB); 172 if (G4Scheduler::Instance()->GetGlobalTime() == G4Scheduler::Instance()->GetStartTime()) { 173 if (reactionData->GetProbability() > G4UniformRand()) { 174 if (!fHasAlreadyReachedNullTime) { 175 fReactants->clear(); 176 fHasAlreadyReachedNullTime = true; 177 } 178 fSampledMinTimeStep = 0.; 179 CheckAndRecordResults(utils); 180 } 181 } 182 } 183 else { 184 G4double tempMinET = GetTimeToEncounter(trackA, *pTrackB); 185 if (tempMinET < 0 || tempMinET > G4Scheduler::Instance()->GetEndTime()) { 186 continue; 187 } 188 if (tempMinET >= fSampledMinTimeStep) { 189 continue; 190 } 191 fSampledMinTimeStep = tempMinET; 192 fReactants->clear(); 193 CheckAndRecordResults(utils); 194 } 195 } 196 } 197 } 198 199 #ifdef G4VERBOSE 200 if (fVerbose != 0) { 201 G4cout << "G4DNAIndependentReactionTimeStepper::CalculateStep will finally " 202 "return :" 203 << G4BestUnit(fSampledMinTimeStep, "Time") << G4endl; 204 205 if (fVerbose > 1) { 206 G4cout << "Selected reactants for trackA: " << pMoleculeA->GetName() << " (" 207 << trackA.GetTrackID() << ") are: "; 208 209 vector<G4Track*>::iterator it; 210 for (it = fReactants->begin(); it != fReactants->end(); it++) { 211 G4Track* trackB = *it; 212 G4cout << GetMolecule(trackB)->GetName() << " (" << trackB->GetTrackID() << ") \t "; 213 } 214 G4cout << G4endl; 215 } 216 } 217 #endif 218 return fSampledMinTimeStep; 219 } 220 221 void G4DNAIndependentReactionTimeStepper::CheckAndRecordResults(const Utils& utils) 222 { 223 if (utils.fTrackB.GetTrackStatus() != fAlive) { 224 return; 225 } 226 227 if (&utils.fTrackB == &utils.fTrackA) { 228 G4ExceptionDescription msg; 229 msg << "A track is reacting with itself" 230 " (which is impossible) ie fpTrackA == trackB" 231 << G4endl; 232 msg << "Molecule A is of type : " << utils.fpMoleculeA->GetName() 233 << " with trackID : " << utils.fTrackA.GetTrackID() 234 << " and B : " << utils.fpMoleculeB->GetName() 235 << " with trackID : " << utils.fTrackB.GetTrackID() << G4endl; 236 G4Exception("G4DNAIndependentReactionTimeStepper::RetrieveResults", 237 "G4DNAIndependentReactionTimeStepper003", FatalErrorInArgument, 238 msg); 239 } 240 241 if (fabs(utils.fTrackB.GetGlobalTime() - utils.fTrackA.GetGlobalTime()) 242 > utils.fTrackA.GetGlobalTime() * (1. - 1. / 100)) 243 { 244 // DEBUG 245 G4ExceptionDescription msg; 246 msg << "The interacting tracks are not synchronized in time" << G4endl; 247 msg << "trackB->GetGlobalTime() != fpTrackA.GetGlobalTime()" << G4endl; 248 249 msg << "fpTrackA : trackID : " << utils.fTrackA.GetTrackID() 250 << "\t Name :" << utils.fpMoleculeA->GetName() 251 << "\t fpTrackA->GetGlobalTime() = " 252 << G4BestUnit(utils.fTrackA.GetGlobalTime(), "Time") << G4endl; 253 254 msg << "trackB : trackID : " << utils.fTrackB.GetTrackID() 255 << "\t Name :" << utils.fpMoleculeB->GetName() 256 << "\t trackB->GetGlobalTime() = " 257 << G4BestUnit(utils.fTrackB.GetGlobalTime(), "Time") << G4endl; 258 259 G4Exception("G4DNAIndependentReactionTimeStepper::RetrieveResults", 260 "G4DNAIndependentReactionTimeStepper004", FatalErrorInArgument, 261 msg); 262 } 263 fReactants->push_back(const_cast<G4Track*>(&utils.fTrackB)); 264 } 265 266 std::unique_ptr<G4ITReactionChange> G4DNAIndependentReactionTimeStepper::FindReaction( 267 G4ITReactionSet* pReactionSet, const G4double& currentStepTime, 268 const G4double& /*previousStepTime*/, const G4bool& /*reachedUserStepTimeLimit*/) 269 { 270 if (pReactionSet == nullptr) { 271 return nullptr; 272 } 273 274 G4ITReactionPerTime& reactionPerTime = pReactionSet->GetReactionsPerTime(); 275 if (reactionPerTime.empty()) { 276 return nullptr; 277 } 278 for (auto reaction_i = reactionPerTime.begin(); reaction_i != reactionPerTime.end(); 279 reaction_i = reactionPerTime.begin()) 280 { 281 if ((*reaction_i)->GetTime() > currentStepTime) { 282 fReactionSet->CleanAllReaction(); 283 return nullptr; 284 } 285 286 G4Track* pTrackA = (*reaction_i)->GetReactants().first; 287 if (pTrackA->GetTrackStatus() == fStopAndKill) { 288 continue; 289 } 290 G4Track* pTrackB = (*reaction_i)->GetReactant(pTrackA); 291 if (pTrackB->GetTrackStatus() == fStopAndKill) { 292 continue; 293 } 294 295 if (pTrackB == pTrackA) { 296 G4ExceptionDescription msg; 297 msg << "The IT reaction process sent back a reaction " 298 "between trackA and trackB. "; 299 msg << "The problem is trackA == trackB"; 300 G4Exception("G4DNAIndependentReactionTimeStepper::FindReaction", 301 "G4DNAIndependentReactionTimeStepper02", FatalErrorInArgument, 302 msg); 303 } 304 pReactionSet->SelectThisReaction(*reaction_i); 305 if (fpReactionProcess != nullptr 306 && fpReactionProcess->TestReactibility(*pTrackA, *pTrackB, currentStepTime, false)) 307 { 308 if ((fSampledPositions.find(pTrackA->GetTrackID()) == fSampledPositions.end() 309 && (fSampledPositions.find(pTrackB->GetTrackID()) == fSampledPositions.end()))) 310 { 311 G4ExceptionDescription msg; 312 msg << "The positions of trackA and trackB have no counted "; 313 G4Exception("G4DNAIndependentReactionTimeStepper::FindReaction", 314 "G4DNAIndependentReactionTimeStepper0001", FatalErrorInArgument, 315 msg); 316 } 317 318 pTrackA->SetPosition(fSampledPositions[pTrackA->GetTrackID()]); 319 pTrackB->SetPosition(fSampledPositions[pTrackB->GetTrackID()]); 320 auto pReactionChange = fpReactionProcess->MakeReaction(*pTrackA, *pTrackB); 321 if (pReactionChange == nullptr) { 322 return nullptr; 323 } 324 return pReactionChange; 325 } 326 } 327 return nullptr; 328 } 329 330 void G4DNAIndependentReactionTimeStepper::SetReactionModel(G4VDNAReactionModel* pReactionModel) 331 { 332 fReactionModel = pReactionModel; 333 } 334 335 G4VDNAReactionModel* G4DNAIndependentReactionTimeStepper::GetReactionModel() 336 { 337 return fReactionModel; 338 } 339 340 void G4DNAIndependentReactionTimeStepper::SetVerbose(G4int flag) 341 { 342 fVerbose = flag; 343 } 344 345 G4double G4DNAIndependentReactionTimeStepper::GetTimeToEncounter(const G4Track& trackA, 346 const G4Track& trackB) 347 { 348 G4double timeToReaction = dynamic_cast<G4DiffusionControlledReactionModel*>(fReactionModel) 349 ->GetTimeToEncounter(trackA, trackB); 350 return timeToReaction; 351 } 352 353 void G4DNAIndependentReactionTimeStepper::SetReactionProcess(G4VITReactionProcess* pReactionProcess) 354 { 355 fpReactionProcess = pReactionProcess; 356 } 357 G4double G4DNAIndependentReactionTimeStepper::CalculateMinTimeStep(G4double /*currentGlobalTime*/, 358 G4double definedMinTimeStep) 359 { 360 G4double fTSTimeStep = DBL_MAX; 361 fCheckedTracks.clear(); 362 363 for (auto pTrack : *fpTrackContainer->GetMainList()) { 364 if (pTrack == nullptr) { 365 G4ExceptionDescription msg; 366 msg << "No track found."; 367 G4Exception("G4DNAIndependentReactionTimeStepper::CalculateMinTimeStep", 368 "G4DNAIndependentReactionTimeStepper006", FatalErrorInArgument, 369 msg); 370 continue; 371 } 372 373 G4TrackStatus trackStatus = pTrack->GetTrackStatus(); 374 if (trackStatus == fStopAndKill || trackStatus == fStopButAlive) { 375 continue; 376 } 377 378 G4double sampledMinTimeStep = CalculateStep(*pTrack, definedMinTimeStep); 379 G4TrackVectorHandle reactants = GetReactants(); 380 381 if (sampledMinTimeStep < fTSTimeStep) { 382 fTSTimeStep = sampledMinTimeStep; 383 if (reactants) { 384 fReactionSet->AddReactions(fTSTimeStep, const_cast<G4Track*>(pTrack), std::move(reactants)); 385 386 fSampledPositions[pTrack->GetTrackID()] = pTrack->GetPosition(); 387 for (const auto& it : *fReactants) { 388 auto pTrackB = it; 389 fSampledPositions[pTrackB->GetTrackID()] = pTrackB->GetPosition(); 390 } 391 ResetReactants(); 392 } 393 } 394 else if (fTSTimeStep == sampledMinTimeStep && G4bool(reactants)) { 395 fReactionSet->AddReactions(fTSTimeStep, const_cast<G4Track*>(pTrack), std::move(reactants)); 396 397 fSampledPositions[pTrack->GetTrackID()] = pTrack->GetPosition(); 398 for (const auto& it : *fReactants) { 399 auto pTrackB = it; 400 fSampledPositions[pTrackB->GetTrackID()] = pTrackB->GetPosition(); 401 } 402 ResetReactants(); 403 } 404 else if (reactants) { 405 ResetReactants(); 406 } 407 } 408 return fTSTimeStep; 409 } 410