| Geant4 Cross Reference |
1 // ******************************************* 1
2 // * License and Disclaimer
3 // *
4 // * The Geant4 software is copyright of th
5 // * the Geant4 Collaboration. It is provided
6 // * conditions of the Geant4 Software License
7 // * LICENSE and available at http://cern.ch/
8 // * include a list of copyright holders.
9 // *
10 // * Neither the authors of this software syst
11 // * institutes,nor the agencies providing fin
12 // * work make any representation or warran
13 // * regarding this software system or assum
14 // * use. Please see the license in the file
15 // * for the full disclaimer and the limitatio
16 // *
17 // * This code implementation is the result
18 // * technical work of the GEANT4 collaboratio
19 // * By using, copying, modifying or distri
20 // * any work based on the software) you ag
21 // * use in resulting scientific publicati
22 // * acceptance of all terms of the Geant4 Sof
23 // *******************************************
24 //
25 /// \file scavenger/src/EmDNAChemistry.cc
26 /// \brief Implementation of the scavenger::Em
27
28 #include "EmDNAChemistry.hh"
29
30 #include "G4DNAChemistryManager.hh"
31 #include "G4DNAElectronHoleRecombination.hh"
32 #include "G4DNAElectronSolvation.hh"
33 #include "G4DNAIRT.hh"
34 #include "G4DNAMolecularDissociation.hh"
35 #include "G4DNAMolecularIRTModel.hh"
36 #include "G4DNAMolecularReactionTable.hh"
37 #include "G4DNASancheExcitationModel.hh"
38 #include "G4DNAVibExcitation.hh"
39 #include "G4DNAWaterDissociationDisplacer.hh"
40 #include "G4DNAWaterExcitationStructure.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4ProcessManager.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4VDNAReactionModel.hh"
45 // particles
46 #include "G4Electron.hh"
47 #include "G4Electron_aq.hh"
48 #include "G4H2O.hh"
49 #include "G4MoleculeTable.hh"
50 #include "G4O2.hh"
51 #include "G4O3.hh"
52 #include "G4OH.hh"
53 #include "G4Oxygen.hh"
54 #include "G4PhysicsListHelper.hh"
55 /****/
56 #include "G4MolecularConfiguration.hh"
57 #include "G4ProcessTable.hh"
58 /****/
59 // factory
60 #include "G4PhysicsConstructorFactory.hh"
61 #include "G4ChemDissociationChannels_option1.h
62 #include "G4ChemDissociationChannels_option1.h
63 namespace scavenger
64 {
65
66 G4_DECLARE_PHYSCONSTR_FACTORY(EmDNAChemistry);
67
68 //....oooOO0OOooo........oooOO0OOooo........oo
69
70 EmDNAChemistry::EmDNAChemistry()
71 : G4VUserChemistryList(true),
72 G4UImessenger(),
73 fpParserDir(new G4UIdirectory("/chem/react
74 fpReactionTableNameCmd(new G4UIcmdWithAStr
75 {
76 G4DNAChemistryManager::Instance()->SetChemis
77 fpParserDir->SetGuidance("Chemistry control"
78 fpReactionTableNameCmd->SetGuidance("Name of
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 void EmDNAChemistry::SetNewValue(G4UIcommand*
84 {
85 if (command == fpReactionTableNameCmd.get())
86 fReactionTableName = newValue;
87 }
88 }
89
90 //....oooOO0OOooo........oooOO0OOooo........oo
91
92 void EmDNAChemistry::ConstructMolecule()
93 {
94 // Create the definition
95 G4ChemDissociationChannels_option1::Construc
96
97 auto G4NO2 = new G4MoleculeDefinition("NO_2"
98 /*mass
99 /*D*/
100 /*char
101 /*elec
102 /*radi
103
104 auto G4NO3 = new G4MoleculeDefinition("NO_3"
105 /*mass
106 /*D*/
107 /*char
108 /*elec
109 /*radi
110 //__________________________________________
111 // Note: Parameters Value changed according
112
113 G4MoleculeTable::Instance()->CreateConfigura
114
115
116
117
118
119 G4MoleculeTable::Instance()->CreateConfigura
120
121
122
123
124
125 G4MoleculeTable::Instance()->CreateConfigura
126
127
128
129
130
131 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oo
134
135 void EmDNAChemistry::ConstructDissociationChan
136 {
137 //-----------------------------------
138 // Get the molecular configuration
139 auto OH = G4MoleculeTable::Instance()->GetCo
140 auto OHm = G4MoleculeTable::Instance()->GetC
141 auto e_aq = G4MoleculeTable::Instance()->Get
142 auto H2 = G4MoleculeTable::Instance()->GetCo
143 auto H3O = G4MoleculeTable::Instance()->GetC
144 auto H = G4MoleculeTable::Instance()->GetCon
145
146 //-------------------------------------
147 // Define the decay channels
148 auto* water = G4H2O::Definition();
149 G4MolecularDissociationChannel* decCh1;
150 G4MolecularDissociationChannel* decCh2;
151
152 auto occ = new G4ElectronOccupancy(*(water->
153
154 ////////////////////////////////////////////
155 // EXCITATIONS
156 ////////////////////////////////////////////
157 G4DNAWaterExcitationStructure waterExcitatio
158 //------------------------------------------
159 //---------------Excitation on the fifth lay
160
161 decCh1 = new G4MolecularDissociationChannel(
162 decCh2 = new G4MolecularDissociationChannel(
163 // Decay 1 : OH + H
164 decCh1->SetEnergy(waterExcitation.Excitation
165 decCh1->SetProbability(0.35);
166 decCh1->SetDisplacementType(G4DNAWaterDissoc
167
168 decCh2->AddProduct(OH);
169 decCh2->AddProduct(H);
170 decCh2->SetProbability(0.65);
171 decCh2->SetDisplacementType(G4DNAWaterDissoc
172
173 // water->AddExcitedState("A^1B_1");
174 occ->RemoveElectron(4, 1); // this is the t
175 occ->AddElectron(5, 1); // the first unoccu
176
177 water->NewConfigurationWithElectronOccupancy
178 water->AddDecayChannel("A^1B_1", decCh1);
179 water->AddDecayChannel("A^1B_1", decCh2);
180
181 //------------------------------------------
182 //---------------Excitation on the fourth la
183 decCh1 = new G4MolecularDissociationChannel(
184 decCh2 = new G4MolecularDissociationChannel(
185 auto decCh3 = new G4MolecularDissociationCha
186
187 // Decay 1 : energy
188 decCh1->SetEnergy(waterExcitation.Excitation
189 decCh1->SetProbability(0.3);
190
191 // Decay 2 : 2OH + H_2
192 decCh2->AddProduct(H2);
193 decCh2->AddProduct(OH);
194 decCh2->AddProduct(OH);
195 decCh2->SetProbability(0.15);
196 decCh2->SetDisplacementType(G4DNAWaterDissoc
197
198 // Decay 3 : OH + H_3Op + e_aq
199 decCh3->AddProduct(OH);
200 decCh3->AddProduct(H3O);
201 decCh3->AddProduct(e_aq);
202 decCh3->SetProbability(0.55);
203 decCh3->SetDisplacementType(G4DNAWaterDissoc
204
205 *occ = *(water->GetGroundStateElectronOccupa
206 occ->RemoveElectron(3); // this is the tran
207 occ->AddElectron(5, 1); // the first unoccu
208
209 water->NewConfigurationWithElectronOccupancy
210 water->AddDecayChannel("B^1A_1", decCh1);
211 water->AddDecayChannel("B^1A_1", decCh2);
212 water->AddDecayChannel("B^1A_1", decCh3);
213
214 //------------------------------------------
215 //-------------------Excitation of 3rd layer
216 decCh1 = new G4MolecularDissociationChannel(
217 decCh2 = new G4MolecularDissociationChannel(
218
219 // Decay channel 1 : : OH + H_3Op + e_aq
220 decCh1->AddProduct(OH);
221 decCh1->AddProduct(H3O);
222 decCh1->AddProduct(e_aq);
223
224 decCh1->SetProbability(0.5);
225 decCh1->SetDisplacementType(G4DNAWaterDissoc
226
227 // Decay channel 2 : energy
228 decCh2->SetEnergy(waterExcitation.Excitation
229 decCh2->SetProbability(0.5);
230
231 // Electronic configuration of this decay
232 *occ = *(water->GetGroundStateElectronOccupa
233 occ->RemoveElectron(2, 1);
234 occ->AddElectron(5, 1);
235
236 // Configure the water molecule
237 water->NewConfigurationWithElectronOccupancy
238 water->AddDecayChannel("Excitation3rdLayer",
239 water->AddDecayChannel("Excitation3rdLayer",
240
241 //------------------------------------------
242 //-------------------Excitation of 2nd layer
243 decCh1 = new G4MolecularDissociationChannel(
244 decCh2 = new G4MolecularDissociationChannel(
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(G4DNAWaterDissoc
253
254 // Decay channel 2 : energy
255 decCh2->SetEnergy(waterExcitation.Excitation
256 decCh2->SetProbability(0.5);
257
258 *occ = *(water->GetGroundStateElectronOccupa
259 occ->RemoveElectron(1, 1);
260 occ->AddElectron(5, 1);
261
262 water->NewConfigurationWithElectronOccupancy
263 water->AddDecayChannel("Excitation2ndLayer",
264 water->AddDecayChannel("Excitation2ndLayer",
265
266 //------------------------------------------
267 //-------------------Excitation of 1st layer
268 decCh1 = new G4MolecularDissociationChannel(
269 decCh2 = new G4MolecularDissociationChannel(
270
271 *occ = *(water->GetGroundStateElectronOccupa
272 occ->RemoveElectron(0, 1);
273 occ->AddElectron(5, 1);
274
275 // Decay Channel 1 : : OH + H_3Op + e_aq
276 decCh1->AddProduct(OH);
277 decCh1->AddProduct(H3O);
278 decCh1->AddProduct(e_aq);
279 decCh1->SetProbability(0.5);
280 decCh1->SetDisplacementType(G4DNAWaterDissoc
281
282 // Decay channel 2 : energy
283 decCh2->SetEnergy(waterExcitation.Excitation
284 decCh2->SetProbability(0.5);
285
286 water->NewConfigurationWithElectronOccupancy
287 water->AddDecayChannel("Excitation1stLayer",
288 water->AddDecayChannel("Excitation1stLayer",
289
290 ////////////////////////////////////////////
291 // IONISATION
292 ////////////////////////////////////////////
293 //------------------------------------------
294 //------------------- Ionisation -----------
295
296 decCh1 = new G4MolecularDissociationChannel(
297
298 // Decay Channel 1 : : OH + H_3Op
299 decCh1->AddProduct(H3O);
300 decCh1->AddProduct(OH);
301 decCh1->SetProbability(1);
302 decCh1->SetDisplacementType(G4DNAWaterDissoc
303
304 *occ = *(water->GetGroundStateElectronOccupa
305 occ->RemoveElectron(4, 1);
306 // this is a ionized h2O with a hole in its
307 water->NewConfigurationWithElectronOccupancy
308 water->AddDecayChannel("Ionisation5", decCh1
309
310 *occ = *(water->GetGroundStateElectronOccupa
311 occ->RemoveElectron(3, 1);
312 water->NewConfigurationWithElectronOccupancy
313 water->AddDecayChannel("Ionisation4", new G4
314
315 *occ = *(water->GetGroundStateElectronOccupa
316 occ->RemoveElectron(2, 1);
317 water->NewConfigurationWithElectronOccupancy
318 water->AddDecayChannel("Ionisation3", new G4
319
320 *occ = *(water->GetGroundStateElectronOccupa
321 occ->RemoveElectron(1, 1);
322 water->NewConfigurationWithElectronOccupancy
323 water->AddDecayChannel("Ionisation2", new G4
324
325 *occ = *(water->GetGroundStateElectronOccupa
326 occ->RemoveElectron(0, 1);
327 water->NewConfigurationWithElectronOccupancy
328 water->AddDecayChannel("Ionisation1", new G4
329
330 ////////////////////////////////////////////
331 // Dissociative Attachment
332 ////////////////////////////////////////////
333 decCh1 = new G4MolecularDissociationChannel(
334
335 // Decay 1 : 2OH + H_2
336 decCh1->AddProduct(H2);
337 decCh1->AddProduct(OHm);
338 decCh1->AddProduct(OH);
339 decCh1->SetProbability(1);
340 decCh1->SetDisplacementType(G4DNAWaterDissoc
341
342 *occ = *(water->GetGroundStateElectronOccupa
343 occ->AddElectron(5, 1); // H_2O^-
344 water->NewConfigurationWithElectronOccupancy
345 water->AddDecayChannel("DissociativeAttachme
346
347 ////////////////////////////////////////////
348 // Electron-hole recombination
349 ////////////////////////////////////////////
350 decCh1 = new G4MolecularDissociationChannel(
351 decCh2 = new G4MolecularDissociationChannel(
352 decCh3 = new G4MolecularDissociationChannel(
353
354 // Decay 1 : 2OH + H_2
355 decCh1->AddProduct(H2);
356 decCh1->AddProduct(OH);
357 decCh1->AddProduct(OH);
358 decCh1->SetProbability(0.15);
359 decCh1->SetDisplacementType(G4DNAWaterDissoc
360
361 // Decay 2 : OH + H
362 decCh2->AddProduct(OH);
363 decCh2->AddProduct(H);
364 decCh2->SetProbability(0.55);
365 decCh2->SetDisplacementType(G4DNAWaterDissoc
366
367 // Decay 3 : relaxation
368 decCh3->SetProbability(0.30);
369
370 const auto pH2Ovib = G4H2O::Definition()->Ne
371 assert(pH2Ovib != nullptr);
372
373 water->AddDecayChannel(pH2Ovib, decCh1);
374 water->AddDecayChannel(pH2Ovib, decCh2);
375 water->AddDecayChannel(pH2Ovib, decCh3);
376
377 delete occ;
378 }
379
380 //....oooOO0OOooo........oooOO0OOooo........oo
381
382 void EmDNAChemistry::ConstructReactionTable(G4
383 {
384 // Construct chemical reactions from user in
385 ParserChemReaction parser;
386 parser.ReadReactionFile(fReactionTableName);
387
388 auto ListReactant1 = parser.GetListReactant1
389 auto ListReactant2 = parser.GetListReactant2
390 auto ListProduct = parser.GetListProduct();
391 auto ListRate = parser.GetListRate();
392
393 const G4int Ntype = 5;
394 for (size_t i = 0; i < Ntype; i++) {
395 for (size_t j = 0; j < ListReactant1[i].si
396 G4MolecularConfiguration* Reactant1 =
397 G4MoleculeTable::Instance()->GetConfig
398 G4MolecularConfiguration* Reactant2 =
399 G4MoleculeTable::Instance()->GetConfig
400
401 auto pReactionData = new G4DNAMolecularR
402
403 for (size_t k = 0; k < ListProduct[i][j]
404 G4MolecularConfiguration* Product =
405 G4MoleculeTable::Instance()->GetConf
406 pReactionData->AddProduct(Product);
407 }
408 // Reaction type II and IV
409 if (i == 1 || i == 3) {
410 pReactionData->SetReactionType(1);
411 }
412
413 theReactionTable->SetReaction(pReactionD
414 }
415 }
416 }
417
418 //....oooOO0OOooo........oooOO0OOooo........oo
419
420 void EmDNAChemistry::ConstructProcess()
421 {
422 auto ph = G4PhysicsListHelper::GetPhysicsLis
423 G4VProcess* process =
424 G4ProcessTable::GetProcessTable()->FindPro
425 if (process) {
426 auto vibExcitation = (G4DNAVibExcitation*)
427 G4VEmModel* model = vibExcitation->EmModel
428 auto sancheExcitationMod = dynamic_cast<G4
429 if (sancheExcitationMod) {
430 sancheExcitationMod->ExtendLowEnergyLimi
431 }
432 }
433 process = G4ProcessTable::GetProcessTable()-
434
435 if (process == nullptr) {
436 ph->RegisterProcess(new G4DNAElectronSolva
437 G4Electron::Definition
438 }
439 auto* theMoleculeTable = G4MoleculeTable::In
440 G4MoleculeDefinitionIterator iterator = theM
441 iterator.reset();
442 while (iterator()) {
443 G4MoleculeDefinition* moleculeDef = iterat
444 if (moleculeDef == G4H2O::Definition()) {
445 moleculeDef->GetProcessManager()->AddRes
446
447 auto dissociationProcess = new G4DNAMole
448 dissociationProcess->SetDisplacer(molecu
449 // dissociationProcess->SetVerboseLevel(
450 moleculeDef->GetProcessManager()->AddRes
451 }
452 }
453 G4DNAChemistryManager::Instance()->Initializ
454 }
455
456 //....oooOO0OOooo........oooOO0OOooo........oo
457
458 void EmDNAChemistry::ConstructTimeStepModel(G4
459 /*
460 {
461 auto irt = new G4DNAMolecularIRTModel();
462 RegisterTimeStepModel(irt, 0);
463 }
464
465 //....oooOO0OOooo........oooOO0OOooo........oo
466
467 } // namespace scavenger
468