00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "G4DNASmoluchowskiReactionModel.hh"
00029 #include "Randomize.hh"
00030 #include "G4Track.hh"
00031 #include "G4DNAMolecularReactionTable.hh"
00032
00033 G4DNASmoluchowskiReactionModel::G4DNASmoluchowskiReactionModel() : G4VDNAReactionModel()
00034 {
00035 fReactionData = 0 ;
00036 }
00037
00038 G4DNASmoluchowskiReactionModel::G4DNASmoluchowskiReactionModel(const G4DNASmoluchowskiReactionModel& __right) :
00039 G4VDNAReactionModel(__right)
00040 {
00041 fReactionData = 0 ;
00042 }
00043
00044 G4DNASmoluchowskiReactionModel& G4DNASmoluchowskiReactionModel::operator=(const G4DNASmoluchowskiReactionModel& right)
00045 {
00046 if(this == &right) return *this;
00047 fReactionData = 0;
00048 return *this;
00049 }
00050
00051 G4DNASmoluchowskiReactionModel::~G4DNASmoluchowskiReactionModel()
00052 {
00053 fReactionData = 0 ;
00054 }
00055
00056 void G4DNASmoluchowskiReactionModel::Initialise(const G4Molecule* __molecule, const G4Track&)
00057 {
00058 fReactionData = fReactionTable->GetReactionData(__molecule);
00059 }
00060
00061 void G4DNASmoluchowskiReactionModel::InitialiseToPrint(const G4Molecule* __molecule)
00062 {
00063 fReactionData = fReactionTable->GetReactionData(__molecule);
00064 }
00065
00066 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4Molecule* mol1,
00067 const G4Molecule* mol2)
00068 {
00069 G4double __output = fReactionTable -> GetReactionData(mol1,mol2)->GetReducedReactionRadius();
00070 return __output ;
00071 }
00072
00073 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4int __i)
00074 {
00075 G4double __output = (*fReactionData)[__i] -> GetReducedReactionRadius();
00076 return __output ;
00077 }
00078
00079 G4bool G4DNASmoluchowskiReactionModel::FindReaction(const G4Track& __trackA,
00080 const G4Track& __trackB,
00081 const G4double __R,
00082 G4double& __r,
00083 const G4bool __alongStepReaction)
00084 {
00085 G4double __postStepSeparation = 0;
00086 bool __do_break = false ;
00087 G4double __R2 = __R*__R ;
00088 int k = 0 ;
00089
00090 for(; k < 3 ; k++)
00091 {
00092 __postStepSeparation += std::pow(__trackA.GetPosition()[k] - __trackB.GetPosition()[k],2);
00093
00094 if(__postStepSeparation > __R2)
00095 {
00096 __do_break = true ;
00097 break ;
00098 }
00099 }
00100
00101 if(__do_break == false)
00102 {
00103
00104
00105 __r = std::sqrt(__postStepSeparation);
00106 return true;
00107 }
00108 else if(__alongStepReaction == true)
00109 {
00110
00111
00112
00113
00114
00115 for(; k < 3 ; k++)
00116 {
00117 __postStepSeparation += std::pow(__trackA.GetPosition()[k] - __trackB.GetPosition()[k],2);
00118 }
00119
00120 __r = (__postStepSeparation = std::sqrt(__postStepSeparation) );
00121
00122 G4Molecule* __moleculeA = GetMolecule(__trackA);
00123 G4Molecule* __moleculeB = GetMolecule(__trackB);
00124
00125 G4double __D = __moleculeA->GetDiffusionCoefficient() + __moleculeB->GetDiffusionCoefficient();
00126
00127 G4ThreeVector __preStepPositionA = __trackA.GetStep()->GetPreStepPoint() ->GetPosition();
00128 G4ThreeVector __preStepPositionB = __trackB.GetStep()->GetPreStepPoint() ->GetPosition();
00129
00130 if(__preStepPositionA == __trackA.GetPosition())
00131 {
00132 G4ExceptionDescription exceptionDescription ;
00133 exceptionDescription << "The molecule : " << __moleculeA->GetName();
00134 exceptionDescription << " did not move since the previous step ";
00135 G4Exception("G4DNASmoluchowskiReactionModel::FindReaction","G4DNASmoluchowskiReactionModel",
00136 FatalErrorInArgument,exceptionDescription);
00137 }
00138
00139 G4double __preStepSeparation = (__preStepPositionA - __preStepPositionB).mag();
00140
00141 G4double __probabiltyOfEncounter = std::exp(-(__preStepSeparation - __R)*(__postStepSeparation - __R)
00142 / (__D* (__trackB.GetStep()->GetDeltaTime())));
00143 G4double __selectedPOE = G4UniformRand();
00144
00145 if(__selectedPOE<=__probabiltyOfEncounter) return true;
00146 }
00147
00148 return false ;
00149 }