Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
34 #include <G4MoleculeFinder.hh>
35 #include "G4PhysicalConstants.hh"
36 #include "G4Electron_aq.hh"
37 #include "G4MoleculeTable.hh"
39 #include "G4H2.hh"
40 #include "G4H2O.hh"
42 #include "G4MoleculeTable.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4VMoleculeCounter.hh"
47 #include "G4Exp.hh"
48 
49 static double onsager_constant = e_squared / (4. * pi * epsilon0 * k_Boltzmann);
50 
51 //------------------------------------------------------------------------------
52 
53 // Parameterisation of dielectric constant vs temperature and density
54 
55 double Y(double density)
56 {
57  return 1. / (1. + 0.0012 / (density * density));
58 }
59 
60 double A(double temperature)
61 {
62  double temp_inverse = 1 / temperature;
63  return 0.7017
64  + 642.0 * temp_inverse
65  - 1.167e5 * temp_inverse * temp_inverse
66  + 9.190e6 * temp_inverse * temp_inverse * temp_inverse;
67 }
68 
69 double B(double temperature)
70 {
71  double temp_inverse = 1 / temperature;
72  return -2.71
73  + 275.4 * temp_inverse
74  + 0.3245e5 * temp_inverse * temp_inverse;
75 }
76 
77 double S(double temp)
78 {
79  double temp_inverse = 1 / temp;
80 
81  return 1.667
82  - 11.41 * temp_inverse
83  - 35260.0 * temp_inverse * temp_inverse;
84 }
85 
86 double C(double temp)
87 {
88  return A(temp) - B(temp) - 3;
89 }
90 
91 double D(double temp)
92 {
93  return B(temp) + 3;
94 }
95 
96 double epsilon(double density, double temperature)
97 {
98  return 1 + G4Exp(std::log(10.)*
99  (Y(density) *
100  (C(temperature) + (S(temperature) - 1)*std::log(density)/std::log(10.))
101  + D(temperature) + std::log(density)/std::log(10.)));
102 }
103 
104 //------------------------------------------------------------------------------
105 
107  G4VITRestDiscreteProcess("G4DNAElectronHoleRecombination", fElectromagnetic)
108 {
109  Create();
110 // G4cout << epsilon(1.0095, 298.) << G4endl;
111 // G4cout << epsilon(1., 293.15) << G4endl;
112 // G4cout << epsilon(0.9277, 423.) << G4endl;
113 // G4cout << epsilon(0.816, 523.) << G4endl;
114 // G4cout << epsilon(0.6, 623) << G4endl;
115 }
116 
118 {
119 }
120 
122 {
123  pParticleChange = &fParticleChange;
124  enableAtRestDoIt = true;
125  enableAlongStepDoIt = false;
126  enablePostStepDoIt = true;
127 
128  SetProcessSubType(60);
129 
131  // ie G4DNAElectronHoleRecombination uses a state class
132  // inheriting from G4ProcessState
133 
134  fIsInitialized = false;
135  fProposesTimeStep = true;
136  fpMoleculeDensity = 0;
137 
138  verboseLevel = 0;
139 }
140 
141 //______________________________________________________________________________
142 
145  const G4Step& /*stepData*/)
146 {
147  fParticleChange.Initialize(track);
150  MakeReaction(track);
151  return &fParticleChange;
152 }
153 
154 //______________________________________________________________________________
155 
157 {
159  G4VITProcess::fpState.reset(new State());
161 }
162 
163 //______________________________________________________________________________
164 
165 void G4DNAElectronHoleRecombination::MakeReaction(const G4Track& track)
166 {
167  fParticleChange.Initialize(track);
168  State* state = fpState->GetState<State>();
169  double random = state->fSampleProba;
170  std::vector<ReactionProfile>& reactants = state->fReactants;
171 
172  G4Track* selected_reactant = 0;
173 
174  for(size_t i = 0; i < reactants.size(); ++i)
175  {
176  if(reactants[i].fElectron->GetTrackStatus() != fAlive) continue;
177  if(reactants[i].fProbability > random)
178  {
179  selected_reactant = reactants[i].fElectron;
180  }
181  break;
182  }
183 
184  // G4cout << "MakeReaction with charge ="
185  // << GetMolecule(track)->GetCharge() << G4endl;
186 
187  if(selected_reactant)
188  {
189  // G4cout << " Will react with TID = " << selected_reactant->GetTrackID()
190  // << G4endl;
191 
194  RemoveAMoleculeAtTime(GetMolecule(track)->GetMolecularConfiguration(),
195  track.GetGlobalTime(),
196  &(track.GetPosition()));
197  GetMolecule(track)->ChangeConfigurationToLabel("H2Ovib");
198 
201  AddAMoleculeAtTime(GetMolecule(track)->GetMolecularConfiguration(),
202  track.GetGlobalTime(),
203  &(track.GetPosition()));
204 
205  // fParticleChange.ProposeTrackStatus(fStopAndKill);
206  fParticleChange.ProposeTrackStatus(fStopButAlive);
207 
208  selected_reactant->SetTrackStatus(fStopAndKill);
209  // G4TrackList::Pop(selected_reactant);
210  // G4ITTrackHolder::Instance()->PushToKill(selected_reactant);
211 
212  }
213  else
214  {
215  fParticleChange.ProposeTrackStatus(fStopButAlive);
216  }
217 }
218 
219 //______________________________________________________________________________
220 
221 G4bool G4DNAElectronHoleRecombination::FindReactant(const G4Track& track)
222 {
223  if(GetMolecule(track)->GetCharge() <= 0)
224  {
225  // G4cout << "La charge est negative ou nulle !! " << G4endl;
226  return false;
227  }
228 
229  const std::vector<double>* densityTable =
231 
232  double temperature = track.GetMaterial()->GetTemperature();
233  double density = (*densityTable)[track.GetMaterial()->GetIndex()] /
234  ( g/(1e-2*m*1e-2*m*1e-2*m) );
235  double eps = epsilon(density, temperature);
236 
237  // G4cout << " temperature = " << temperature << G4endl;
238  // G4cout << " density = " << density << G4endl;
239  // G4cout << " eps = " << eps << G4endl;
240 
241  double onsager_radius = onsager_constant * 1. / (temperature * eps);
242 
244 
247  e_aq.GetMoleculeID(),
248  10. * onsager_radius);
249 
250  // double distance = -1.;
251  // double probability = -1.;
252 
253  if(results == 0 || results->GetSize() == 0)
254  {
255  // G4cout << "rien trouve a moins de 10 rc" << G4endl;
256  return false;
257  }
258 
259  results->Sort();
260 
261  State* state = fpState->GetState<State>();
262  std::vector<ReactionProfile>& reactants = state->fReactants;
263  state->fSampleProba = G4UniformRand();
264 
265  reactants.resize(results->GetSize());
266 
267  for(size_t i = 0; results->End() == false; results->Next(), ++i)
268  {
269  reactants[i].fElectron = results->GetItem<G4IT>()->GetTrack();
270  reactants[i].fDistance = std::sqrt(results->GetDistanceSqr());
271 
272  if(reactants[i].fDistance != 0)
273  {
274  reactants[i].fProbability = 1.
275  - G4Exp(-onsager_radius / reactants[i].fDistance);
276  }
277  else
278  {
279  reactants[i].fProbability = 1.;
280  }
281 
282  // G4cout << "dis = "
283  // << reactants[i].fDistance << " prob = " << reactants[i].fProbability << G4endl;
284  }
285 
286  if(results->GetSize() != 0 && reactants.empty())
287  {
288  G4cout << "Size is = " << results->GetSize() << G4endl;
289  abort();
290  }
291 
292  if(reactants.empty()) return false;
293 
294  // G4cout << " reactants[0].fDistance =" << reactants[0].fDistance
295  // << " onsager_radius = " << onsager_radius << "\t";
296  //
297  // G4cout << " reactants[0].fProbability =" << reactants[0].fProbability
298  // << "state->fSampleProba = " << state->fSampleProba << "\t";
299 
300  if(reactants[0].fProbability > state->fSampleProba) return true;
301  return false;
302 }
303 
304 //______________________________________________________________________________
305 
306 G4bool
309 {
311  {
313  ->GetMoleculeDefinition("H2O", false);
314 
315  if(H2O) // if this condition does not hold => process cannot be applied
316  {
318  G4H2O::Definition()->NewConfiguration("H2Ovib");
319 
320  assert(vib != 0);
321 
323  ->GetConfiguration("H2", false);
325  ->GetConfiguration("OH", false);
327  ->GetConfiguration("H", false);
328 
329  double probaRemaining = 1.;
330 
331  if(OH || H2)
332  {
334  new G4MolecularDissociationChannel("H2Ovib_DissociativeDecay1");
335  if(H2)
336  {
337  diss2->AddProduct(H2);
338  }
339  if(OH)
340  {
341  diss2->AddProduct(OH);
342  diss2->AddProduct(OH);
343  }
344 
345  double proba = 0.15;
346  diss2->SetProbability(proba);
347  probaRemaining -= proba;
349  B1A1_DissociationDecay);
350  G4H2O::Definition()->AddDecayChannel(vib, diss2);
351  }
352 
353  if(OH || H)
354  {
356  new G4MolecularDissociationChannel("H2Ovib_DissociativeDecay2");
357  if(OH)
358  {
359  diss3->AddProduct(OH);
360  }
361  if(H)
362  {
363  diss3->AddProduct(H);
364  }
365  double proba = 0.55;
366  diss3->SetProbability(proba);
367  probaRemaining -= proba;
369  A1B1_DissociationDecay);
370  G4H2O::Definition()->AddDecayChannel(vib, diss3);
371  }
372 
374  new G4MolecularDissociationChannel("H2Ovib_NonDissociative");
375  channel1->SetProbability(probaRemaining);
376  G4H2O::Definition()->AddDecayChannel(vib, channel1);
377  }
378  }
379 
380  if(particle.GetParticleName() == "H2O") return true;
381  return false;
382 }
383 
384 //______________________________________________________________________________
385 
387  G4double,
389 {
390  // G4cout << "track ID = " << track.GetTrackID()
391  // << " G4DNAElectronHoleRecombination --> ";
392  if(FindReactant(track))
393  {
394  // G4cout << " SELECTED " ;
395  // G4cout << G4endl;
396  return 0;
397  }
398 
399  // G4cout << " NOT SELECTED " << G4endl;
400  return DBL_MAX;
401 }
402 
403 //______________________________________________________________________________
404 
407 {
408 // G4cout << "track ID = " << track.GetTrackID()
409 // << " G4DNAElectronHoleRecombination --> ";
410  if(FindReactant(track))
411  {
412 // G4cout << " SELECTED " ;
413 // G4cout << G4endl;
414  return 0;
415  }
416 
417 // G4cout << " NOT SELECTED " << G4endl;
418  return DBL_MAX;
419 }
420 
422  const G4Step& step)
423 {
424  return AtRestDoIt(track, step);
425 }
void SetTrackStatus(const G4TrackStatus aTrackStatus)
static double onsager_constant
Definition: G4IT.hh:88
const std::vector< double > * GetDensityTableFor(const G4Material *) const
static G4Electron_aq * Definition()
virtual void ClearNumberOfInteractionLengthLeft()
double Y(double density)
Identical to G4VRestDiscreteProcess with dependency from G4VITProcess.
G4MoleculeDefinition * GetMoleculeDefinition(const G4String &, bool mustExist=true)
G4int verboseLevel
Definition: G4VProcess.hh:368
static G4H2O * Definition()
Definition: G4H2O.cc:46
double S(double temp)
void AddDecayChannel(const G4MolecularConfiguration *molConf, const G4MolecularDissociationChannel *channel)
size_t GetIndex() const
Definition: G4Material.hh:262
const G4ThreeVector & GetPosition() const
#define State(X)
static G4ITFinder * Instance()
static const G4double eps
double C(double temp)
double B(double temperature)
G4KDTreeResultHandle FindNearestInRange(const T *point, int key, G4double)
virtual G4bool IsApplicable(const G4ParticleDefinition &)
G4bool enableAtRestDoIt
Definition: G4VProcess.hh:350
const G4String & GetParticleName() const
void SetInstantiateProcessState(G4bool flag)
void ChangeConfigurationToLabel(const G4String &label)
Definition: G4Molecule.cc:611
virtual void ClearInteractionTimeLeft()
G4MolecularConfiguration * NewConfiguration(const G4String &excitedStateLabel)
virtual void StartTracking(G4Track *)
Definition: G4VProcess.cc:101
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
static constexpr double m
Definition: G4SIunits.hh:129
float k_Boltzmann
Definition: hepunit.py:299
void AddProduct(const G4Molecule *, G4double=0)
bool G4bool
Definition: G4Types.hh:79
virtual G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
static G4MoleculeTable * Instance()
virtual void StartTracking(G4Track *)
Definition: G4VITProcess.cc:86
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
G4shared_ptr< G4ProcessState > fpState
Definition: G4Step.hh:76
G4bool enablePostStepDoIt
Definition: G4VProcess.hh:352
G4double GetGlobalTime() const
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:76
static G4VMoleculeCounter * Instance()
G4Material * GetMaterial() const
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
virtual void Initialize(const G4Track &)
static G4DNAMolecularMaterial * Instance()
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
virtual G4double GetMeanLifeTime(const G4Track &aTrack, G4ForceCondition *condition)
double D(double temp)
G4bool fProposesTimeStep
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
G4bool enableAlongStepDoIt
Definition: G4VProcess.hh:351
G4bool IsMasterThread()
Definition: G4Threading.cc:146
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
static G4bool InUse()
G4ForceCondition
#define DBL_MAX
Definition: templates.hh:83
double epsilon(double density, double temperature)
G4MolecularConfiguration * GetConfiguration(const G4String &, bool mustExist=true)