ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EmDNAChemistry_option1.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EmDNAChemistry_option1.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 //
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  9.46e-9 * (m2/s));
127  CreateConfiguration("OHm", // just a tag to store and retrieve from
128  // G4MoleculeTable
130  -1, // charge
131  5.3e-9 * (m2 / s));
132  OHm->SetMass(17.0079 * g / Avogadro * c_squared);
135  2.2e-9 * (m2/s));
142  4.8e-9 * (m2/s));
145  2.3e-9 * (m2/s));
146 }
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149 
151 {
152  //-----------------------------------
153  //Get the molecular configuration
166 
167  //-------------------------------------
168  //Define the decay channels
172 
174  *(water->GetGroundStateElectronOccupancy()));
175 
177  // EXCITATIONS //
179  G4DNAWaterExcitationStructure waterExcitation;
180  //--------------------------------------------------------
181  //---------------Excitation on the fifth layer------------
182 
183  decCh1 = new G4MolecularDissociationChannel("A^1B_1_Relaxation");
184  decCh2 = new G4MolecularDissociationChannel("A^1B_1_DissociativeDecay");
185  //Decay 1 : OH + H
186  decCh1->SetEnergy(waterExcitation.ExcitationEnergy(0));
187  decCh1->SetProbability(0.35);
188  decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::NoDisplacement);
189 
190  decCh2->AddProduct(OH);
191  decCh2->AddProduct(H);
192  decCh2->SetProbability(0.65);
193  decCh2->SetDisplacementType(
194  G4DNAWaterDissociationDisplacer::A1B1_DissociationDecay);
195 
196 // water->AddExcitedState("A^1B_1");
197  occ->RemoveElectron(4, 1); // this is the transition form ground state to
198  occ->AddElectron(5, 1); // the first unoccupied orbital: A^1B_1
199 
200  water->NewConfigurationWithElectronOccupancy("A^1B_1", *occ);
201  water->AddDecayChannel("A^1B_1", decCh1);
202  water->AddDecayChannel("A^1B_1", decCh2);
203 
204  //--------------------------------------------------------
205  //---------------Excitation on the fourth layer-----------
206  decCh1 = new G4MolecularDissociationChannel("B^1A_1_Relaxation_Channel");
207  decCh2 = new G4MolecularDissociationChannel("B^1A_1_DissociativeDecay");
209  "B^1A_1_AutoIonisation_Channel");
210 
211  //Decay 1 : energy
212  decCh1->SetEnergy(waterExcitation.ExcitationEnergy(1));
213  decCh1->SetProbability(0.3);
214 
215  //Decay 2 : 2OH + H_2
216  decCh2->AddProduct(H2);
217  decCh2->AddProduct(OH);
218  decCh2->AddProduct(OH);
219  decCh2->SetProbability(0.15);
220  decCh2->SetDisplacementType(
221  G4DNAWaterDissociationDisplacer::B1A1_DissociationDecay);
222 
223  //Decay 3 : OH + H_3Op + e_aq
224  decCh3->AddProduct(OH);
225  decCh3->AddProduct(H3O);
226  decCh3->AddProduct(e_aq);
227  decCh3->SetProbability(0.55);
228  decCh3->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
229 
230  *occ = *(water->GetGroundStateElectronOccupancy());
231  occ->RemoveElectron(3); // this is the transition form ground state to
232  occ->AddElectron(5, 1); // the first unoccupied orbital: B^1A_1
233 
234  water->NewConfigurationWithElectronOccupancy("B^1A_1", *occ);
235  water->AddDecayChannel("B^1A_1", decCh1);
236  water->AddDecayChannel("B^1A_1", decCh2);
237  water->AddDecayChannel("B^1A_1", decCh3);
238 
239  //-------------------------------------------------------
240  //-------------------Excitation of 3rd layer-----------------
241  decCh1 = new G4MolecularDissociationChannel(
242  "Excitation3rdLayer_AutoIonisation_Channel");
243  decCh2 = new G4MolecularDissociationChannel(
244  "Excitation3rdLayer_Relaxation_Channel");
245 
246  //Decay channel 1 : : OH + H_3Op + e_aq
247  decCh1->AddProduct(OH);
248  decCh1->AddProduct(H3O);
249  decCh1->AddProduct(e_aq);
250 
251  decCh1->SetProbability(0.5);
252  decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
253 
254  //Decay channel 2 : energy
255  decCh2->SetEnergy(waterExcitation.ExcitationEnergy(2));
256  decCh2->SetProbability(0.5);
257 
258  //Electronic configuration of this decay
259  *occ = *(water->GetGroundStateElectronOccupancy());
260  occ->RemoveElectron(2, 1);
261  occ->AddElectron(5, 1);
262 
263  //Configure the water molecule
264  water->NewConfigurationWithElectronOccupancy("Excitation3rdLayer", *occ);
265  water->AddDecayChannel("Excitation3rdLayer", decCh1);
266  water->AddDecayChannel("Excitation3rdLayer", decCh2);
267 
268  //-------------------------------------------------------
269  //-------------------Excitation of 2nd layer-----------------
270  decCh1 = new G4MolecularDissociationChannel(
271  "Excitation2ndLayer_AutoIonisation_Channel");
272  decCh2 = new G4MolecularDissociationChannel(
273  "Excitation2ndLayer_Relaxation_Channel");
274 
275  //Decay Channel 1 : : OH + H_3Op + e_aq
276  decCh1->AddProduct(OH);
277  decCh1->AddProduct(H3O);
278  decCh1->AddProduct(e_aq);
279 
280  decCh1->SetProbability(0.5);
281  decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
282 
283  //Decay channel 2 : energy
284  decCh2->SetEnergy(waterExcitation.ExcitationEnergy(3));
285  decCh2->SetProbability(0.5);
286 
287  *occ = *(water->GetGroundStateElectronOccupancy());
288  occ->RemoveElectron(1, 1);
289  occ->AddElectron(5, 1);
290 
291  water->NewConfigurationWithElectronOccupancy("Excitation2ndLayer", *occ);
292  water->AddDecayChannel("Excitation2ndLayer", decCh1);
293  water->AddDecayChannel("Excitation2ndLayer", decCh2);
294 
295  //-------------------------------------------------------
296  //-------------------Excitation of 1st layer-----------------
297  decCh1 = new G4MolecularDissociationChannel(
298  "Excitation1stLayer_AutoIonisation_Channel");
299  decCh2 = new G4MolecularDissociationChannel(
300  "Excitation1stLayer_Relaxation_Channel");
301 
302  *occ = *(water->GetGroundStateElectronOccupancy());
303  occ->RemoveElectron(0, 1);
304  occ->AddElectron(5, 1);
305 
306  //Decay Channel 1 : : OH + H_3Op + e_aq
307  decCh1->AddProduct(OH);
308  decCh1->AddProduct(H3O);
309  decCh1->AddProduct(e_aq);
310  decCh1->SetProbability(0.5);
311  decCh1->SetDisplacementType(G4DNAWaterDissociationDisplacer::AutoIonisation);
312 
313  //Decay channel 2 : energy
314  decCh2->SetEnergy(waterExcitation.ExcitationEnergy(4));
315  decCh2->SetProbability(0.5);
316 
317  water->NewConfigurationWithElectronOccupancy("Excitation1stLayer", *occ);
318  water->AddDecayChannel("Excitation1stLayer", decCh1);
319  water->AddDecayChannel("Excitation1stLayer", decCh2);
320 
322  // IONISATION //
324  //--------------------------------------------------------
325  //------------------- Ionisation -------------------------
326 
327  decCh1 = new G4MolecularDissociationChannel("Ionisation_Channel");
328 
329  //Decay Channel 1 : : OH + H_3Op
330  decCh1->AddProduct(H3O);
331  decCh1->AddProduct(OH);
332  decCh1->SetProbability(1);
333  decCh1->SetDisplacementType(
334  G4DNAWaterDissociationDisplacer::Ionisation_DissociationDecay);
335 
336  *occ = *(water->GetGroundStateElectronOccupancy());
337  occ->RemoveElectron(4, 1);
338  // this is a ionized h2O with a hole in its last orbital
339  water->NewConfigurationWithElectronOccupancy("Ionisation5", *occ);
340  water->AddDecayChannel("Ionisation5",
341  decCh1);
342 
343  *occ = *(water->GetGroundStateElectronOccupancy());
344  occ->RemoveElectron(3, 1);
345  water->NewConfigurationWithElectronOccupancy("Ionisation4", *occ);
346  water->AddDecayChannel("Ionisation4",
347  new G4MolecularDissociationChannel(*decCh1));
348 
349  *occ = *(water->GetGroundStateElectronOccupancy());
350  occ->RemoveElectron(2, 1);
351  water->NewConfigurationWithElectronOccupancy("Ionisation3", *occ);
352  water->AddDecayChannel("Ionisation3",
353  new G4MolecularDissociationChannel(*decCh1));
354 
355  *occ = *(water->GetGroundStateElectronOccupancy());
356  occ->RemoveElectron(1, 1);
357  water->NewConfigurationWithElectronOccupancy("Ionisation2", *occ);
358  water->AddDecayChannel("Ionisation2",
359  new G4MolecularDissociationChannel(*decCh1));
360 
361  *occ = *(water->GetGroundStateElectronOccupancy());
362  occ->RemoveElectron(0, 1);
363  water->NewConfigurationWithElectronOccupancy("Ionisation1", *occ);
364  water->AddDecayChannel("Ionisation1",
365  new G4MolecularDissociationChannel(*decCh1));
366 
368  // Dissociative Attachment //
370  decCh1 = new G4MolecularDissociationChannel("DissociativeAttachment");
371 
372  //Decay 1 : 2OH + H_2
373  decCh1->AddProduct(H2);
374  decCh1->AddProduct(OHm);
375  decCh1->AddProduct(OH);
376  decCh1->SetProbability(1);
378  DissociativeAttachment);
379 
380  *occ = *(water->GetGroundStateElectronOccupancy());
381  occ->AddElectron(5, 1); // H_2O^-
382  water->NewConfigurationWithElectronOccupancy("DissociativeAttachment", *occ);
383  water->AddDecayChannel("DissociativeAttachment", decCh1);
384 
385  delete occ;
386 }
387 
388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
389 
391  theReactionTable)
392 {
393  //-----------------------------------
394  //Get the molecular configuration
409 
410  //------------------------------------------------------------------
411  // e_aq + e_aq + 2H2O -> H2 + 2OH-
412  G4DNAMolecularReactionData* reactionData =
413  new G4DNAMolecularReactionData(0.636e10 * (1e-3 * m3 / (mole * s)), e_aq, e_aq);
414  reactionData->AddProduct(OHm);
415  reactionData->AddProduct(OHm);
416  reactionData->AddProduct(H2);
417  theReactionTable->SetReaction(reactionData);
418  //------------------------------------------------------------------
419  // e_aq + *OH -> OH-
420  reactionData = new G4DNAMolecularReactionData(
421  2.95e10 * (1e-3 * m3 / (mole * s)), e_aq, OH);
422  reactionData->AddProduct(OHm);
423  theReactionTable->SetReaction(reactionData);
424  //------------------------------------------------------------------
425  // e_aq + H* + H2O -> H2 + OH-
426  reactionData = new G4DNAMolecularReactionData(
427  2.50e10 * (1e-3 * m3 / (mole * s)), e_aq, H);
428  reactionData->AddProduct(OHm);
429  reactionData->AddProduct(H2);
430  theReactionTable->SetReaction(reactionData);
431  //------------------------------------------------------------------
432  // e_aq + H3O+ -> H* + H2O
433  reactionData = new G4DNAMolecularReactionData(
434  2.11e10 * (1e-3 * m3 / (mole * s)), e_aq, H3Op);
435  reactionData->AddProduct(H);
436  theReactionTable->SetReaction(reactionData);
437  //------------------------------------------------------------------
438  // e_aq + H2O2 -> OH- + *OH
439  reactionData = new G4DNAMolecularReactionData(
440  1.10e10 * (1e-3 * m3 / (mole * s)), e_aq, H2O2);
441  reactionData->AddProduct(OHm);
442  reactionData->AddProduct(OH);
443  theReactionTable->SetReaction(reactionData);
444  //------------------------------------------------------------------
445  // *OH + *OH -> H2O2
446  reactionData = new G4DNAMolecularReactionData(
447  0.55e10 * (1e-3 * m3 / (mole * s)), OH, OH);
448  reactionData->AddProduct(H2O2);
449  theReactionTable->SetReaction(reactionData);
450  //------------------------------------------------------------------
451  // *OH + *H -> H2O
452  theReactionTable->SetReaction(1.55e10 * (1e-3 * m3 / (mole * s)), OH, H);
453  //------------------------------------------------------------------
454  // *H + *H -> H2
455  reactionData = new G4DNAMolecularReactionData(
456  0.503e10 * (1e-3 * m3 / (mole * s)), H, H);
457  reactionData->AddProduct(H2);
458  theReactionTable->SetReaction(reactionData);
459  //------------------------------------------------------------------
460  // H3O+ + OH- -> 2H2O
461  theReactionTable->SetReaction(1.13e11 * (1e-3 * m3 / (mole * s)), H3Op, OHm);
462  //------------------------------------------------------------------
463 }
464 
465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
466 
468 {
470 
471  //===============================================================
472  // Extend vibrational to low energy
473  // Anyway, solvation of electrons is taken into account from 7.4 eV
474  // So below this threshold, for now, no accurate modeling is done
475  //
476  G4VProcess* process =
478  FindProcess("e-_G4DNAVibExcitation", "e-");
479 
480  if (process)
481  {
482  G4DNAVibExcitation* vibExcitation = (G4DNAVibExcitation*) process;
483  G4VEmModel* model = vibExcitation->EmModel();
484  G4DNASancheExcitationModel* sancheExcitationMod =
485  dynamic_cast<G4DNASancheExcitationModel*>(model);
486  if(sancheExcitationMod)
487  {
488  sancheExcitationMod->ExtendLowEnergyLimit(0.025 * eV);
489  }
490  }
491 
492  //===============================================================
493  // *** Electron Solvatation ***
494  //
495  process =
497  FindProcess("e-_G4DNAElectronSolvation", "e-");
498 
499  if (process == 0)
500  {
501  ph->RegisterProcess(
502  new G4DNAElectronSolvation("e-_G4DNAElectronSolvation"),
504  }
505 
506  //===============================================================
507  // Define processes for molecules
508  //
509  G4MoleculeTable* theMoleculeTable = G4MoleculeTable::Instance();
511  theMoleculeTable->GetDefintionIterator();
512  iterator.reset();
513  while (iterator())
514  {
515  G4MoleculeDefinition* moleculeDef = iterator.value();
516 
517  if (moleculeDef != G4H2O::Definition())
518  {
519  // G4cout << "Brownian motion added for : "
520  // << moleculeDef->GetName() << G4endl;
522  // brown->SetVerboseLevel(4);
523  ph->RegisterProcess(brown, moleculeDef);
524  }
525  else
526  {
527  moleculeDef->GetProcessManager()
529  G4DNAMolecularDissociation* dissociationProcess =
530  new G4DNAMolecularDissociation("H2O_DNAMolecularDecay");
531  dissociationProcess->SetDisplacer(
532  moleculeDef, new G4DNAWaterDissociationDisplacer);
533  dissociationProcess->SetVerboseLevel(1);
534 // ph->RegisterProcess(dissociationProcess, moleculeDef);
535 
536  moleculeDef->GetProcessManager()
537  ->AddRestProcess(dissociationProcess, 1);
538  }
539  /*
540  * Warning : end of particles and processes are needed by
541  * EM Physics builders
542  */
543  }
544 
546 }
547 
548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
549 
551  reactionTable)
552 {
553 
554  //=========================================
555  // Diffusion controlled reaction model
556  //=========================================
562  G4VDNAReactionModel* reactionRadiusComputer =
564  reactionTable->PrintTable(reactionRadiusComputer);
565 
571  G4DNAMolecularStepByStepModel* stepByStep =
573  stepByStep->SetReactionModel(reactionRadiusComputer);
574 // ((G4DNAMoleculeEncounterStepper*) stepByStep->GetTimeStepper())->
575 // SetVerbose(5);
576 
577  RegisterTimeStepModel(stepByStep, 0);
578 }