Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4DNASmoluchowskiReactionModel.cc 85244 2014-10-27 08:24:13Z gcosmo $ 26 // 27 // 27 #include "G4DNASmoluchowskiReactionModel.hh" 28 #include "G4DNASmoluchowskiReactionModel.hh" 28 #include "Randomize.hh" 29 #include "Randomize.hh" 29 #include "G4Track.hh" 30 #include "G4Track.hh" 30 #include "G4DNAMolecularReactionTable.hh" 31 #include "G4DNAMolecularReactionTable.hh" 31 #include "G4UnitsTable.hh" 32 #include "G4UnitsTable.hh" 32 #include "G4Molecule.hh" << 33 #include "G4Exp.hh" << 34 33 35 G4DNASmoluchowskiReactionModel::G4DNASmoluchow << 34 G4DNASmoluchowskiReactionModel::G4DNASmoluchowskiReactionModel() : >> 35 G4VDNAReactionModel() >> 36 { >> 37 fReactionData = 0; >> 38 } >> 39 >> 40 G4DNASmoluchowskiReactionModel::G4DNASmoluchowskiReactionModel(const G4DNASmoluchowskiReactionModel& __right) : >> 41 G4VDNAReactionModel(__right) >> 42 { >> 43 fReactionData = 0; >> 44 } 36 45 37 G4DNASmoluchowskiReactionModel::~G4DNASmolucho << 46 G4DNASmoluchowskiReactionModel& G4DNASmoluchowskiReactionModel::operator=(const G4DNASmoluchowskiReactionModel& right) >> 47 { >> 48 if (this == &right) return *this; >> 49 fReactionData = 0; >> 50 return *this; >> 51 } >> 52 >> 53 G4DNASmoluchowskiReactionModel::~G4DNASmoluchowskiReactionModel() >> 54 { >> 55 fReactionData = 0; >> 56 } 38 57 39 void G4DNASmoluchowskiReactionModel::Initialis << 58 void G4DNASmoluchowskiReactionModel::Initialise(const G4Molecule* __molecule, 40 59 const G4Track&) 41 { 60 { 42 fpReactionData = fpReactionTable->GetReact << 61 fReactionData = fReactionTable->GetReactionData(__molecule); 43 } 62 } 44 63 45 void G4DNASmoluchowskiReactionModel::Initialis << 64 void G4DNASmoluchowskiReactionModel::InitialiseToPrint(const G4Molecule* __molecule) 46 { 65 { 47 fpReactionData = fpReactionTable->GetReact << 66 fReactionData = fReactionTable->GetReactionData(__molecule); 48 } 67 } 49 68 50 G4double G4DNASmoluchowskiReactionModel::GetRe << 69 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4Molecule* mol1, 51 << 70 const G4Molecule* mol2) 52 { 71 { 53 G4double __output = fpReactionTable->GetRe << 72 G4double __output = fReactionTable->GetReactionData(mol1, mol2) 54 return __output; << 73 ->GetReducedReactionRadius(); >> 74 return __output; 55 } 75 } 56 76 57 G4double G4DNASmoluchowskiReactionModel::GetRe << 77 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4int __i) 58 { 78 { 59 G4double __output = (*fpReactionData)[__i] << 79 G4double __output = (*fReactionData)[__i]->GetReducedReactionRadius(); 60 return __output; << 80 return __output; 61 } 81 } 62 82 63 G4bool G4DNASmoluchowskiReactionModel::FindRea 83 G4bool G4DNASmoluchowskiReactionModel::FindReaction(const G4Track& __trackA, 64 84 const G4Track& __trackB, 65 << 85 const G4double __R, 66 << 86 G4double& __r, 67 87 const G4bool __alongStepReaction) 68 { 88 { 69 const G4double R2 = __reactionRadius * __r << 89 G4double __postStepSeparation = 0; 70 G4double postStepSeparation = 0.; << 90 bool __do_break = false; 71 bool do_break = false; << 91 G4double __R2 = __R * __R; 72 int k = 0; << 92 int k = 0; >> 93 >> 94 for (; k < 3; k++) >> 95 { >> 96 __postStepSeparation += std::pow( >> 97 __trackA.GetPosition()[k] - __trackB.GetPosition()[k], 2); 73 98 74 for (; k < 3; ++k) << 99 if (__postStepSeparation > __R2) 75 { 100 { 76 postStepSeparation += std::pow( << 101 __do_break = true; 77 __trackA.GetPosition()[k] - __trac << 102 break; 78 << 79 if (postStepSeparation > R2) << 80 { << 81 do_break = true; << 82 break; << 83 } << 84 } 103 } >> 104 } >> 105 >> 106 if (__do_break == false) >> 107 { >> 108 // The loop was not break >> 109 // => __r^2 < __R^2 >> 110 __r = std::sqrt(__postStepSeparation); >> 111 return true; >> 112 } >> 113 else if (__alongStepReaction == true) >> 114 { >> 115 //G4cout << "alongStepReaction==true" << G4endl; >> 116 //Along step cheack and >> 117 // the loop has break 85 118 86 if (!do_break) << 119 // Continue loop >> 120 for (; k < 3; k++) 87 { 121 { 88 // The loop was not break << 122 __postStepSeparation += std::pow( 89 // => r^2 < R^2 << 123 __trackA.GetPosition()[k] - __trackB.GetPosition()[k], 2); 90 __separationDistance = std::sqrt(postS << 91 return true; << 92 } 124 } 93 if (__alongStepReaction) << 125 // Use Green approach : the Brownian bridge 94 { << 126 __r = (__postStepSeparation = std::sqrt(__postStepSeparation)); 95 //Along step check and the loop has br << 127 >> 128 G4Molecule* __moleculeA = GetMolecule(__trackA); >> 129 G4Molecule* __moleculeB = GetMolecule(__trackB); >> 130 >> 131 G4double __D = __moleculeA->GetDiffusionCoefficient() >> 132 + __moleculeB->GetDiffusionCoefficient(); 96 133 97 // Continue loop << 134 G4ThreeVector __preStepPositionA = __trackA.GetStep()->GetPreStepPoint() 98 for (; k < 3; ++k) << 135 ->GetPosition(); 99 { << 136 G4ThreeVector __preStepPositionB = __trackB.GetStep()->GetPreStepPoint() 100 postStepSeparation += std::pow( << 137 ->GetPosition(); 101 __trackA.GetPosition()[k] - __ << 138 102 } << 139 if (__preStepPositionA == __trackA.GetPosition()) 103 // Use Green approach : the Brownian b << 140 { 104 __separationDistance = (postStepSepara << 141 G4ExceptionDescription exceptionDescription; 105 << 142 exceptionDescription << "The molecule : " << __moleculeA->GetName(); 106 auto pMoleculeA = GetMolecule(__trackA << 143 exceptionDescription << " with track ID :" << __trackA.GetTrackID(); 107 auto pMoleculeB = GetMolecule(__trackB << 144 exceptionDescription << " did not move since the previous step." << G4endl; 108 << 145 exceptionDescription << "Current position : " 109 G4double D = pMoleculeA->GetDiffusionC << 146 << G4BestUnit(__trackA.GetPosition(), "Length") 110 + pMoleculeB->GetDiffusio << 147 << G4endl; 111 << 148 exceptionDescription << "Previous position : " 112 const auto& preStepPositionA = __track << 149 << G4BestUnit(__preStepPositionA, "Length") << G4endl; 113 const auto& preStepPositionB = __track << 150 G4Exception("G4DNASmoluchowskiReactionModel::FindReaction", 114 << 151 "G4DNASmoluchowskiReactionModel", FatalErrorInArgument, 115 G4double preStepSeparation = (preStepP << 152 exceptionDescription); 116 << 117 //=================================== << 118 // Brownian bridge << 119 G4double __probabiltyOfEncounter = G4E << 120 << 121 << 122 ); << 123 G4double __selectedPOE = G4UniformRand << 124 << 125 if (__selectedPOE <= __probabiltyOfEnc << 126 { << 127 return true; << 128 } << 129 //=================================== << 130 } 153 } 131 154 132 return false; << 155 G4double __preStepSeparation = >> 156 (__preStepPositionA - __preStepPositionB).mag(); >> 157 >> 158 G4double __probabiltyOfEncounter = std::exp( >> 159 -(__preStepSeparation - __R) * (__postStepSeparation - __R) / (__D >> 160 * (__trackB.GetStep()->GetDeltaTime()))); >> 161 G4double __selectedPOE = G4UniformRand(); >> 162 >> 163 if (__selectedPOE <= __probabiltyOfEncounter) return true; >> 164 } >> 165 >> 166 return false; 133 } 167 } 134 168