ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EmDNAChemistry.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EmDNAChemistry.cc
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 
31 #include "G4DNAWaterDissociationDisplacer.hh"
32 #include "G4DNAChemistryManager.hh"
33 #include "G4DNAWaterExcitationStructure.hh"
34 #include "G4ProcessManager.hh"
35 
36 #include "G4DNAGenericIonsManager.hh"
37 
38 // *** Processes and models for Geant4-DNA
39 
40 #include "G4DNAElectronSolvation.hh"
41 
42 #include "G4DNAAttachment.hh"
43 #include "G4DNAVibExcitation.hh"
44 
45 #include "G4DNAElastic.hh"
46 #include "G4DNAChampionElasticModel.hh"
47 #include "G4DNAScreenedRutherfordElasticModel.hh"
48 #include "G4DNAUeharaScreenedRutherfordElasticModel.hh"
49 
50 #include "G4DNAMolecularDissociation.hh"
51 #include "G4DNABrownianTransportation.hh"
52 #include "G4DNAMolecularReactionTable.hh"
53 #include "G4DNAMolecularStepByStepModel.hh"
54 #include "G4VDNAReactionModel.hh"
55 #include "G4DNASmoluchowskiReactionModel.hh"
56 
57 #include "G4DNAElectronHoleRecombination.hh"
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 /****/
78 #include "G4DNAMoleculeEncounterStepper.hh"
79 #include "G4ProcessVector.hh"
80 #include "G4ProcessTable.hh"
81 #include "G4DNASecondOrderReaction.hh"
82 #include "G4MolecularConfiguration.hh"
83 /****/
84 
85 // factory
86 #include "G4PhysicsConstructorFactory.hh"
87 
88 G4_DECLARE_PHYSCONSTR_FACTORY(G4EmDNAChemistry);
89 
90 #include "G4Threading.hh"
91 
92 G4EmDNAChemistry::G4EmDNAChemistry() :
93  G4VUserChemistryList(true)
94 {
95  G4DNAChemistryManager::Instance()->SetChemistryList(this);
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 
100 G4EmDNAChemistry::~G4EmDNAChemistry()
101 {
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
106 void G4EmDNAChemistry::ConstructMolecule()
107 {
108  //-----------------------------------
109  G4Electron::Definition(); // safety
110 
111  //-----------------------------------
112  // Create the definition
113  G4H2O::Definition();
114  G4Hydrogen::Definition();
115  G4H3O::Definition();
116  G4OH::Definition();
117  G4Electron_aq::Definition();
118  G4H2O2::Definition();
119  G4H2::Definition();
120 
121  //____________________________________________________________________________
122 
123  G4MoleculeTable::Instance()->CreateConfiguration("H3Op", G4H3O::Definition());
124  G4MolecularConfiguration* OHm = G4MoleculeTable::Instance()->
125  CreateConfiguration("OHm", // just a tag to store and retrieve from
126  // G4MoleculeTable
127  G4OH::Definition(),
128  -1, // charge
129  5.0e-9 * (m2 / s));
130  OHm->SetMass(17.0079 * g / Avogadro * c_squared);
131  G4MoleculeTable::Instance()->CreateConfiguration("OH", G4OH::Definition());
132  G4MoleculeTable::Instance()->CreateConfiguration("e_aq",
133  G4Electron_aq::Definition());
134  G4MoleculeTable::Instance()->CreateConfiguration("H",
135  G4Hydrogen::Definition());
136  G4MoleculeTable::Instance()->CreateConfiguration("H2", G4H2::Definition());
137  G4MoleculeTable::Instance()->CreateConfiguration("H2O2", G4H2O2::Definition());
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 
142 void G4EmDNAChemistry::ConstructDissociationChannels()
143 {
144  //-----------------------------------
145  //Get the molecular configuration
146  G4MolecularConfiguration* OH =
147  G4MoleculeTable::Instance()->GetConfiguration("OH");
148  G4MolecularConfiguration* OHm =
149  G4MoleculeTable::Instance()->GetConfiguration("OHm");
150  G4MolecularConfiguration* e_aq =
151  G4MoleculeTable::Instance()->GetConfiguration("e_aq");
152  G4MolecularConfiguration* H2 =
153  G4MoleculeTable::Instance()->GetConfiguration("H2");
154  G4MolecularConfiguration* H3O =
155  G4MoleculeTable::Instance()->GetConfiguration("H3Op");
156  G4MolecularConfiguration* H =
157  G4MoleculeTable::Instance()->GetConfiguration("H");
158 
159  //-------------------------------------
160  //Define the decay channels
161  G4MoleculeDefinition* water = G4H2O::Definition();
162  G4MolecularDissociationChannel* decCh1;
163  G4MolecularDissociationChannel* decCh2;
164 
165  G4ElectronOccupancy* occ = new G4ElectronOccupancy(
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");
200  G4MolecularDissociationChannel* decCh3 = new G4MolecularDissociationChannel(
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);
369  decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::
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 
382 void G4EmDNAChemistry::ConstructReactionTable(G4DNAMolecularReactionTable*
383  theReactionTable)
384 {
385  //-----------------------------------
386  //Get the molecular configuration
387  G4MolecularConfiguration* OH =
388  G4MoleculeTable::Instance()->GetConfiguration("OH");
389  G4MolecularConfiguration* OHm =
390  G4MoleculeTable::Instance()->GetConfiguration("OHm");
391  G4MolecularConfiguration* e_aq =
392  G4MoleculeTable::Instance()->GetConfiguration("e_aq");
393  G4MolecularConfiguration* H2 =
394  G4MoleculeTable::Instance()->GetConfiguration("H2");
395  G4MolecularConfiguration* H3Op =
396  G4MoleculeTable::Instance()->GetConfiguration("H3Op");
397  G4MolecularConfiguration* H =
398  G4MoleculeTable::Instance()->GetConfiguration("H");
399  G4MolecularConfiguration* H2O2 =
400  G4MoleculeTable::Instance()->GetConfiguration("H2O2");
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 
459 void G4EmDNAChemistry::ConstructProcess()
460 {
461  auto pPhysicsListHelper = G4PhysicsListHelper::GetPhysicsListHelper();
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* pProcess = G4ProcessTable::GetProcessTable()->
469  FindProcess("e-_G4DNAVibExcitation", "e-");
470 
471  if (pProcess != nullptr)
472  {
473  G4DNAVibExcitation* pVibExcitation = (G4DNAVibExcitation*) pProcess;
474  G4VEmModel* pModel = pVibExcitation->EmModel();
475  G4DNASancheExcitationModel* pSancheExcitationMod =
476  dynamic_cast<G4DNASancheExcitationModel*>(pModel);
477  if(pSancheExcitationMod != nullptr)
478  {
479  pSancheExcitationMod->ExtendLowEnergyLimit(0.025 * eV);
480  }
481  }
482 
483  //===============================================================
484  // Electron Solvatation
485  //
486  pProcess = G4ProcessTable::GetProcessTable()->FindProcess("e-_G4DNAElectronSolvation", "e-");
487 
488  if (pProcess == nullptr)
489  {
490  pPhysicsListHelper->RegisterProcess(new G4DNAElectronSolvation("e-_G4DNAElectronSolvation"),
491  G4Electron::Definition());
492  }
493 
494  //===============================================================
495  // Define processes for molecules
496  //
497  G4MoleculeTable* pMoleculeTable = G4MoleculeTable::Instance();
498  G4MoleculeDefinitionIterator iterator = pMoleculeTable->GetDefintionIterator();
499  iterator.reset();
500  while (iterator())
501  {
502  G4MoleculeDefinition* pMoleculeDef = iterator.value();
503 
504  if (pMoleculeDef != G4H2O::Definition())
505  {
506  G4DNABrownianTransportation* pBrownianTransport = new G4DNABrownianTransportation();
507  pPhysicsListHelper->RegisterProcess(pBrownianTransport, pMoleculeDef);
508  }
509  else
510  {
511  pMoleculeDef->GetProcessManager()->AddRestProcess(new G4DNAElectronHoleRecombination(), 2);
512  G4DNAMolecularDissociation* pDissociationProcess = new G4DNAMolecularDissociation("H2O_DNAMolecularDecay");
513  pDissociationProcess->SetDisplacer(pMoleculeDef, new G4DNAWaterDissociationDisplacer);
514  pDissociationProcess->SetVerboseLevel(1);
515 
516  pMoleculeDef->GetProcessManager()->AddRestProcess(pDissociationProcess, 1);
517  }
518  }
519 
520  G4DNAChemistryManager::Instance()->Initialize();
521 }
522 
523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
524 
525 void G4EmDNAChemistry::ConstructTimeStepModel(G4DNAMolecularReactionTable*
526  reactionTable)
527 {
528  G4VDNAReactionModel* reactionRadiusComputer = new G4DNASmoluchowskiReactionModel();
529  reactionTable->PrintTable(reactionRadiusComputer);
530 
531  G4DNAMolecularStepByStepModel* stepByStep = new G4DNAMolecularStepByStepModel();
532  stepByStep->SetReactionModel(reactionRadiusComputer);
533 // ((G4DNAMoleculeEncounterStepper*) stepByStep->GetTimeStepper())->
534 // SetVerbose(5);
535 
536  RegisterTimeStepModel(stepByStep, 0);
537 }