Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EmDNAChemistry.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 #include "G4EmDNAChemistry.hh"
27 
28 #include "G4PhysicalConstants.hh"
29 #include "G4SystemOfUnits.hh"
30 
32 #include "G4DNAChemistryManager.hh"
34 #include "G4ProcessManager.hh"
35 
37 
38 // *** Processes and models for Geant4-DNA
39 
41 
42 #include "G4DNAAttachment.hh"
43 #include "G4DNAVibExcitation.hh"
44 
45 #include "G4DNAElastic.hh"
49 
54 #include "G4VDNAReactionModel.hh"
56 
58 
59 // particles
60 
61 #include "G4Electron.hh"
62 #include "G4Proton.hh"
63 #include "G4GenericIon.hh"
64 
65 #include "G4MoleculeTable.hh"
66 #include "G4H2O.hh"
67 #include "G4H2.hh"
68 #include "G4Hydrogen.hh"
69 #include "G4OH.hh"
70 #include "G4H3O.hh"
71 #include "G4Electron_aq.hh"
72 #include "G4H2O2.hh"
73 
74 #include "G4PhysicsListHelper.hh"
75 #include "G4BuilderType.hh"
76 
77 /****/
79 #include "G4ProcessVector.hh"
80 #include "G4ProcessTable.hh"
83 /****/
84 
85 // factory
87 
89 
90 #include "G4Threading.hh"
91 
94 {
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 
101 {
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
107 {
108  //-----------------------------------
109  G4Electron::Definition(); // safety
110 
111  //-----------------------------------
112  // Create the definition
120 
121  //____________________________________________________________________________
122 
125  CreateConfiguration("OHm", // just a tag to store and retrieve from
126  // G4MoleculeTable
128  -1, // charge
129  5.0e-9 * (m2 / s));
130  OHm->SetMass(17.0079 * g / Avogadro * c_squared);
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 
143 {
144  //-----------------------------------
145  //Get the molecular configuration
158 
159  //-------------------------------------
160  //Define the decay channels
164 
166  *(water->GetGroundStateElectronOccupancy()));
167 
169  // EXCITATIONS //
171  G4DNAWaterExcitationStructure waterExcitation;
172  //--------------------------------------------------------
173  //---------------Excitation on the fifth layer------------
174 
175  decCh1 = new G4MolecularDissociationChannel("A^1B_1_Relaxation");
176  decCh2 = new G4MolecularDissociationChannel("A^1B_1_DissociativeDecay");
177  //Decay 1 : OH + H
178  decCh1->SetEnergy(waterExcitation.ExcitationEnergy(0));
179  decCh1->SetProbability(0.35);
180  decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::NoDisplacement);
181 
182  decCh2->AddProduct(OH);
183  decCh2->AddProduct(H);
184  decCh2->SetProbability(0.65);
185  decCh2->SetDisplacementType(
186  G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
187 
188 // water->AddExcitedState("A^1B_1");
189  occ->RemoveElectron(4, 1); // this is the transition form ground state to
190  occ->AddElectron(5, 1); // the first unoccupied orbital: A^1B_1
191 
192  water->NewConfigurationWithElectronOccupancy("A^1B_1", *occ);
193  water->AddDecayChannel("A^1B_1", decCh1);
194  water->AddDecayChannel("A^1B_1", decCh2);
195 
196  //--------------------------------------------------------
197  //---------------Excitation on the fourth layer-----------
198  decCh1 = new G4MolecularDissociationChannel("B^1A_1_Relaxation_Channel");
199  decCh2 = new G4MolecularDissociationChannel("B^1A_1_DissociativeDecay");
201  "B^1A_1_AutoIonisation_Channel");
202 
203  //Decay 1 : energy
204  decCh1->SetEnergy(waterExcitation.ExcitationEnergy(1));
205  decCh1->SetProbability(0.3);
206 
207  //Decay 2 : 2OH + H_2
208  decCh2->AddProduct(H2);
209  decCh2->AddProduct(OH);
210  decCh2->AddProduct(OH);
211  decCh2->SetProbability(0.15);
212  decCh2->SetDisplacementType(
213  G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay);
214 
215  //Decay 3 : OH + H_3Op + e_aq
216  decCh3->AddProduct(OH);
217  decCh3->AddProduct(H3O);
218  decCh3->AddProduct(e_aq);
219  decCh3->SetProbability(0.55);
220  decCh3->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
221 
222  *occ = *(water->GetGroundStateElectronOccupancy());
223  occ->RemoveElectron(3); // this is the transition form ground state to
224  occ->AddElectron(5, 1); // the first unoccupied orbital: B^1A_1
225 
226  water->NewConfigurationWithElectronOccupancy("B^1A_1", *occ);
227  water->AddDecayChannel("B^1A_1", decCh1);
228  water->AddDecayChannel("B^1A_1", decCh2);
229  water->AddDecayChannel("B^1A_1", decCh3);
230 
231  //-------------------------------------------------------
232  //-------------------Excitation of 3rd layer-----------------
233  decCh1 = new G4MolecularDissociationChannel(
234  "Excitation3rdLayer_AutoIonisation_Channel");
235  decCh2 = new G4MolecularDissociationChannel(
236  "Excitation3rdLayer_Relaxation_Channel");
237 
238  //Decay channel 1 : : OH + H_3Op + e_aq
239  decCh1->AddProduct(OH);
240  decCh1->AddProduct(H3O);
241  decCh1->AddProduct(e_aq);
242 
243  decCh1->SetProbability(0.5);
244  decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
245 
246  //Decay channel 2 : energy
247  decCh2->SetEnergy(waterExcitation.ExcitationEnergy(2));
248  decCh2->SetProbability(0.5);
249 
250  //Electronic configuration of this decay
251  *occ = *(water->GetGroundStateElectronOccupancy());
252  occ->RemoveElectron(2, 1);
253  occ->AddElectron(5, 1);
254 
255  //Configure the water molecule
256  water->NewConfigurationWithElectronOccupancy("Excitation3rdLayer", *occ);
257  water->AddDecayChannel("Excitation3rdLayer", decCh1);
258  water->AddDecayChannel("Excitation3rdLayer", decCh2);
259 
260  //-------------------------------------------------------
261  //-------------------Excitation of 2nd layer-----------------
262  decCh1 = new G4MolecularDissociationChannel(
263  "Excitation2ndLayer_AutoIonisation_Channel");
264  decCh2 = new G4MolecularDissociationChannel(
265  "Excitation2ndLayer_Relaxation_Channel");
266 
267  //Decay Channel 1 : : OH + H_3Op + e_aq
268  decCh1->AddProduct(OH);
269  decCh1->AddProduct(H3O);
270  decCh1->AddProduct(e_aq);
271 
272  decCh1->SetProbability(0.5);
273  decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
274 
275  //Decay channel 2 : energy
276  decCh2->SetEnergy(waterExcitation.ExcitationEnergy(3));
277  decCh2->SetProbability(0.5);
278 
279  *occ = *(water->GetGroundStateElectronOccupancy());
280  occ->RemoveElectron(1, 1);
281  occ->AddElectron(5, 1);
282 
283  water->NewConfigurationWithElectronOccupancy("Excitation2ndLayer", *occ);
284  water->AddDecayChannel("Excitation2ndLayer", decCh1);
285  water->AddDecayChannel("Excitation2ndLayer", decCh2);
286 
287  //-------------------------------------------------------
288  //-------------------Excitation of 1st layer-----------------
289  decCh1 = new G4MolecularDissociationChannel(
290  "Excitation1stLayer_AutoIonisation_Channel");
291  decCh2 = new G4MolecularDissociationChannel(
292  "Excitation1stLayer_Relaxation_Channel");
293 
294  *occ = *(water->GetGroundStateElectronOccupancy());
295  occ->RemoveElectron(0, 1);
296  occ->AddElectron(5, 1);
297 
298  //Decay Channel 1 : : OH + H_3Op + e_aq
299  decCh1->AddProduct(OH);
300  decCh1->AddProduct(H3O);
301  decCh1->AddProduct(e_aq);
302  decCh1->SetProbability(0.5);
303  decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
304 
305  //Decay channel 2 : energy
306  decCh2->SetEnergy(waterExcitation.ExcitationEnergy(4));
307  decCh2->SetProbability(0.5);
308 
309  water->NewConfigurationWithElectronOccupancy("Excitation1stLayer", *occ);
310  water->AddDecayChannel("Excitation1stLayer", decCh1);
311  water->AddDecayChannel("Excitation1stLayer", decCh2);
312 
314  // IONISATION //
316  //--------------------------------------------------------
317  //------------------- Ionisation -------------------------
318 
319  decCh1 = new G4MolecularDissociationChannel("Ionisation_Channel");
320 
321  //Decay Channel 1 : : OH + H_3Op
322  decCh1->AddProduct(H3O);
323  decCh1->AddProduct(OH);
324  decCh1->SetProbability(1);
325  decCh1->SetDisplacementType(
326  G4DNAWaterDissociationDisplacer::Ionisation_DissociationDecay);
327 
328  *occ = *(water->GetGroundStateElectronOccupancy());
329  occ->RemoveElectron(4, 1);
330  // this is a ionized h2O with a hole in its last orbital
331  water->NewConfigurationWithElectronOccupancy("Ionisation5", *occ);
332  water->AddDecayChannel("Ionisation5",
333  decCh1);
334 
335  *occ = *(water->GetGroundStateElectronOccupancy());
336  occ->RemoveElectron(3, 1);
337  water->NewConfigurationWithElectronOccupancy("Ionisation4", *occ);
338  water->AddDecayChannel("Ionisation4",
339  new G4MolecularDissociationChannel(*decCh1));
340 
341  *occ = *(water->GetGroundStateElectronOccupancy());
342  occ->RemoveElectron(2, 1);
343  water->NewConfigurationWithElectronOccupancy("Ionisation3", *occ);
344  water->AddDecayChannel("Ionisation3",
345  new G4MolecularDissociationChannel(*decCh1));
346 
347  *occ = *(water->GetGroundStateElectronOccupancy());
348  occ->RemoveElectron(1, 1);
349  water->NewConfigurationWithElectronOccupancy("Ionisation2", *occ);
350  water->AddDecayChannel("Ionisation2",
351  new G4MolecularDissociationChannel(*decCh1));
352 
353  *occ = *(water->GetGroundStateElectronOccupancy());
354  occ->RemoveElectron(0, 1);
355  water->NewConfigurationWithElectronOccupancy("Ionisation1", *occ);
356  water->AddDecayChannel("Ionisation1",
357  new G4MolecularDissociationChannel(*decCh1));
358 
360  // Dissociative Attachment //
362  decCh1 = new G4MolecularDissociationChannel("DissociativeAttachment");
363 
364  //Decay 1 : 2OH + H_2
365  decCh1->AddProduct(H2);
366  decCh1->AddProduct(OHm);
367  decCh1->AddProduct(OH);
368  decCh1->SetProbability(1);
370  DissociativeAttachment);
371 
372  *occ = *(water->GetGroundStateElectronOccupancy());
373  occ->AddElectron(5, 1); // H_2O^-
374  water->NewConfigurationWithElectronOccupancy("DissociativeAttachment", *occ);
375  water->AddDecayChannel("DissociativeAttachment", decCh1);
376 
377  delete occ;
378 }
379 
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
381 
383  theReactionTable)
384 {
385  //-----------------------------------
386  //Get the molecular configuration
401 
402  //------------------------------------------------------------------
403  // e_aq + e_aq + 2H2O -> H2 + 2OH-
404  G4DNAMolecularReactionData* reactionData =
405  new G4DNAMolecularReactionData(0.5e10 * (1e-3 * m3 / (mole * s)), e_aq, e_aq);
406  reactionData->AddProduct(OHm);
407  reactionData->AddProduct(OHm);
408  reactionData->AddProduct(H2);
409  theReactionTable->SetReaction(reactionData);
410  //------------------------------------------------------------------
411  // e_aq + *OH -> OH-
412  reactionData = new G4DNAMolecularReactionData(
413  2.95e10 * (1e-3 * m3 / (mole * s)), e_aq, OH);
414  reactionData->AddProduct(OHm);
415  theReactionTable->SetReaction(reactionData);
416  //------------------------------------------------------------------
417  // e_aq + H* + H2O -> H2 + OH-
418  reactionData = new G4DNAMolecularReactionData(
419  2.65e10 * (1e-3 * m3 / (mole * s)), e_aq, H);
420  reactionData->AddProduct(OHm);
421  reactionData->AddProduct(H2);
422  theReactionTable->SetReaction(reactionData);
423  //------------------------------------------------------------------
424  // e_aq + H3O+ -> H* + H2O
425  reactionData = new G4DNAMolecularReactionData(
426  2.11e10 * (1e-3 * m3 / (mole * s)), e_aq, H3Op);
427  reactionData->AddProduct(H);
428  theReactionTable->SetReaction(reactionData);
429  //------------------------------------------------------------------
430  // e_aq + H2O2 -> OH- + *OH
431  reactionData = new G4DNAMolecularReactionData(
432  1.41e10 * (1e-3 * m3 / (mole * s)), e_aq, H2O2);
433  reactionData->AddProduct(OHm);
434  reactionData->AddProduct(OH);
435  theReactionTable->SetReaction(reactionData);
436  //------------------------------------------------------------------
437  // *OH + *OH -> H2O2
438  reactionData = new G4DNAMolecularReactionData(
439  0.44e10 * (1e-3 * m3 / (mole * s)), OH, OH);
440  reactionData->AddProduct(H2O2);
441  theReactionTable->SetReaction(reactionData);
442  //------------------------------------------------------------------
443  // *OH + *H -> H2O
444  theReactionTable->SetReaction(1.44e10 * (1e-3 * m3 / (mole * s)), OH, H);
445  //------------------------------------------------------------------
446  // *H + *H -> H2
447  reactionData = new G4DNAMolecularReactionData(
448  1.20e10 * (1e-3 * m3 / (mole * s)), H, H);
449  reactionData->AddProduct(H2);
450  theReactionTable->SetReaction(reactionData);
451  //------------------------------------------------------------------
452  // H3O+ + OH- -> 2H2O
453  theReactionTable->SetReaction(1.43e11 * (1e-3 * m3 / (mole * s)), H3Op, OHm);
454  //------------------------------------------------------------------
455 }
456 
457 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
458 
460 {
462 
463  //===============================================================
464  // Extend vibrational to low energy
465  // Anyway, solvation of electrons is taken into account from 7.4 eV
466  // So below this threshold, for now, no accurate modeling is done
467  //
468  G4VProcess* process =
470  FindProcess("e-_G4DNAVibExcitation", "e-");
471 
472  if (process)
473  {
474  G4DNAVibExcitation* vibExcitation = (G4DNAVibExcitation*) process;
475  G4VEmModel* model = vibExcitation->EmModel();
476  G4DNASancheExcitationModel* sancheExcitationMod =
477  dynamic_cast<G4DNASancheExcitationModel*>(model);
478  if(sancheExcitationMod)
479  {
480  sancheExcitationMod->ExtendLowEnergyLimit(0.025 * eV);
481  }
482  }
483 
484  //===============================================================
485  // *** Electron Solvatation ***
486  //
487  process =
489  FindProcess("e-_G4DNAElectronSolvation", "e-");
490 
491  if (process == 0)
492  {
493  ph->RegisterProcess(
494  new G4DNAElectronSolvation("e-_G4DNAElectronSolvation"),
496  }
497 
498  //===============================================================
499  // Define processes for molecules
500  //
501  G4MoleculeTable* theMoleculeTable = G4MoleculeTable::Instance();
503  theMoleculeTable->GetDefintionIterator();
504  iterator.reset();
505  while (iterator())
506  {
507  G4MoleculeDefinition* moleculeDef = iterator.value();
508 
509  if (moleculeDef != G4H2O::Definition())
510  {
511  // G4cout << "Brownian motion added for : "
512  // << moleculeDef->GetName() << G4endl;
514  // brown->SetVerboseLevel(4);
515  ph->RegisterProcess(brown, moleculeDef);
516  }
517  else
518  {
519  moleculeDef->GetProcessManager()
521  G4DNAMolecularDissociation* dissociationProcess =
522  new G4DNAMolecularDissociation("H2O_DNAMolecularDecay");
523  dissociationProcess->SetDecayDisplacer(
524  moleculeDef, new G4DNAWaterDissociationDisplacer);
525  dissociationProcess->SetVerboseLevel(1);
526 // ph->RegisterProcess(dissociationProcess, moleculeDef);
527 
528  moleculeDef->GetProcessManager()
529  ->AddRestProcess(dissociationProcess, 1);
530  }
531  /*
532  * Warning : end of particles and processes are needed by
533  * EM Physics builders
534  */
535  }
536 
538 }
539 
540 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
541 
543  reactionTable)
544 {
545 
546  //=========================================
547  // Diffusion controlled reaction model
548  //=========================================
554  G4VDNAReactionModel* reactionRadiusComputer =
556  reactionTable->PrintTable(reactionRadiusComputer);
557 
563  G4DNAMolecularStepByStepModel* stepByStep =
565  stepByStep->SetReactionModel(reactionRadiusComputer);
566 // ((G4DNAMoleculeEncounterStepper*) stepByStep->GetTimeStepper())->
567 // SetVerbose(5);
568 
569  RegisterTimeStepModel(stepByStep, 0);
570 }
virtual void ConstructDissociationChannels()
static G4Electron_aq * Definition()
static G4H2O * Definition()
Definition: G4H2O.cc:46
virtual void ConstructReactionTable(G4DNAMolecularReactionTable *reactionTable)
void AddDecayChannel(const G4MolecularConfiguration *molConf, const G4MolecularDissociationChannel *channel)
G4VEmModel * EmModel(G4int index=1) const
virtual void ConstructProcess()
virtual void ConstructTimeStepModel(G4DNAMolecularReactionTable *reactionTable)
static G4Electron * Definition()
Definition: G4Electron.cc:49
virtual ~G4EmDNAChemistry()
static constexpr double m3
Definition: G4SIunits.hh:131
float Avogadro
Definition: hepunit.py:253
void SetReaction(G4double observedReactionRate, G4MolecularConfiguration *reactive1, G4MolecularConfiguration *reactive2)
const XML_Char * s
Definition: expat.h:262
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
void SetChemistryList(G4VUserChemistryList *)
virtual void ConstructMolecule()
void AddProduct(const G4Molecule *, G4double=0)
static G4H3O * Definition()
Definition: G4H3O.cc:47
G4bool RegisterProcess(G4VProcess *process, G4ParticleDefinition *particle)
G4int AddElectron(G4int orbit, G4int number=1)
static G4MoleculeTable * Instance()
void PrintTable(G4VDNAReactionModel *=0)
static constexpr double eV
Definition: G4SIunits.hh:215
void AddProduct(G4MolecularConfiguration *molecule)
G4MolecularConfiguration * CreateConfiguration(const G4String &userIdentifier, const G4MoleculeDefinition *molDef, const G4String &configurationLabel, const G4ElectronOccupancy &eOcc)
static G4DNAChemistryManager * Instance()
G4MoleculeDefinitionIterator GetDefintionIterator()
G4ProcessManager * GetProcessManager() const
static G4PhysicsListHelper * GetPhysicsListHelper()
static G4H2O2 * Definition()
Definition: G4H2O2.cc:46
static G4OH * Definition()
Definition: G4OH.cc:46
G4int AddRestProcess(G4VProcess *aProcess, G4int ord=ordDefault)
void SetDecayDisplacer(const G4ParticleDefinition *, G4VMolecularDecayDisplacer *)
const G4ElectronOccupancy * GetGroundStateElectronOccupancy() const
G4MolecularConfiguration * NewConfigurationWithElectronOccupancy(const G4String &excitedStateLabel, const G4ElectronOccupancy &, double decayTime=0.)
static G4H2 * Definition()
Definition: G4H2.cc:46
void SetReactionModel(G4VDNAReactionModel *)
const XML_Char XML_Content * model
Definition: expat.h:151
static constexpr double mole
Definition: G4SIunits.hh:286
static constexpr double m2
Definition: G4SIunits.hh:130
static G4ProcessTable * GetProcessTable()
void SetVerboseLevel(G4int value)
Definition: G4VProcess.hh:437
G4int RemoveElectron(G4int orbit, G4int number=1)
#define G4_DECLARE_PHYSCONSTR_FACTORY(physics_constructor)
void RegisterTimeStepModel(G4VITStepModel *timeStepModel, double startingTime=0)
static G4Hydrogen * Definition()
Definition: G4Hydrogen.cc:46
G4MolecularConfiguration * GetConfiguration(const G4String &, bool mustExist=true)