Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4EmDNAChemistry_option1.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//
27
29#include "G4SystemOfUnits.hh"
30
34#include "G4ProcessManager.hh"
35
37
38// *** Processes and models for Geant4-DNA
39
41
42#include "G4DNAAttachment.hh"
43#include "G4DNAVibExcitation.hh"
45
46#include "G4DNAElastic.hh"
50
57
59
60// particles
61
62#include "G4Electron.hh"
63#include "G4Proton.hh"
64#include "G4GenericIon.hh"
65
66#include "G4MoleculeTable.hh"
67#include "G4H2O.hh"
68#include "G4H2.hh"
69#include "G4Hydrogen.hh"
70#include "G4OH.hh"
71#include "G4H3O.hh"
72#include "G4Electron_aq.hh"
73#include "G4H2O2.hh"
74
76#include "G4BuilderType.hh"
77
78/****/
80#include "G4ProcessVector.hh"
81#include "G4ProcessTable.hh"
84/****/
85
86// factory
88
90
91#include "G4Threading.hh"
92
95{
97}
98
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
102{
103}
104
105//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106
108{
109 //-----------------------------------
110 G4Electron::Definition(); // safety
111
112 //-----------------------------------
113 // Create the definition
120 G4H2::Definition();
121
122 //____________________________________________________________________________
123
126 9.46e-9 * (m2/s));
128 CreateConfiguration("OHm", // just a tag to store and retrieve from
129 // G4MoleculeTable
131 -1, // charge
132 5.3e-9 * (m2 / s));
133 OHm->SetMass(17.0079 * g / Avogadro * c_squared);
136 2.2e-9 * (m2/s));
141 G4MoleculeTable::Instance()->CreateConfiguration("H2", G4H2::Definition());
143 4.8e-9 * (m2/s));
146 2.3e-9 * (m2/s));
147}
148
149//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
150
152{
153 //-----------------------------------
154 //Get the molecular configuration
167
168 //-------------------------------------
169 //Define the decay channels
173
176
177 //////////////////////////////////////////////////////////
178 // EXCITATIONS //
179 //////////////////////////////////////////////////////////
180 G4DNAWaterExcitationStructure waterExcitation;
181 //--------------------------------------------------------
182 //---------------Excitation on the fifth layer------------
183
184 decCh1 = new G4MolecularDissociationChannel("A^1B_1_Relaxation");
185 decCh2 = new G4MolecularDissociationChannel("A^1B_1_DissociativeDecay");
186 //Decay 1 : OH + H
187 decCh1->SetEnergy(waterExcitation.ExcitationEnergy(0));
188 decCh1->SetProbability(0.35);
189 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::NoDisplacement);
190
191 decCh2->AddProduct(OH);
192 decCh2->AddProduct(H);
193 decCh2->SetProbability(0.65);
194 decCh2->SetDisplacementType(
195 G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
196
197// water->AddExcitedState("A^1B_1");
198 occ->RemoveElectron(4, 1); // this is the transition form ground state to
199 occ->AddElectron(5, 1); // the first unoccupied orbital: A^1B_1
200
201 water->NewConfigurationWithElectronOccupancy("A^1B_1", *occ);
202 water->AddDecayChannel("A^1B_1", decCh1);
203 water->AddDecayChannel("A^1B_1", decCh2);
204
205 //--------------------------------------------------------
206 //---------------Excitation on the fourth layer-----------
207 decCh1 = new G4MolecularDissociationChannel("B^1A_1_Relaxation_Channel");
208 decCh2 = new G4MolecularDissociationChannel("B^1A_1_DissociativeDecay");
210 "B^1A_1_AutoIonisation_Channel");
211
212 //Decay 1 : energy
213 decCh1->SetEnergy(waterExcitation.ExcitationEnergy(1));
214 decCh1->SetProbability(0.3);
215
216 //Decay 2 : 2OH + H_2
217 decCh2->AddProduct(H2);
218 decCh2->AddProduct(OH);
219 decCh2->AddProduct(OH);
220 decCh2->SetProbability(0.15);
221 decCh2->SetDisplacementType(
222 G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay);
223
224 //Decay 3 : OH + H_3Op + e_aq
225 decCh3->AddProduct(OH);
226 decCh3->AddProduct(H3O);
227 decCh3->AddProduct(e_aq);
228 decCh3->SetProbability(0.55);
229 decCh3->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
230
231 *occ = *(water->GetGroundStateElectronOccupancy());
232 occ->RemoveElectron(3); // this is the transition form ground state to
233 occ->AddElectron(5, 1); // the first unoccupied orbital: B^1A_1
234
235 water->NewConfigurationWithElectronOccupancy("B^1A_1", *occ);
236 water->AddDecayChannel("B^1A_1", decCh1);
237 water->AddDecayChannel("B^1A_1", decCh2);
238 water->AddDecayChannel("B^1A_1", decCh3);
239
240 //-------------------------------------------------------
241 //-------------------Excitation of 3rd layer-----------------
243 "Excitation3rdLayer_AutoIonisation_Channel");
245 "Excitation3rdLayer_Relaxation_Channel");
246
247 //Decay channel 1 : : OH + H_3Op + e_aq
248 decCh1->AddProduct(OH);
249 decCh1->AddProduct(H3O);
250 decCh1->AddProduct(e_aq);
251
252 decCh1->SetProbability(0.5);
253 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
254
255 //Decay channel 2 : energy
256 decCh2->SetEnergy(waterExcitation.ExcitationEnergy(2));
257 decCh2->SetProbability(0.5);
258
259 //Electronic configuration of this decay
260 *occ = *(water->GetGroundStateElectronOccupancy());
261 occ->RemoveElectron(2, 1);
262 occ->AddElectron(5, 1);
263
264 //Configure the water molecule
265 water->NewConfigurationWithElectronOccupancy("Excitation3rdLayer", *occ);
266 water->AddDecayChannel("Excitation3rdLayer", decCh1);
267 water->AddDecayChannel("Excitation3rdLayer", decCh2);
268
269 //-------------------------------------------------------
270 //-------------------Excitation of 2nd layer-----------------
272 "Excitation2ndLayer_AutoIonisation_Channel");
274 "Excitation2ndLayer_Relaxation_Channel");
275
276 //Decay Channel 1 : : OH + H_3Op + e_aq
277 decCh1->AddProduct(OH);
278 decCh1->AddProduct(H3O);
279 decCh1->AddProduct(e_aq);
280
281 decCh1->SetProbability(0.5);
282 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
283
284 //Decay channel 2 : energy
285 decCh2->SetEnergy(waterExcitation.ExcitationEnergy(3));
286 decCh2->SetProbability(0.5);
287
288 *occ = *(water->GetGroundStateElectronOccupancy());
289 occ->RemoveElectron(1, 1);
290 occ->AddElectron(5, 1);
291
292 water->NewConfigurationWithElectronOccupancy("Excitation2ndLayer", *occ);
293 water->AddDecayChannel("Excitation2ndLayer", decCh1);
294 water->AddDecayChannel("Excitation2ndLayer", decCh2);
295
296 //-------------------------------------------------------
297 //-------------------Excitation of 1st layer-----------------
299 "Excitation1stLayer_AutoIonisation_Channel");
301 "Excitation1stLayer_Relaxation_Channel");
302
303 *occ = *(water->GetGroundStateElectronOccupancy());
304 occ->RemoveElectron(0, 1);
305 occ->AddElectron(5, 1);
306
307 //Decay Channel 1 : : OH + H_3Op + e_aq
308 decCh1->AddProduct(OH);
309 decCh1->AddProduct(H3O);
310 decCh1->AddProduct(e_aq);
311 decCh1->SetProbability(0.5);
312 decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
313
314 //Decay channel 2 : energy
315 decCh2->SetEnergy(waterExcitation.ExcitationEnergy(4));
316 decCh2->SetProbability(0.5);
317
318 water->NewConfigurationWithElectronOccupancy("Excitation1stLayer", *occ);
319 water->AddDecayChannel("Excitation1stLayer", decCh1);
320 water->AddDecayChannel("Excitation1stLayer", decCh2);
321
322 /////////////////////////////////////////////////////////
323 // IONISATION //
324 /////////////////////////////////////////////////////////
325 //--------------------------------------------------------
326 //------------------- Ionisation -------------------------
327
328 decCh1 = new G4MolecularDissociationChannel("Ionisation_Channel");
329
330 //Decay Channel 1 : : OH + H_3Op
331 decCh1->AddProduct(H3O);
332 decCh1->AddProduct(OH);
333 decCh1->SetProbability(1);
334 decCh1->SetDisplacementType(
335 G4DNAWaterDissociationDisplacer::Ionisation_DissociationDecay);
336
337 *occ = *(water->GetGroundStateElectronOccupancy());
338 occ->RemoveElectron(4, 1);
339 // this is a ionized h2O with a hole in its last orbital
340 water->NewConfigurationWithElectronOccupancy("Ionisation5", *occ);
341 water->AddDecayChannel("Ionisation5",
342 decCh1);
343
344 *occ = *(water->GetGroundStateElectronOccupancy());
345 occ->RemoveElectron(3, 1);
346 water->NewConfigurationWithElectronOccupancy("Ionisation4", *occ);
347 water->AddDecayChannel("Ionisation4",
348 new G4MolecularDissociationChannel(*decCh1));
349
350 *occ = *(water->GetGroundStateElectronOccupancy());
351 occ->RemoveElectron(2, 1);
352 water->NewConfigurationWithElectronOccupancy("Ionisation3", *occ);
353 water->AddDecayChannel("Ionisation3",
354 new G4MolecularDissociationChannel(*decCh1));
355
356 *occ = *(water->GetGroundStateElectronOccupancy());
357 occ->RemoveElectron(1, 1);
358 water->NewConfigurationWithElectronOccupancy("Ionisation2", *occ);
359 water->AddDecayChannel("Ionisation2",
360 new G4MolecularDissociationChannel(*decCh1));
361
362 *occ = *(water->GetGroundStateElectronOccupancy());
363 occ->RemoveElectron(0, 1);
364 water->NewConfigurationWithElectronOccupancy("Ionisation1", *occ);
365 water->AddDecayChannel("Ionisation1",
366 new G4MolecularDissociationChannel(*decCh1));
367
368 //////////////////////////////////////////////////////////
369 // Dissociative Attachment //
370 //////////////////////////////////////////////////////////
371 decCh1 = new G4MolecularDissociationChannel("DissociativeAttachment");
372
373 //Decay 1 : 2OH + H_2
374 decCh1->AddProduct(H2);
375 decCh1->AddProduct(OHm);
376 decCh1->AddProduct(OH);
377 decCh1->SetProbability(1);
379 DissociativeAttachment);
380
381 *occ = *(water->GetGroundStateElectronOccupancy());
382 occ->AddElectron(5, 1); // H_2O^-
383 water->NewConfigurationWithElectronOccupancy("DissociativeAttachment", *occ);
384 water->AddDecayChannel("DissociativeAttachment", decCh1);
385
386 //////////////////////////////////////////////////////////
387 // Electron-hole recombination //
388 //////////////////////////////////////////////////////////
389 decCh1 = new G4MolecularDissociationChannel("H2Ovib_DissociationDecay1");
390 decCh2 = new G4MolecularDissociationChannel("H2Ovib_DissociationDecay2");
391 decCh3 = new G4MolecularDissociationChannel("H2Ovib_DissociationDecay3");
392
393 //Decay 1 : 2OH + H_2
394 decCh1->AddProduct(H2);
395 decCh1->AddProduct(OH);
396 decCh1->AddProduct(OH);
397 decCh1->SetProbability(0.15);
399 B1A1_DissociationDecay);
400
401 //Decay 2 : OH + H
402 decCh2->AddProduct(OH);
403 decCh2->AddProduct(H);
404 decCh2->SetProbability(0.55);
406 A1B1_DissociationDecay);
407
408 //Decay 3 : relaxation
409 decCh3->SetProbability(0.30);
410
411 const auto pH2Ovib = G4H2O::Definition()->NewConfiguration("H2Ovib");
412 assert(pH2Ovib != nullptr);
413
414 water->AddDecayChannel(pH2Ovib, decCh1);
415 water->AddDecayChannel(pH2Ovib, decCh2);
416 water->AddDecayChannel(pH2Ovib, decCh3);
417
418 delete occ;
419}
420
421//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
422
424 theReactionTable)
425{
426 //-----------------------------------
427 //Get the molecular configuration
442
443 //------------------------------------------------------------------
444 // e_aq + e_aq + 2H2O -> H2 + 2OH-
445 G4DNAMolecularReactionData* reactionData =
446 new G4DNAMolecularReactionData(0.636e10 * (1e-3 * m3 / (mole * s)), e_aq, e_aq);
447 reactionData->AddProduct(OHm);
448 reactionData->AddProduct(OHm);
449 reactionData->AddProduct(H2);
450 theReactionTable->SetReaction(reactionData);
451 //------------------------------------------------------------------
452 // e_aq + *OH -> OH-
453 reactionData = new G4DNAMolecularReactionData(
454 2.95e10 * (1e-3 * m3 / (mole * s)), e_aq, OH);
455 reactionData->AddProduct(OHm);
456 theReactionTable->SetReaction(reactionData);
457 //------------------------------------------------------------------
458 // e_aq + H* + H2O -> H2 + OH-
459 reactionData = new G4DNAMolecularReactionData(
460 2.50e10 * (1e-3 * m3 / (mole * s)), e_aq, H);
461 reactionData->AddProduct(OHm);
462 reactionData->AddProduct(H2);
463 theReactionTable->SetReaction(reactionData);
464 //------------------------------------------------------------------
465 // e_aq + H3O+ -> H* + H2O
466 reactionData = new G4DNAMolecularReactionData(
467 2.11e10 * (1e-3 * m3 / (mole * s)), e_aq, H3Op);
468 reactionData->AddProduct(H);
469 theReactionTable->SetReaction(reactionData);
470 //------------------------------------------------------------------
471 // e_aq + H2O2 -> OH- + *OH
472 reactionData = new G4DNAMolecularReactionData(
473 1.10e10 * (1e-3 * m3 / (mole * s)), e_aq, H2O2);
474 reactionData->AddProduct(OHm);
475 reactionData->AddProduct(OH);
476 theReactionTable->SetReaction(reactionData);
477 //------------------------------------------------------------------
478 // *OH + *OH -> H2O2
479 reactionData = new G4DNAMolecularReactionData(
480 0.55e10 * (1e-3 * m3 / (mole * s)), OH, OH);
481 reactionData->AddProduct(H2O2);
482 theReactionTable->SetReaction(reactionData);
483 //------------------------------------------------------------------
484 // *OH + *H -> H2O
485 theReactionTable->SetReaction(1.55e10 * (1e-3 * m3 / (mole * s)), OH, H);
486 //------------------------------------------------------------------
487 // *H + *H -> H2
488 reactionData = new G4DNAMolecularReactionData(
489 0.503e10 * (1e-3 * m3 / (mole * s)), H, H);
490 reactionData->AddProduct(H2);
491 theReactionTable->SetReaction(reactionData);
492 //------------------------------------------------------------------
493 // H3O+ + OH- -> 2H2O
494 theReactionTable->SetReaction(1.13e11 * (1e-3 * m3 / (mole * s)), H3Op, OHm);
495 //------------------------------------------------------------------
496}
497
498//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
499
501{
503
504 //===============================================================
505 // Extend vibrational to low energy
506 // Anyway, solvation of electrons is taken into account from 7.4 eV
507 // So below this threshold, for now, no accurate modeling is done
508 //
509 G4VProcess* process =
511 FindProcess("e-_G4DNAVibExcitation", "e-");
512
513 if (process)
514 {
515 G4DNAVibExcitation* vibExcitation = (G4DNAVibExcitation*) process;
516 G4VEmModel* model = vibExcitation->EmModel();
517 G4DNASancheExcitationModel* sancheExcitationMod =
518 dynamic_cast<G4DNASancheExcitationModel*>(model);
519 if(sancheExcitationMod)
520 {
521 sancheExcitationMod->ExtendLowEnergyLimit(0.025 * eV);
522 }
523 }
524
525 //===============================================================
526 // *** Electron Solvatation ***
527 //
528 process =
530 FindProcess("e-_G4DNAElectronSolvation", "e-");
531
532 if (process == 0)
533 {
534 ph->RegisterProcess(
535 new G4DNAElectronSolvation("e-_G4DNAElectronSolvation"),
537 }
538
539 //===============================================================
540 // Define processes for molecules
541 //
542 G4MoleculeTable* theMoleculeTable = G4MoleculeTable::Instance();
544 theMoleculeTable->GetDefintionIterator();
545 iterator.reset();
546 while (iterator())
547 {
548 G4MoleculeDefinition* moleculeDef = iterator.value();
549
550 if (moleculeDef != G4H2O::Definition())
551 {
552 // G4cout << "Brownian motion added for : "
553 // << moleculeDef->GetName() << G4endl;
555 // brown->SetVerboseLevel(4);
556 ph->RegisterProcess(brown, moleculeDef);
557 }
558 else
559 {
560 moleculeDef->GetProcessManager()
562 G4DNAMolecularDissociation* dissociationProcess =
563 new G4DNAMolecularDissociation("H2O_DNAMolecularDecay");
564 dissociationProcess->SetDisplacer(
565 moleculeDef, new G4DNAWaterDissociationDisplacer);
566 dissociationProcess->SetVerboseLevel(1);
567// ph->RegisterProcess(dissociationProcess, moleculeDef);
568
569 moleculeDef->GetProcessManager()
570 ->AddRestProcess(dissociationProcess, 1);
571 }
572 /*
573 * Warning : end of particles and processes are needed by
574 * EM Physics builders
575 */
576 }
577
579}
580
581//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
582
584 reactionTable)
585{
586
587 //=========================================
588 // Diffusion controlled reaction model
589 //=========================================
590 /**
591 * The reaction model defines how to compute the reaction range between
592 * molecules
593 */
594
595 G4VDNAReactionModel* reactionRadiusComputer =
597 reactionTable->PrintTable(reactionRadiusComputer);
598
599 /**
600 * The StepByStep model tells the step manager how to behave before and
601 * after each step, how to compute the time steps.
602 */
603
606 stepByStep->SetReactionModel(reactionRadiusComputer);
607// ((G4DNAMoleculeEncounterStepper*) stepByStep->GetTimeStepper())->
608// SetVerbose(5);
609
610 RegisterTimeStepModel(stepByStep, 0);
611}
#define G4_DECLARE_PHYSCONSTR_FACTORY(physics_constructor)
static G4DNAChemistryManager * Instance()
void SetChemistryList(G4VUserChemistryList &)
void SetDisplacer(Species *, Displacer *)
void PrintTable(G4VDNAReactionModel *=0)
void SetReaction(G4double observedReactionRate, Reactant *reactive1, Reactant *reactive2)
void SetReactionModel(G4VDNAReactionModel *)
G4int AddElectron(G4int orbit, G4int number=1)
G4int RemoveElectron(G4int orbit, G4int number=1)
static G4Electron_aq * Definition()
static G4Electron * Definition()
Definition: G4Electron.cc:48
virtual void ConstructReactionTable(G4DNAMolecularReactionTable *reactionTable)
virtual void ConstructTimeStepModel(G4DNAMolecularReactionTable *reactionTable)
static G4H2O2 * Definition()
Definition: G4H2O2.cc:45
static G4H2O * Definition()
Definition: G4H2O.cc:42
static G4H3O * Definition()
Definition: G4H3O.cc:46
static G4Hydrogen * Definition()
Definition: G4Hydrogen.cc:45
void AddProduct(Product *, G4double displacement=0.)
const G4ElectronOccupancy * GetGroundStateElectronOccupancy() const
void AddDecayChannel(const G4MolecularConfiguration *molConf, const G4MolecularDissociationChannel *channel)
G4MolecularConfiguration * NewConfiguration(const G4String &excitedStateLabel)
G4MolecularConfiguration * NewConfigurationWithElectronOccupancy(const G4String &excitedStateLabel, const G4ElectronOccupancy &, double decayTime=0.)
G4MolecularConfiguration * GetConfiguration(const G4String &, bool mustExist=true)
G4MolecularConfiguration * CreateConfiguration(const G4String &userIdentifier, const G4MoleculeDefinition *molDef, const G4String &configurationLabel, const G4ElectronOccupancy &eOcc)
G4MoleculeDefinitionIterator GetDefintionIterator()
static G4MoleculeTable * Instance()
static G4OH * Definition()
Definition: G4OH.cc:45
G4ProcessManager * GetProcessManager() const
G4bool RegisterProcess(G4VProcess *process, G4ParticleDefinition *particle)
static G4PhysicsListHelper * GetPhysicsListHelper()
G4int AddRestProcess(G4VProcess *aProcess, G4int ord=ordDefault)
static G4ProcessTable * GetProcessTable()
G4VEmModel * EmModel(size_t index=0) const
void SetVerboseLevel(G4int value)
Definition: G4VProcess.hh:416
void RegisterTimeStepModel(G4VITStepModel *timeStepModel, G4double startingTime=0)