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 // 26 // 27 // Author: Mathieu Karamitros (kara@cenbg.in2p3.fr) 28 // 29 // WARNING : This class is released as a prototype. 30 // It might strongly evolve or even disapear in the next releases. 31 // 32 // History: 33 // ----------- 34 // 10 Oct 2011 M.Karamitros created 35 // 36 // ------------------------------------------------------------------- 37 38 #include "G4DNAMolecularReaction.hh" 39 #include "G4DNAMolecularReactionTable.hh" 40 #include "G4VDNAReactionModel.hh" 41 #include "G4MolecularConfiguration.hh" 42 #include "G4Molecule.hh" 43 #include "G4MoleculeFinder.hh" 44 #include "G4ITReactionChange.hh" 45 #include "G4ITReaction.hh" 46 47 #include "G4ITTrackHolder.hh" 48 49 G4DNAMolecularReaction::G4DNAMolecularReaction() 50 : 51 fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)) 52 , fpReactionModel(nullptr) 53 { 54 } 55 56 G4DNAMolecularReaction::G4DNAMolecularReaction(G4VDNAReactionModel* pReactionModel) 57 : G4DNAMolecularReaction() 58 { 59 fpReactionModel = pReactionModel; 60 } 61 62 G4bool G4DNAMolecularReaction::TestReactibility(const G4Track &trackA, 63 const G4Track &trackB, 64 double currentStepTime, 65 bool userStepTimeLimit) /*const*/ 66 { 67 const auto pMoleculeA = GetMolecule(trackA)->GetMolecularConfiguration(); 68 const auto pMoleculeB = GetMolecule(trackB)->GetMolecularConfiguration(); 69 70 const G4double reactionRadius = fpReactionModel->GetReactionRadius(pMoleculeA, pMoleculeB); 71 72 G4double separationDistance = -1.; 73 74 if (currentStepTime == 0.) 75 { 76 userStepTimeLimit = false; 77 } 78 79 G4bool output = fpReactionModel->FindReaction(trackA, trackB, reactionRadius, 80 separationDistance, userStepTimeLimit); 81 return output; 82 } 83 84 std::unique_ptr<G4ITReactionChange> G4DNAMolecularReaction::MakeReaction(const G4Track &trackA, 85 const G4Track &trackB) 86 { 87 std::unique_ptr<G4ITReactionChange> pChanges(new G4ITReactionChange()); 88 pChanges->Initialize(trackA, trackB); 89 90 const auto pMoleculeA = GetMolecule(trackA)->GetMolecularConfiguration(); 91 const auto pMoleculeB = GetMolecule(trackB)->GetMolecularConfiguration(); 92 93 const auto pReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB); 94 95 const G4int nbProducts = pReactionData->GetNbProducts(); 96 97 if (nbProducts != 0) 98 { 99 const G4double D1 = pMoleculeA->GetDiffusionCoefficient(); 100 const G4double D2 = pMoleculeB->GetDiffusionCoefficient(); 101 const G4double sqrD1 = D1 == 0. ? 0. : std::sqrt(D1); 102 const G4double sqrD2 = D2 == 0. ? 0. : std::sqrt(D2); 103 const G4double inv_numerator = 1./(sqrD1 + sqrD2); 104 const G4ThreeVector reactionSite = sqrD2 * inv_numerator * trackA.GetPosition() 105 + sqrD1 * inv_numerator * trackB.GetPosition(); 106 107 for (G4int j = 0; j < nbProducts; ++j) 108 { 109 auto pProduct = new G4Molecule(pReactionData->GetProduct(j)); 110 auto pProductTrack = pProduct->BuildTrack(trackA.GetGlobalTime(), reactionSite); 111 112 pProductTrack->SetTrackStatus(fAlive); 113 114 G4ITTrackHolder::Instance()->Push(pProductTrack); 115 116 pChanges->AddSecondary(pProductTrack); 117 G4MoleculeFinder::Instance()->Push(pProductTrack); 118 } 119 } 120 121 pChanges->KillParents(true); 122 return pChanges; 123 } 124 125 void G4DNAMolecularReaction::SetReactionModel(G4VDNAReactionModel* pReactionModel) 126 { 127 fpReactionModel = pReactionModel; 128 } 129 130 std::vector<std::unique_ptr<G4ITReactionChange>> G4DNAMolecularReaction::FindReaction( 131 G4ITReactionSet* pReactionSet, 132 const double currentStepTime, 133 const double /*fGlobalTime*/, 134 const bool reachedUserStepTimeLimit) 135 { 136 std::vector<std::unique_ptr<G4ITReactionChange>> fReactionInfo; 137 fReactionInfo.clear(); 138 139 if (pReactionSet == nullptr) 140 { 141 return fReactionInfo; 142 } 143 144 G4ITReactionPerTrackMap& reactionPerTrackMap = pReactionSet->GetReactionMap(); 145 for (auto tracks_i = reactionPerTrackMap.cbegin(); 146 tracks_i != reactionPerTrackMap.cend(); 147 tracks_i = reactionPerTrackMap.cbegin()) 148 { 149 G4Track* pTrackA = tracks_i->first; 150 if (pTrackA->GetTrackStatus() == fStopAndKill) 151 { 152 continue; 153 } 154 155 G4ITReactionPerTrackPtr reactionPerTrack = tracks_i->second; 156 G4ITReactionList& reactionList = reactionPerTrack->GetReactionList(); 157 158 assert(reactionList.cbegin() != reactionList.cend()); 159 160 for (auto it = reactionList.cbegin(); it != reactionList.cend(); it = reactionList.cbegin()) 161 { 162 G4ITReactionPtr reaction(*it); 163 G4Track* pTrackB = reaction->GetReactant(pTrackA); 164 if (pTrackB->GetTrackStatus() == fStopAndKill) 165 { 166 continue; 167 } 168 169 if (pTrackB == pTrackA) 170 { 171 G4ExceptionDescription exceptionDescription; 172 exceptionDescription 173 << "The IT reaction process sent back a reaction between trackA and trackB. "; 174 exceptionDescription << "The problem is trackA == trackB"; 175 G4Exception("G4ITModelProcessor::FindReaction", 176 "ITModelProcessor005", 177 FatalErrorInArgument, 178 exceptionDescription); 179 } 180 181 pReactionSet->SelectThisReaction(std::move(reaction)); 182 183 if (TestReactibility(*pTrackA, *pTrackB, currentStepTime, reachedUserStepTimeLimit)) 184 { 185 auto pReactionChange = MakeReaction(*pTrackA, *pTrackB); 186 187 if (pReactionChange) 188 { 189 fReactionInfo.push_back(std::move(pReactionChange)); 190 break; 191 } 192 } 193 } 194 } 195 196 pReactionSet->CleanAllReaction(); 197 return fReactionInfo; 198 } 199