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 // 20/2/2019 26 // Author: HoangTRAN 27 28 #include "G4DNAIndependentReactionTimeStepper. 29 30 #include "G4ChemicalMoleculeFinder.hh" 31 #include "G4DNAChemistryManager.hh" 32 #include "G4DNAMakeReaction.hh" 33 #include "G4DNAMolecularReactionTable.hh" 34 #include "G4DiffusionControlledReactionModel.h 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::Ut 48 : fTrackA(trackA), fTrackB(trackB) 49 { 50 fpMoleculeA = GetMolecule(trackA); 51 fpMoleculeB = GetMolecule(trackA); 52 } 53 54 G4DNAIndependentReactionTimeStepper::G4DNAInde 55 { 56 fReactionSet->SortByTime(); 57 } 58 59 void G4DNAIndependentReactionTimeStepper::Prep 60 { 61 G4VITTimeStepComputer::Prepare(); 62 fSampledPositions.clear(); 63 BuildChemicalMoleculeFinder() 64 } 65 66 void G4DNAIndependentReactionTimeStepper::Init 67 { 68 if (fReactants != nullptr) { 69 fReactants.reset(); 70 } 71 fSampledMinTimeStep = DBL_MAX; 72 fHasAlreadyReachedNullTime = false; 73 } 74 75 G4double G4DNAIndependentReactionTimeStepper:: 76 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 << "G4DNAIndependentReactionTimeSte 89 G4cout << "Check done for molecule : " << 90 << ") " << G4endl; 91 } 92 #endif 93 94 auto pMolConfA = pMoleculeA->GetMolecularCon 95 96 const auto pReactantList = fMolecularReactio 97 98 if (pReactantList == nullptr) { 99 if(fVerbose > 1) { 100 G4ExceptionDescription msg; 101 msg << "G4DNAIndependentReactionTimeStep 102 "for the reaction because the mol 103 << pMoleculeA->GetName() << " does n 104 << G4endl; 105 G4Exception("G4DNAIndependentReactionTim 106 "G4DNAIndependentReactionTim 107 } 108 return DBL_MAX; 109 } 110 111 auto nbReactives = (G4int)pReactantList->siz 112 113 if (nbReactives == 0) { 114 if(fVerbose != 0){ 115 G4ExceptionDescription msg; 116 msg << "G4DNAIndependentReactionTimeStep 117 "return infinity " 118 "for the reaction because the mol 119 << pMoleculeA->GetName() << " does n 120 << "This message can also result fro 121 "the reaction table." 122 << G4endl; 123 G4Exception("G4DNAIndependentReactionTim 124 "G4DNAIndependentReactionTim 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 * p 135 fRCutOff = G4IRTUtils::GetRCutOff(); 136 //________________________________________ 137 // Retrieve reaction range 138 const G4double Reff = fReactionModel->GetR 139 std::vector<std::pair<G4TrackList::iterato 140 resultIndices.clear(); 141 G4ChemicalMoleculeFinder::Instance()->Find 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 exceptionDescri 154 exceptionDescription << "No trackB no 155 G4Exception( 156 "G4DNAIndependentReactionTimeStepper 157 "::CalculateStep()", 158 "G4DNAIndependentReactionTimeStepper 159 } 160 else { 161 if (fCheckedTracks.find(pTrackB->GetTr 162 continue; 163 } 164 165 Utils utils(trackA, *pTrackB); 166 167 auto pMolB = GetMolecule(pTrackB); 168 auto pMolConfB = pMolB->GetMolecularCo 169 G4double distance = (trackA.GetPositio 170 if (distance * distance < Reff * Reff) 171 auto reactionData = fMolecularReacti 172 if (G4Scheduler::Instance()->GetGlob 173 if (reactionData->GetProbability() 174 if (!fHasAlreadyReachedNullTime) 175 fReactants->clear(); 176 fHasAlreadyReachedNullTime = t 177 } 178 fSampledMinTimeStep = 0.; 179 CheckAndRecordResults(utils); 180 } 181 } 182 } 183 else { 184 G4double tempMinET = GetTimeToEncoun 185 if (tempMinET < 0 || tempMinET > G4S 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 << "G4DNAIndependentReactionTimeSte 202 "return :" 203 << G4BestUnit(fSampledMinTimeStep, 204 205 if (fVerbose > 1) { 206 G4cout << "Selected reactants for trackA 207 << trackA.GetTrackID() << ") are: 208 209 vector<G4Track*>::iterator it; 210 for (it = fReactants->begin(); it != fRe 211 G4Track* trackB = *it; 212 G4cout << GetMolecule(trackB)->GetName 213 } 214 G4cout << G4endl; 215 } 216 } 217 #endif 218 return fSampledMinTimeStep; 219 } 220 221 void G4DNAIndependentReactionTimeStepper::Chec 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 imposs 231 << G4endl; 232 msg << "Molecule A is of type : " << utils 233 << " with trackID : " 234 << " and B : " << uti 235 << " with trackID : " 236 G4Exception("G4DNAIndependentReactionTimeS 237 "G4DNAIndependentReactionTimeS 238 msg); 239 } 240 241 if (fabs(utils.fTrackB.GetGlobalTime() - uti 242 > utils.fTrackA.GetGlobalTime() * (1. - 243 { 244 // DEBUG 245 G4ExceptionDescription msg; 246 msg << "The interacting tracks are not syn 247 msg << "trackB->GetGlobalTime() != fpTrack 248 249 msg << "fpTrackA : trackID : " << utils.fT 250 << "\t Name :" << uti 251 << "\t fpTrackA->GetG 252 << G4BestUnit(utils.f 253 254 msg << "trackB : trackID : " << utils.fTra 255 << "\t Name :" << uti 256 << "\t trackB->GetGlo 257 << G4BestUnit(utils.f 258 259 G4Exception("G4DNAIndependentReactionTimeS 260 "G4DNAIndependentReactionTimeS 261 msg); 262 } 263 fReactants->push_back(const_cast<G4Track*>(& 264 } 265 266 std::unique_ptr<G4ITReactionChange> G4DNAIndep 267 G4ITReactionSet* pReactionSet, const G4doubl 268 const G4double& /*previousStepTime*/, const 269 { 270 if (pReactionSet == nullptr) { 271 return nullptr; 272 } 273 274 G4ITReactionPerTime& reactionPerTime = pReac 275 if (reactionPerTime.empty()) { 276 return nullptr; 277 } 278 for (auto reaction_i = reactionPerTime.begin 279 reaction_i = reactionPerTime.begin()) 280 { 281 if ((*reaction_i)->GetTime() > currentStep 282 fReactionSet->CleanAllReaction(); 283 return nullptr; 284 } 285 286 G4Track* pTrackA = (*reaction_i)->GetReact 287 if (pTrackA->GetTrackStatus() == fStopAndK 288 continue; 289 } 290 G4Track* pTrackB = (*reaction_i)->GetReact 291 if (pTrackB->GetTrackStatus() == fStopAndK 292 continue; 293 } 294 295 if (pTrackB == pTrackA) { 296 G4ExceptionDescription msg; 297 msg << "The IT reaction process sent bac 298 "between trackA 299 msg << "The problem is trackA == trackB" 300 G4Exception("G4DNAIndependentReactionTim 301 "G4DNAIndependentReactionTim 302 msg); 303 } 304 pReactionSet->SelectThisReaction(*reaction 305 if (fpReactionProcess != nullptr 306 && fpReactionProcess->TestReactibility 307 { 308 if ((fSampledPositions.find(pTrackA->Get 309 && (fSampledPositions.find(pTrackB- 310 { 311 G4ExceptionDescription msg; 312 msg << "The positions of trackA and tr 313 G4Exception("G4DNAIndependentReactionT 314 "G4DNAIndependentReactionT 315 msg); 316 } 317 318 pTrackA->SetPosition(fSampledPositions[p 319 pTrackB->SetPosition(fSampledPositions[p 320 auto pReactionChange = fpReactionProcess 321 if (pReactionChange == nullptr) { 322 return nullptr; 323 } 324 return pReactionChange; 325 } 326 } 327 return nullptr; 328 } 329 330 void G4DNAIndependentReactionTimeStepper::SetR 331 { 332 fReactionModel = pReactionModel; 333 } 334 335 G4VDNAReactionModel* G4DNAIndependentReactionT 336 { 337 return fReactionModel; 338 } 339 340 void G4DNAIndependentReactionTimeStepper::SetV 341 { 342 fVerbose = flag; 343 } 344 345 G4double G4DNAIndependentReactionTimeStepper:: 346 347 { 348 G4double timeToReaction = dynamic_cast<G4Dif 349 ->GetTimeToEncou 350 return timeToReaction; 351 } 352 353 void G4DNAIndependentReactionTimeStepper::SetR 354 { 355 fpReactionProcess = pReactionProcess; 356 } 357 G4double G4DNAIndependentReactionTimeStepper:: 358 359 { 360 G4double fTSTimeStep = DBL_MAX; 361 fCheckedTracks.clear(); 362 363 for (auto pTrack : *fpTrackContainer->GetMai 364 if (pTrack == nullptr) { 365 G4ExceptionDescription msg; 366 msg << "No track found."; 367 G4Exception("G4DNAIndependentReactionTim 368 "G4DNAIndependentReactionTim 369 msg); 370 continue; 371 } 372 373 G4TrackStatus trackStatus = pTrack->GetTra 374 if (trackStatus == fStopAndKill || trackSt 375 continue; 376 } 377 378 G4double sampledMinTimeStep = CalculateSte 379 G4TrackVectorHandle reactants = GetReactan 380 381 if (sampledMinTimeStep < fTSTimeStep) { 382 fTSTimeStep = sampledMinTimeStep; 383 if (reactants) { 384 fReactionSet->AddReactions(fTSTimeStep 385 386 fSampledPositions[pTrack->GetTrackID() 387 for (const auto& it : *fReactants) { 388 auto pTrackB = it; 389 fSampledPositions[pTrackB->GetTrackI 390 } 391 ResetReactants(); 392 } 393 } 394 else if (fTSTimeStep == sampledMinTimeStep 395 fReactionSet->AddReactions(fTSTimeStep, 396 397 fSampledPositions[pTrack->GetTrackID()] 398 for (const auto& it : *fReactants) { 399 auto pTrackB = it; 400 fSampledPositions[pTrackB->GetTrackID( 401 } 402 ResetReactants(); 403 } 404 else if (reactants) { 405 ResetReactants(); 406 } 407 } 408 return fTSTimeStep; 409 } 410