Geant4_10
G4DNASmoluchowskiReactionModel.cc
Go to the documentation of this file.
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 // $Id: G4DNASmoluchowskiReactionModel.cc 70171 2013-05-24 13:34:18Z gcosmo $
27 //
29 #include "Randomize.hh"
30 #include "G4Track.hh"
32 #include "G4UnitsTable.hh"
33 
35 {
36  fReactionData = 0 ;
37 }
38 
40  G4VDNAReactionModel(__right)
41 {
42  fReactionData = 0 ;
43 }
44 
45 G4DNASmoluchowskiReactionModel& G4DNASmoluchowskiReactionModel::operator=(const G4DNASmoluchowskiReactionModel& right)
46 {
47  if(this == &right) return *this;
48  fReactionData = 0;
49  return *this;
50 }
51 
53 {
54  fReactionData = 0 ;
55 }
56 
58 {
59  fReactionData = fReactionTable->GetReactionData(__molecule);
60 }
61 
63 {
64  fReactionData = fReactionTable->GetReactionData(__molecule);
65 }
66 
68  const G4Molecule* mol2)
69 {
70  G4double __output = fReactionTable -> GetReactionData(mol1,mol2)->GetReducedReactionRadius();
71  return __output ;
72 }
73 
75 {
76  G4double __output = (*fReactionData)[__i] -> GetReducedReactionRadius();
77  return __output ;
78 }
79 
81  const G4Track& __trackB,
82  const G4double __R,
83  G4double& __r,
84  const G4bool __alongStepReaction)
85 {
86  G4double __postStepSeparation = 0;
87  bool __do_break = false ;
88  G4double __R2 = __R*__R ;
89  int k = 0 ;
90 
91  for(; k < 3 ; k++)
92  {
93  __postStepSeparation += std::pow(__trackA.GetPosition()[k] - __trackB.GetPosition()[k],2);
94 
95  if(__postStepSeparation > __R2)
96  {
97  __do_break = true ;
98  break ;
99  }
100  }
101 
102  if(__do_break == false)
103  {
104  // The loop was not break
105  // => __r^2 < __R^2
106  __r = std::sqrt(__postStepSeparation);
107  return true;
108  }
109  else if(__alongStepReaction == true)
110  {
111  //G4cout << "alongStepReaction==true" << G4endl;
112  //Along step cheack and
113  // the loop has break
114 
115  // Continue loop
116  for(; k < 3 ; k++)
117  {
118  __postStepSeparation += std::pow(__trackA.GetPosition()[k] - __trackB.GetPosition()[k],2);
119  }
120  // Use Green approach : the Brownian bridge
121  __r = (__postStepSeparation = std::sqrt(__postStepSeparation) );
122 
123  G4Molecule* __moleculeA = GetMolecule(__trackA);
124  G4Molecule* __moleculeB = GetMolecule(__trackB);
125 
126  G4double __D = __moleculeA->GetDiffusionCoefficient() + __moleculeB->GetDiffusionCoefficient();
127 
128  G4ThreeVector __preStepPositionA = __trackA.GetStep()->GetPreStepPoint() ->GetPosition();
129  G4ThreeVector __preStepPositionB = __trackB.GetStep()->GetPreStepPoint() ->GetPosition();
130 
131  if(__preStepPositionA == __trackA.GetPosition())
132  {
133  G4ExceptionDescription exceptionDescription ;
134  exceptionDescription << "The molecule : " << __moleculeA->GetName();
135  exceptionDescription << " with track ID :" << __trackA.GetTrackID();
136  exceptionDescription << " did not move since the previous step." << G4endl;
137  exceptionDescription << "Current position : " << G4BestUnit(__trackA.GetPosition(),"Length") << G4endl;
138  exceptionDescription << "Previous position : " << G4BestUnit(__preStepPositionA,"Length") << G4endl;
139  G4Exception("G4DNASmoluchowskiReactionModel::FindReaction","G4DNASmoluchowskiReactionModel",
140  FatalErrorInArgument,exceptionDescription);
141  }
142 
143  G4double __preStepSeparation = (__preStepPositionA - __preStepPositionB).mag();
144 
145  G4double __probabiltyOfEncounter = std::exp(-(__preStepSeparation - __R)*(__postStepSeparation - __R)
146  / (__D* (__trackB.GetStep()->GetDeltaTime())));
147  G4double __selectedPOE = G4UniformRand();
148 
149  if(__selectedPOE<=__probabiltyOfEncounter) return true;
150  }
151 
152  return false ;
153 }
virtual G4bool FindReaction(const G4Track &, const G4Track &, const G4double, G4double &, const G4bool)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
const G4ThreeVector & GetPosition() const
G4double GetDiffusionCoefficient() const
Definition: G4Molecule.cc:389
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
const G4Step * GetStep() const
int G4int
Definition: G4Types.hh:78
virtual void Initialise(const G4Molecule *, const G4Track &)
const G4DNAMolecularReactionTable * fReactionTable
G4StepPoint * GetPreStepPoint() const
const G4String & GetName() const
Definition: G4Molecule.cc:243
#define G4UniformRand()
Definition: Randomize.hh:87
const G4ThreeVector & GetPosition() const
bool G4bool
Definition: G4Types.hh:79
virtual void InitialiseToPrint(const G4Molecule *)
G4double GetDeltaTime() const
G4int GetTrackID() const
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:67
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const G4DNAMolecularReactionData * GetReactionData(const G4Molecule *, const G4Molecule *) const
virtual G4double GetReactionRadius(const G4Molecule *, const G4Molecule *)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76