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