Geant4-11
G4DNAElectronHoleRecombination.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/*
27 * G4DNAElectronHoleRecombination.cc
28 *
29 * Created on: Jun 17, 2015
30 * Author: mkaramit
31 *
32 */
33
36#include "G4MoleculeFinder.hh"
37#include "G4Molecule.hh"
39#include "G4Electron_aq.hh"
40#include "G4H2O.hh"
43#include "G4SystemOfUnits.hh"
44#include "G4VMoleculeCounter.hh"
45#include "G4Exp.hh"
47
49
50//------------------------------------------------------------------------------
51
52// Parameterisation of dielectric constant vs temperature and density
53
55{
56 return 1. / (1. + 0.0012 / (density * density));
57}
58
59G4double A(G4double temperature)
60{
61 G4double temp_inverse = 1 / temperature;
62 return 0.7017
63 + 642.0 * temp_inverse
64 - 1.167e5 * temp_inverse * temp_inverse
65 + 9.190e6 * temp_inverse * temp_inverse * temp_inverse;
66}
67
68G4double B(G4double temperature)
69{
70 G4double temp_inverse = 1 / temperature;
71 return -2.71
72 + 275.4 * temp_inverse
73 + 0.3245e5 * temp_inverse * temp_inverse;
74}
75
77{
78 G4double temp_inverse = 1 / temp;
79
80 return 1.667
81 - 11.41 * temp_inverse
82 - 35260.0 * temp_inverse * temp_inverse;
83}
84
86{
87 return A(temp) - B(temp) - 3;
88}
89
91{
92 return B(temp) + 3;
93}
94
95G4double epsilon(G4double density, G4double temperature)
96{
97 return 1 + G4Exp(std::log(10.) *
98 (Y(density) *
99 (C(temperature) + (S(temperature) - 1) * std::log(density) / std::log(10.))
100 + D(temperature) + std::log(density) / std::log(10.)));
101}
102
103//------------------------------------------------------------------------------
104
106 : G4VITRestDiscreteProcess("G4DNAElectronHoleRecombination",
108{
109 Create();
110}
111
113
115{
117 enableAtRestDoIt = true;
118 enableAlongStepDoIt = false;
119 enablePostStepDoIt = true;
120
122
124 // ie G4DNAElectronHoleRecombination uses a state class
125 // inheriting from G4ProcessState
126
127 fIsInitialized = false;
128 fProposesTimeStep = true;
129 fpMoleculeDensity = nullptr;
130
131 verboseLevel = 0;
132}
133
134//______________________________________________________________________________
135
138 const G4Step& /*stepData*/)
139{
143 MakeReaction(track);
144 return &fParticleChange;
145}
146
147//______________________________________________________________________________
148
150{
152 G4VITProcess::fpState.reset(new State());
154}
155
156//______________________________________________________________________________
157
159{
161 auto pState = fpState->GetState<State>();
162 G4double random = pState->fSampleProba;
163 std::vector<ReactantInfo>& reactants = pState->fReactants;
164
165 G4Track* pSelectedReactant = nullptr;
166
167 for (const auto& reactantInfo : reactants)
168 {
169 if (reactantInfo.fElectron->GetTrackStatus() != fAlive)
170 {
171 continue;
172 }
173 if (reactantInfo.fProbability > random)
174 {
175 pSelectedReactant = reactantInfo.fElectron;
176 }
177 break;
178 }
179
180 if (pSelectedReactant)
181 {
183 {
185 RemoveAMoleculeAtTime(GetMolecule(track)->GetMolecularConfiguration(),
186 track.GetGlobalTime(),
187 &(track.GetPosition()));
188 }
189 GetMolecule(track)->ChangeConfigurationToLabel("H2Ovib");
190
192 {
194 AddAMoleculeAtTime(GetMolecule(track)->GetMolecularConfiguration(),
195 track.GetGlobalTime(),
196 &(track.GetPosition()));
197 }
198
199 // fParticleChange.ProposeTrackStatus(fStopAndKill);
201
202 pSelectedReactant->SetTrackStatus(fStopAndKill);
203 // G4TrackList::Pop(pSelectedReactant);
204 // G4ITTrackHolder::Instance()->PushToKill(pSelectedReactant);
205 }
206 else
207 {
209 }
210}
211
212//______________________________________________________________________________
213
215{
216 if (GetMolecule(track)->GetCharge() <= 0)
217 {
218 return false;
219 }
220
221 const auto pDensityTable =
223
224 G4double temperature = track.GetMaterial()->GetTemperature();
225 G4double density = (*pDensityTable)[track.GetMaterial()->GetIndex()] / (g / (1e-2 * m * 1e-2 * m * 1e-2 * m));
226 G4double eps = epsilon(density, temperature);
227
228 G4double onsagerRadius = onsager_constant * 1. / (temperature * eps);
229
231
232 auto pState = fpState->GetState<State>();
233 std::vector<ReactantInfo>& reactants = pState->fReactants;
234 pState->fSampleProba = G4UniformRand();
235
236 //Updated : Hoang Tran: added Octree finder
237 if(G4ChemicalMoleculeFinder::Instance()->IsOctreeUsed())
238 {
239 if(!G4ChemicalMoleculeFinder::Instance()->IsOctreeBuilt())
240 {
243 }
244 std::vector<std::pair<G4TrackList::iterator,G4double>> resultIndices;
245 resultIndices.clear();
246
248 FindNearest(track,
249 e_aq.GetMoleculeID(),
250 10. * onsagerRadius,
251 resultIndices,
252 true);
253
254 if(resultIndices.empty())
255 {
256 return false;
257 }
258 reactants.resize(resultIndices.size());
259 unsigned int i = 0;
260 for(auto& it : resultIndices)
261 {
262 reactants[i].fElectron = *(std::get<0>(it));
263 reactants[i].fDistance = (reactants[i].fElectron->GetPosition() -
264 track.GetPosition()).mag();
265 if (reactants[i].fDistance != 0)
266 {
267 reactants[i].fProbability = 1. - G4Exp(-onsagerRadius /
268 reactants[i].fDistance);
269 }
270 else
271 {
272 reactants[i].fProbability = 1.;
273 }
274 i++;
275 }
276 }
277 else
278 {
281 e_aq.GetMoleculeID(),
282 10. * onsagerRadius);
283
284 if (results == 0 || results->GetSize() == 0)
285 {
286 return false;
287 }
288
289 results->Sort();
290 reactants.resize(results->GetSize());
291
292 for (size_t i = 0; !results->End(); results->Next(), ++i)
293 {
294 reactants[i].fElectron = results->GetItem<G4IT>()->GetTrack();
295 reactants[i].fDistance = std::sqrt(results->GetDistanceSqr());
296
297 if (reactants[i].fDistance != 0)
298 {
299 reactants[i].fProbability = 1. - G4Exp(-onsagerRadius / reactants[i].fDistance);
300 }
301 else
302 {
303 reactants[i].fProbability = 1.;
304 }
305 }
306 }
307 return reactants.empty() ? false : reactants[0].fProbability > pState->fSampleProba;
308}
309
310//______________________________________________________________________________
311
312G4bool
314IsApplicable(const G4ParticleDefinition& particle)
315{
316 if (&particle != G4H2O::DefinitionIfExists())
317 {
318 return false;
319 }
320
321 return true;
322}
323
324//______________________________________________________________________________
325
327 G4double,
329{
330 if (FindReactant(track))
331 {
332 return 0.;
333 }
334
335 return DBL_MAX;
336}
337
338//______________________________________________________________________________
339
342{
343 if (FindReactant(track))
344 {
345 return 0.;
346 }
347 return DBL_MAX;
348}
349
351 const G4Step& step)
352{
353 return AtRestDoIt(track, step);
354}
#define BuildChemicalMoleculeFinder()
G4double epsilon(G4double density, G4double temperature)
G4double A(G4double temperature)
G4double C(G4double temp)
G4double S(G4double temp)
G4double B(G4double temperature)
G4double D(G4double temp)
G4double Y(G4double density)
static G4double onsager_constant
static const G4double eps
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
G4ForceCondition
#define State(X)
@ fLowEnergyTransportation
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:76
@ fElectromagnetic
static constexpr double m
Definition: G4SIunits.hh:109
static constexpr double g
Definition: G4SIunits.hh:168
static constexpr double pi
Definition: G4SIunits.hh:55
@ fAlive
@ fStopAndKill
@ fStopButAlive
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
#define G4UniformRand()
Definition: Randomize.hh:52
G4double GetMeanLifeTime(const G4Track &aTrack, G4ForceCondition *condition) override
G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &) override
G4bool IsApplicable(const G4ParticleDefinition &) override
const std::vector< double > * fpMoleculeDensity
G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &) override
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) override
static G4DNAMolecularMaterial * Instance()
const std::vector< G4double > * GetDensityTableFor(const G4Material *) const
Retrieve a table of volumetric mass densities (mass per unit volume) in the G4 unit system for chosen...
static G4Electron_aq * Definition()
static G4H2O * DefinitionIfExists()
Definition: G4H2O.cc:80
G4KDTreeResultHandle FindNearestInRange(const T *point, int key, G4double)
static G4ITFinder * Instance()
Definition: G4IT.hh:88
G4double GetTemperature() const
Definition: G4Material.hh:178
size_t GetIndex() const
Definition: G4Material.hh:256
void ChangeConfigurationToLabel(const G4String &label)
Definition: G4Molecule.cc:546
G4int GetMoleculeID() const
Definition: G4Molecule.cc:467
static G4OctreeFinder * Instance()
void SetOctreeBuilt(G4bool used)
virtual void Initialize(const G4Track &)
Definition: G4Step.hh:62
void SetTrackStatus(const G4TrackStatus aTrackStatus)
const G4ThreeVector & GetPosition() const
G4double GetGlobalTime() const
G4Material * GetMaterial() const
virtual void StartTracking(G4Track *)
Definition: G4VITProcess.cc:85
void SetInstantiateProcessState(G4bool flag)
G4shared_ptr< G4ProcessState > fpState
virtual void ClearInteractionTimeLeft()
virtual void ClearNumberOfInteractionLengthLeft()
G4bool fProposesTimeStep
Identical to G4VRestDiscreteProcess with dependency from G4VITProcess.
static G4VMoleculeCounter * Instance()
void ProposeTrackStatus(G4TrackStatus status)
G4bool enableAtRestDoIt
Definition: G4VProcess.hh:359
G4int verboseLevel
Definition: G4VProcess.hh:356
G4bool enableAlongStepDoIt
Definition: G4VProcess.hh:360
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:406
virtual void StartTracking(G4Track *)
Definition: G4VProcess.cc:87
G4bool enablePostStepDoIt
Definition: G4VProcess.hh:361
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:321
float k_Boltzmann
Definition: hepunit.py:298
#define DBL_MAX
Definition: templates.hh:62