Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // ------------------------------------------- 27 // -------------------------------------------------------------- 28 // GEANT 4 - Underground Dark Matter Detecto 28 // GEANT 4 - Underground Dark Matter Detector Advanced Example 29 // 29 // 30 // For information related to this code c 30 // For information related to this code contact: Alex Howard 31 // e-mail: alexander.howard@cern.ch 31 // e-mail: alexander.howard@cern.ch 32 // ------------------------------------------- 32 // -------------------------------------------------------------- 33 // Comments 33 // Comments 34 // 34 // 35 // Underground Advanced 35 // Underground Advanced 36 // by A. Howard and H. Araujo 36 // by A. Howard and H. Araujo 37 // (27th November 2001) 37 // (27th November 2001) 38 // 38 // 39 // PhysicsList program 39 // PhysicsList program 40 // 40 // 41 // Modified: 41 // Modified: 42 // 42 // 43 // 14-02-03 Fix bugs in msc and hIon instancia 43 // 14-02-03 Fix bugs in msc and hIon instanciation + cut per region 44 // 44 // 45 // 05-02-05 AH - changes to G4Decay - added is 45 // 05-02-05 AH - changes to G4Decay - added is not short lived protection 46 // and redefined particles to allow n 46 // and redefined particles to allow non-static creation 47 // i.e. changed construction to G4Mes 47 // i.e. changed construction to G4MesonConstructor, G4BaryonConstructor 48 // << 49 // 23-10-09 LP - migrated EM physics from the << 50 // the new G4Livermore model implemen << 51 // 48 // 52 // ------------------------------------------- 49 // -------------------------------------------------------------- 53 50 54 #include <iomanip> << 55 << 56 #include "DMXPhysicsList.hh" 51 #include "DMXPhysicsList.hh" 57 52 58 #include "globals.hh" 53 #include "globals.hh" 59 #include "G4SystemOfUnits.hh" << 60 #include "G4ProcessManager.hh" 54 #include "G4ProcessManager.hh" 61 #include "G4ProcessVector.hh" 55 #include "G4ProcessVector.hh" 62 56 63 #include "G4ParticleDefinition.hh" 57 #include "G4ParticleDefinition.hh" 64 #include "G4ParticleWithCuts.hh" 58 #include "G4ParticleWithCuts.hh" 65 #include "G4ParticleTypes.hh" 59 #include "G4ParticleTypes.hh" 66 #include "G4ParticleTable.hh" 60 #include "G4ParticleTable.hh" 67 61 68 #include "G4ios.hh" 62 #include "G4ios.hh" 69 #include "G4UserLimits.hh" << 63 #include <iomanip> 70 << 71 #include "G4MesonConstructor.hh" << 72 #include "G4BaryonConstructor.hh" << 73 #include "G4IonConstructor.hh" << 74 #include "G4ShortLivedConstructor.hh" << 75 << 76 #include "DMXMaxTimeCuts.hh" << 77 #include "DMXMinEkineCuts.hh" << 78 #include "G4StepLimiter.hh" << 79 << 80 // gamma << 81 #include "G4PhotoElectricEffect.hh" << 82 #include "G4LivermorePhotoElectricModel.hh" << 83 << 84 #include "G4ComptonScattering.hh" << 85 #include "G4LivermoreComptonModel.hh" << 86 << 87 #include "G4GammaConversion.hh" << 88 #include "G4BetheHeitler5DModel.hh" << 89 << 90 #include "G4RayleighScattering.hh" << 91 #include "G4LivermoreRayleighModel.hh" << 92 << 93 // e- << 94 #include "G4eMultipleScattering.hh" << 95 << 96 #include "G4eIonisation.hh" << 97 #include "G4LivermoreIonisationModel.hh" << 98 << 99 #include "G4eBremsstrahlung.hh" << 100 #include "G4UniversalFluctuation.hh" << 101 << 102 // e+ << 103 #include "G4eIonisation.hh" << 104 #include "G4eBremsstrahlung.hh" << 105 #include "G4eplusAnnihilation.hh" << 106 << 107 // alpha and GenericIon and deuterons, triton, << 108 //muon: << 109 #include "G4MuIonisation.hh" << 110 #include "G4MuBremsstrahlung.hh" << 111 #include "G4MuPairProduction.hh" << 112 #include "G4MuMultipleScattering.hh" << 113 #include "G4MuonMinusCapture.hh" << 114 << 115 //OTHERS: << 116 #include "G4hIonisation.hh" << 117 #include "G4hMultipleScattering.hh" << 118 #include "G4hBremsstrahlung.hh" << 119 #include "G4ionIonisation.hh" << 120 #include "G4IonParametrisedLossModel.hh" << 121 #include "G4NuclearStopping.hh" << 122 << 123 //em process options to allow msc step-limitat << 124 #include "G4EmParameters.hh" << 125 #include "G4VAtomDeexcitation.hh" << 126 #include "G4UAtomicDeexcitation.hh" << 127 #include "G4LossTableManager.hh" << 128 << 129 #include "G4Scintillation.hh" << 130 #include "G4OpAbsorption.hh" << 131 #include "G4OpBoundaryProcess.hh" << 132 #include "G4OpticalParameters.hh" << 133 << 134 // Elastic processes: << 135 #include "G4HadronElasticProcess.hh" << 136 #include "G4ChipsElasticModel.hh" << 137 #include "G4ElasticHadrNucleusHE.hh" << 138 << 139 // Inelastic processes: << 140 #include "G4HadronInelasticProcess.hh" << 141 << 142 // High energy FTFP model and Bertini cascade << 143 #include "G4FTFModel.hh" << 144 #include "G4LundStringFragmentation.hh" << 145 #include "G4ExcitedStringDecay.hh" << 146 #include "G4PreCompoundModel.hh" << 147 #include "G4GeneratorPrecompoundInterface.hh" << 148 #include "G4TheoFSGenerator.hh" << 149 #include "G4CascadeInterface.hh" << 150 << 151 // Cross sections << 152 #include "G4VCrossSectionDataSet.hh" << 153 #include "G4CrossSectionDataSetRegistry.hh" << 154 << 155 #include "G4CrossSectionElastic.hh" << 156 #include "G4CrossSectionInelastic.hh" << 157 #include "G4BGGPionElasticXS.hh" << 158 #include "G4BGGPionInelasticXS.hh" << 159 #include "G4AntiNuclElastic.hh" << 160 << 161 #include "G4CrossSectionInelastic.hh" << 162 #include "G4BGGNucleonInelasticXS.hh" << 163 #include "G4BGGNucleonElasticXS.hh" << 164 #include "G4NeutronInelasticXS.hh" << 165 #include "G4NeutronElasticXS.hh" << 166 #include "G4ComponentAntiNuclNuclearXS.hh" << 167 #include "G4ComponentGGNuclNuclXsc.hh" << 168 #include "G4ComponentGGHadronNucleusXsc.hh" << 169 << 170 #include "G4HadronElastic.hh" << 171 #include "G4NeutronCaptureProcess.hh" << 172 << 173 // Neutron high-precision models: <20 MeV << 174 #include "G4ParticleHPElastic.hh" << 175 #include "G4ParticleHPElasticData.hh" << 176 #include "G4ParticleHPCapture.hh" << 177 #include "G4ParticleHPCaptureData.hh" << 178 #include "G4ParticleHPInelastic.hh" << 179 #include "G4ParticleHPInelasticData.hh" << 180 << 181 // Stopping processes << 182 #include "G4HadronStoppingProcess.hh" << 183 #include "G4HadronicAbsorptionBertini.hh" << 184 #include "G4HadronicAbsorptionFritiof.hh" << 185 64 186 #include "G4HadronicParameters.hh" << 65 #include "G4UserLimits.hh" 187 66 188 #include "G4Decay.hh" << 189 #include "G4RadioactiveDecay.hh" << 190 #include "G4PhysicsListHelper.hh" << 191 #include "G4NuclideTable.hh" << 192 #include "G4NuclearLevelData.hh" << 193 67 194 // Constructor /////////////////////////////// 68 // Constructor ///////////////////////////////////////////////////////////// 195 DMXPhysicsList::DMXPhysicsList() : G4VUserPhys 69 DMXPhysicsList::DMXPhysicsList() : G4VUserPhysicsList() 196 { 70 { 197 71 198 defaultCutValue = 1.0*micrometer; // 72 defaultCutValue = 1.0*micrometer; // 199 cutForGamma = defaultCutValue; 73 cutForGamma = defaultCutValue; 200 cutForElectron = 1.0*nanometer; 74 cutForElectron = 1.0*nanometer; 201 cutForPositron = defaultCutValue; 75 cutForPositron = defaultCutValue; 202 76 203 VerboseLevel = 1; 77 VerboseLevel = 1; 204 OpVerbLevel = 0; 78 OpVerbLevel = 0; 205 79 206 //set a finer grid of the physic tables in o << 207 //former LowEnergy models have 200 bins up t << 208 G4EmParameters* param = G4EmParameters::Inst << 209 param->SetMaxEnergy(100*GeV); << 210 param->SetNumberOfBinsPerDecade(20); << 211 param->SetMscStepLimitType(fMinimal); << 212 param->SetFluo(true); << 213 param->SetPixe(true); << 214 param->SetAuger(true); << 215 << 216 G4EmParameters::Instance()->AddPhysics("Worl << 217 G4DeexPrecoParameters* deex = G4NuclearLevel << 218 deex->SetStoreICLevelData(true); << 219 deex->SetMaxLifeTime(G4NuclideTable::GetInst << 220 /std::log(2.)); << 221 SetVerboseLevel(VerboseLevel); 80 SetVerboseLevel(VerboseLevel); 222 } 81 } 223 82 >> 83 224 // Destructor //////////////////////////////// 84 // Destructor ////////////////////////////////////////////////////////////// 225 DMXPhysicsList::~DMXPhysicsList() 85 DMXPhysicsList::~DMXPhysicsList() 226 {;} 86 {;} 227 87 >> 88 228 // Construct Particles /////////////////////// 89 // Construct Particles ///////////////////////////////////////////////////// 229 void DMXPhysicsList::ConstructParticle() 90 void DMXPhysicsList::ConstructParticle() 230 { 91 { >> 92 231 // In this method, static member functions s 93 // In this method, static member functions should be called 232 // for all particles which you want to use. 94 // for all particles which you want to use. 233 // This ensures that objects of these partic 95 // This ensures that objects of these particle types will be 234 // created in the program. 96 // created in the program. 235 97 236 ConstructMyBosons(); 98 ConstructMyBosons(); 237 ConstructMyLeptons(); 99 ConstructMyLeptons(); 238 ConstructMyHadrons(); 100 ConstructMyHadrons(); 239 ConstructMyShortLiveds(); 101 ConstructMyShortLiveds(); >> 102 240 } 103 } 241 104 >> 105 242 // construct Bosons:////////////////////////// 106 // construct Bosons:///////////////////////////////////////////////////// 243 void DMXPhysicsList::ConstructMyBosons() 107 void DMXPhysicsList::ConstructMyBosons() 244 { 108 { 245 // pseudo-particles 109 // pseudo-particles 246 G4Geantino::GeantinoDefinition(); 110 G4Geantino::GeantinoDefinition(); 247 G4ChargedGeantino::ChargedGeantinoDefinition 111 G4ChargedGeantino::ChargedGeantinoDefinition(); 248 112 249 // gamma 113 // gamma 250 G4Gamma::GammaDefinition(); 114 G4Gamma::GammaDefinition(); 251 115 252 //OpticalPhotons 116 //OpticalPhotons 253 G4OpticalPhoton::OpticalPhotonDefinition(); 117 G4OpticalPhoton::OpticalPhotonDefinition(); 254 118 255 } 119 } 256 120 >> 121 257 // construct Leptons:///////////////////////// 122 // construct Leptons:///////////////////////////////////////////////////// 258 void DMXPhysicsList::ConstructMyLeptons() 123 void DMXPhysicsList::ConstructMyLeptons() 259 { 124 { 260 // leptons 125 // leptons 261 G4Electron::ElectronDefinition(); 126 G4Electron::ElectronDefinition(); 262 G4Positron::PositronDefinition(); 127 G4Positron::PositronDefinition(); 263 G4MuonPlus::MuonPlusDefinition(); 128 G4MuonPlus::MuonPlusDefinition(); 264 G4MuonMinus::MuonMinusDefinition(); 129 G4MuonMinus::MuonMinusDefinition(); 265 130 266 G4NeutrinoE::NeutrinoEDefinition(); 131 G4NeutrinoE::NeutrinoEDefinition(); 267 G4AntiNeutrinoE::AntiNeutrinoEDefinition(); 132 G4AntiNeutrinoE::AntiNeutrinoEDefinition(); 268 G4NeutrinoMu::NeutrinoMuDefinition(); 133 G4NeutrinoMu::NeutrinoMuDefinition(); 269 G4AntiNeutrinoMu::AntiNeutrinoMuDefinition() 134 G4AntiNeutrinoMu::AntiNeutrinoMuDefinition(); 270 } 135 } 271 136 >> 137 >> 138 #include "G4MesonConstructor.hh" >> 139 #include "G4BaryonConstructor.hh" >> 140 #include "G4IonConstructor.hh" >> 141 272 // construct Hadrons:///////////////////////// 142 // construct Hadrons:///////////////////////////////////////////////////// 273 void DMXPhysicsList::ConstructMyHadrons() 143 void DMXPhysicsList::ConstructMyHadrons() 274 { 144 { 275 // mesons 145 // mesons 276 G4MesonConstructor mConstructor; 146 G4MesonConstructor mConstructor; 277 mConstructor.ConstructParticle(); 147 mConstructor.ConstructParticle(); 278 148 279 // baryons 149 // baryons 280 G4BaryonConstructor bConstructor; 150 G4BaryonConstructor bConstructor; 281 bConstructor.ConstructParticle(); 151 bConstructor.ConstructParticle(); 282 152 283 // ions 153 // ions 284 G4IonConstructor iConstructor; 154 G4IonConstructor iConstructor; 285 iConstructor.ConstructParticle(); 155 iConstructor.ConstructParticle(); >> 156 286 } 157 } 287 158 >> 159 288 // construct Shortliveds:///////////////////// 160 // construct Shortliveds:///////////////////////////////////////////////////// 289 void DMXPhysicsList::ConstructMyShortLiveds() 161 void DMXPhysicsList::ConstructMyShortLiveds() 290 { 162 { 291 G4ShortLivedConstructor slConstructor; << 163 // ShortLiveds 292 slConstructor.ConstructParticle(); << 164 ; 293 } 165 } 294 166 >> 167 >> 168 >> 169 295 // Construct Processes /////////////////////// 170 // Construct Processes ////////////////////////////////////////////////////// 296 void DMXPhysicsList::ConstructProcess() 171 void DMXPhysicsList::ConstructProcess() 297 { 172 { >> 173 298 AddTransportation(); 174 AddTransportation(); 299 175 300 ConstructEM(); 176 ConstructEM(); 301 177 302 ConstructOp(); 178 ConstructOp(); 303 179 304 ConstructHad(); 180 ConstructHad(); 305 181 306 ConstructGeneral(); 182 ConstructGeneral(); >> 183 307 } 184 } 308 185 >> 186 309 // Transportation //////////////////////////// 187 // Transportation /////////////////////////////////////////////////////////// >> 188 #include "DMXMaxTimeCuts.hh" >> 189 #include "DMXMinEkineCuts.hh" >> 190 #include "G4StepLimiter.hh" >> 191 310 void DMXPhysicsList::AddTransportation() { 192 void DMXPhysicsList::AddTransportation() { 311 193 312 G4VUserPhysicsList::AddTransportation(); 194 G4VUserPhysicsList::AddTransportation(); 313 195 314 auto particleIterator=GetParticleIterator(); << 196 theParticleIterator->reset(); 315 particleIterator->reset(); << 197 while( (*theParticleIterator)() ){ 316 while( (*particleIterator)() ){ << 198 G4ParticleDefinition* particle = theParticleIterator->value(); 317 G4ParticleDefinition* particle = particleI << 318 G4ProcessManager* pmanager = particle->Get 199 G4ProcessManager* pmanager = particle->GetProcessManager(); 319 G4String particleName = particle->GetParti 200 G4String particleName = particle->GetParticleName(); 320 // time cuts for ONLY neutrons: 201 // time cuts for ONLY neutrons: 321 if(particleName == "neutron") 202 if(particleName == "neutron") 322 pmanager->AddDiscreteProcess(new DMXMaxT 203 pmanager->AddDiscreteProcess(new DMXMaxTimeCuts()); 323 // Energy cuts to kill charged (embedded i 204 // Energy cuts to kill charged (embedded in method) particles: 324 pmanager->AddDiscreteProcess(new DMXMinEki 205 pmanager->AddDiscreteProcess(new DMXMinEkineCuts()); 325 206 326 // Step limit applied to all particles: 207 // Step limit applied to all particles: 327 pmanager->AddDiscreteProcess(new G4StepLim << 208 pmanager->AddProcess(new G4StepLimiter, -1,-1,1); >> 209 328 } 210 } 329 } 211 } 330 212 >> 213 331 // Electromagnetic Processes ///////////////// 214 // Electromagnetic Processes //////////////////////////////////////////////// 332 // all charged particles 215 // all charged particles >> 216 #include "G4MultipleScattering.hh" >> 217 >> 218 // gamma >> 219 #include "G4LowEnergyRayleigh.hh" >> 220 #include "G4LowEnergyPhotoElectric.hh" >> 221 #include "G4LowEnergyCompton.hh" >> 222 #include "G4LowEnergyGammaConversion.hh" >> 223 >> 224 >> 225 // e- >> 226 #include "G4LowEnergyIonisation.hh" >> 227 #include "G4LowEnergyBremsstrahlung.hh" >> 228 >> 229 // e+ >> 230 #include "G4eIonisation.hh" >> 231 #include "G4eBremsstrahlung.hh" >> 232 #include "G4eplusAnnihilation.hh" >> 233 >> 234 >> 235 // alpha and GenericIon and deuterons, triton, He3: >> 236 //hIonisation #include "G4hLowEnergyIonisation.hh" -> moved to G4hIonisation >> 237 #include "G4EnergyLossTables.hh" >> 238 // hLowEnergyIonisation uses Ziegler 1988 as the default >> 239 >> 240 >> 241 //muon: >> 242 #include "G4MuIonisation.hh" >> 243 #include "G4MuBremsstrahlung.hh" >> 244 #include "G4MuPairProduction.hh" >> 245 #include "G4MuonMinusCaptureAtRest.hh" >> 246 >> 247 //OTHERS: >> 248 #include "G4hIonisation.hh" // standard hadron ionisation >> 249 >> 250 //em process options to allow msc step-limitation to be switched off >> 251 #include "G4EmProcessOptions.hh" >> 252 333 void DMXPhysicsList::ConstructEM() { 253 void DMXPhysicsList::ConstructEM() { 334 << 335 G4LossTableManager* man = G4LossTableManager << 336 man->SetAtomDeexcitation(new G4UAtomicDeexci << 337 254 338 G4EmParameters* em_params = G4EmParameters:: << 255 // processes: >> 256 >> 257 G4LowEnergyPhotoElectric* lowePhot = new G4LowEnergyPhotoElectric(); >> 258 G4LowEnergyIonisation* loweIon = new G4LowEnergyIonisation(); >> 259 G4LowEnergyBremsstrahlung* loweBrem = new G4LowEnergyBremsstrahlung(); >> 260 >> 261 // note LowEIon uses proton as basis for its data-base, therefore >> 262 // cannot specify different LowEnergyIonisation models for different >> 263 // particles, but can change model globally for Ion, Alpha and Proton. >> 264 >> 265 >> 266 //fluorescence apply specific cut for fluorescence from photons, electrons >> 267 //and bremsstrahlung photons: >> 268 G4double fluorcut = 250*eV; >> 269 lowePhot->SetCutForLowEnSecPhotons(fluorcut); >> 270 loweIon->SetCutForLowEnSecPhotons(fluorcut); >> 271 loweBrem->SetCutForLowEnSecPhotons(fluorcut); 339 272 340 auto particleIterator=GetParticleIterator(); << 273 // setting tables explicitly for electronic stopping power 341 particleIterator->reset(); << 274 // ahadronLowEIon->SetElectronicStoppingPowerModel 342 while( (*particleIterator)() ){ << 275 // (G4GenericIon::GenericIonDefinition(), "ICRU_R49p") ; 343 G4ParticleDefinition* particle = particleI << 276 // ahadronLowEIon->SetElectronicStoppingPowerModel >> 277 // (G4Proton::ProtonDefinition(), "ICRU_R49p") ; >> 278 >> 279 // Switch off the Barkas and Bloch corrections >> 280 // ahadronLowEIon->SetBarkasOff(); >> 281 >> 282 >> 283 theParticleIterator->reset(); >> 284 while( (*theParticleIterator)() ){ >> 285 G4ParticleDefinition* particle = theParticleIterator->value(); 344 G4ProcessManager* pmanager = particle->Get 286 G4ProcessManager* pmanager = particle->GetProcessManager(); 345 G4String particleName = particle->GetParti 287 G4String particleName = particle->GetParticleName(); 346 G4String particleType = particle->GetParti 288 G4String particleType = particle->GetParticleType(); 347 G4double charge = particle->GetPDGCharge() 289 G4double charge = particle->GetPDGCharge(); 348 290 349 if (particleName == "gamma") 291 if (particleName == "gamma") 350 { 292 { 351 //gamma 293 //gamma 352 G4RayleighScattering* theRayleigh = new G4Ra << 294 pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh()); 353 pmanager->AddDiscreteProcess(theRayleigh); << 295 pmanager->AddDiscreteProcess(lowePhot); 354 << 296 pmanager->AddDiscreteProcess(new G4LowEnergyCompton()); 355 G4PhotoElectricEffect* thePhotoElectricEffec << 297 pmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion()); 356 thePhotoElectricEffect->SetEmModel(new G4Liv << 357 pmanager->AddDiscreteProcess(thePhotoElectri << 358 << 359 G4ComptonScattering* theComptonScattering = << 360 theComptonScattering->SetEmModel(new G4Liver << 361 pmanager->AddDiscreteProcess(theComptonScatt << 362 << 363 G4GammaConversion* theGammaConversion = new << 364 theGammaConversion->SetEmModel(new G4BetheHe << 365 pmanager->AddDiscreteProcess(theGammaConvers << 366 << 367 } 298 } 368 else if (particleName == "e-") 299 else if (particleName == "e-") 369 { 300 { 370 //electron 301 //electron 371 // process ordering: AddProcess(name, at res 302 // process ordering: AddProcess(name, at rest, along step, post step) 372 // Multiple scattering << 303 // -1 = not implemented, then ordering 373 G4eMultipleScattering* msc = new G4eMultiple << 304 G4MultipleScattering* aMultipleScattering = new G4MultipleScattering(); 374 em_params->SetMscStepLimitType(fUseDistanceT << 305 pmanager->AddProcess(aMultipleScattering, -1, 1, 1); 375 pmanager->AddProcess(msc,-1, 1, -1); << 306 pmanager->AddProcess(loweIon, -1, 2, 2); 376 << 307 pmanager->AddProcess(loweBrem, -1,-1, 3); 377 // Ionisation << 378 G4eIonisation* eIonisation = new G4eIonisati << 379 G4VEmModel* theIoniLiv = new G4Livermo << 380 theIoniLiv->SetHighEnergyLimit(0.1*MeV << 381 eIonisation->AddEmModel(0, theIoniLiv, << 382 em_params->SetStepFunction(0.2, 100*um); //i << 383 pmanager->AddProcess(eIonisation,-1, 2, 1); << 384 << 385 // Bremsstrahlung << 386 G4eBremsstrahlung* eBremsstrahlung = new G4e << 387 pmanager->AddProcess(eBremsstrahlung, -1,-3, << 388 } 308 } 389 else if (particleName == "e+") 309 else if (particleName == "e+") 390 { 310 { 391 //positron << 311 //positron 392 G4eMultipleScattering* msc = new G4eMultiple << 312 G4MultipleScattering* aMultipleScattering = new G4MultipleScattering(); 393 msc->SetStepLimitType(fUseDistanceToBoundary << 313 pmanager->AddProcess(aMultipleScattering, -1, 1, 1); 394 pmanager->AddProcess(msc,-1, 1, -1); << 314 pmanager->AddProcess(new G4eIonisation(), -1, 2, 2); 395 << 315 pmanager->AddProcess(new G4eBremsstrahlung(), -1,-1, 3); 396 // Ionisation << 316 pmanager->AddProcess(new G4eplusAnnihilation(),0,-1, 4); 397 G4eIonisation* eIonisation = new G4eIonisati << 398 // eIonisation->SetStepFunction(0.2, 100*um) << 399 pmanager->AddProcess(eIonisation, << 400 << 401 //Bremsstrahlung (use default, no low-energy << 402 pmanager->AddProcess(new G4eBremsstrahlung() << 403 << 404 //Annihilation << 405 pmanager->AddProcess(new G4eplusAnnihilation << 406 } 317 } 407 else if( particleName == "mu+" || 318 else if( particleName == "mu+" || 408 particleName == "mu-" ) 319 particleName == "mu-" ) 409 { 320 { 410 //muon 321 //muon 411 pmanager->AddProcess(new G4MuMultipleScatter << 322 G4MultipleScattering* aMultipleScattering = new G4MultipleScattering(); 412 pmanager->AddProcess(new G4MuIonisation(), << 323 pmanager->AddProcess(aMultipleScattering, -1, 1, 1); 413 pmanager->AddProcess(new G4MuBremsstrahlung( << 324 pmanager->AddProcess(new G4MuIonisation(), -1, 2, 2); 414 pmanager->AddProcess(new G4MuPairProduction( << 325 pmanager->AddProcess(new G4MuBremsstrahlung(), -1,-1, 3); >> 326 pmanager->AddProcess(new G4MuPairProduction(), -1,-1, 4); 415 if( particleName == "mu-" ) 327 if( particleName == "mu-" ) 416 pmanager->AddProcess(new G4MuonMinusCaptur << 328 pmanager->AddProcess(new G4MuonMinusCaptureAtRest(), 0,-1,-1); 417 } 329 } 418 else if (particleName == "proton" || << 330 else if (particleName == "proton" || 419 particleName == "pi+" || << 331 particleName == "alpha" || 420 particleName == "pi-") << 421 { << 422 //multiple scattering << 423 pmanager->AddProcess(new G4hMultipleScatteri << 424 << 425 //ionisation << 426 G4hIonisation* hIonisation = new G4hIonisati << 427 em_params->SetStepFunctionMuHad(0.2, 50*um); << 428 pmanager->AddProcess(hIonisation, << 429 << 430 //bremmstrahlung << 431 pmanager->AddProcess(new G4hBremsstrahlung, << 432 } << 433 else if(particleName == "alpha" || << 434 particleName == "deuteron" || 332 particleName == "deuteron" || 435 particleName == "triton" || 333 particleName == "triton" || 436 particleName == "He3") << 334 particleName == "He3" || 437 { << 335 particleName == "GenericIon" || 438 //multiple scattering << 336 (particleType == "nucleus" && charge != 0)) 439 pmanager->AddProcess(new G4hMultipleScatteri << 440 << 441 //ionisation << 442 G4ionIonisation* ionIoni = new G4ionIonisati << 443 em_params->SetStepFunctionLightIons(0.1, 1*C << 444 pmanager->AddProcess(ionIoni, << 445 pmanager->AddProcess(new G4NuclearStopping() << 446 } << 447 else if (particleName == "GenericIon") << 448 { 337 { 449 // OBJECT may be dynamically created as eith 338 // OBJECT may be dynamically created as either a GenericIon or nucleus 450 // G4Nucleus exists and therefore has partic 339 // G4Nucleus exists and therefore has particle type nucleus 451 // genericIon: 340 // genericIon: 452 << 341 G4MultipleScattering* aMultipleScattering = new G4MultipleScattering(); 453 //multiple scattering << 342 //hIonisation G4hLowEnergyIonisation* ahadronLowEIon = new G4hLowEnergyIonisation(); 454 pmanager->AddProcess(new G4hMultipleScatteri << 343 G4hIonisation* ahadronIon = new G4hIonisation(); 455 << 344 pmanager->AddProcess(aMultipleScattering,-1,1,1); 456 //ionisation << 345 //hIonisation pmanager->AddProcess(ahadronLowEIon,-1,2,2); 457 G4ionIonisation* ionIoni = new G4ionIonisati << 346 pmanager->AddProcess(ahadronIon,-1,2,2); 458 em_params->SetStepFunctionIons(0.1, 1*CLHEP: << 347 // ahadronLowEIon->SetNuclearStoppingOff() ; 459 pmanager->AddProcess(ionIoni, << 348 // ahadronLowEIon->SetNuclearStoppingPowerModel("ICRU_R49") ; 460 pmanager->AddProcess(new G4NuclearStopping() << 349 // ahadronLowEIon->SetNuclearStoppingOn() ; >> 350 >> 351 //fluorescence switch off for hadrons (for now) PIXE: >> 352 //hIonisation ahadronLowEIon->SetFluorescence(false); 461 } 353 } 462 << 463 else if ((!particle->IsShortLived()) && 354 else if ((!particle->IsShortLived()) && 464 (charge != 0.0) && 355 (charge != 0.0) && 465 (particle->GetParticleName() != "charge 356 (particle->GetParticleName() != "chargedgeantino")) 466 { 357 { 467 //all others charged particles except geanti 358 //all others charged particles except geantino 468 G4hMultipleScattering* aMultipleScatte << 359 G4MultipleScattering* aMultipleScattering = new G4MultipleScattering(); >> 360 //hIonisation G4hLowEnergyIonisation* ahadronLowEIon = new G4hLowEnergyIonisation(); 469 G4hIonisation* ahadronIon = new G4hIon 361 G4hIonisation* ahadronIon = new G4hIonisation(); 470 << 362 pmanager->AddProcess(aMultipleScattering,-1,1,1); 471 //multiple scattering << 363 //hIonisation pmanager->AddProcess(ahadronLowEIon, -1,2,2); 472 pmanager->AddProcess(aMultipleScattering,-1, << 364 pmanager->AddProcess(ahadronIon, -1,2,2); 473 << 365 // pmanager->AddProcess(new G4hIonisation(), -1,2,2); 474 //ionisation << 475 pmanager->AddProcess(ahadronIon, -1,2, << 476 } 366 } >> 367 477 } 368 } >> 369 >> 370 // turn off msc step-limitation - especially as electron cut 1nm >> 371 G4EmProcessOptions opt; >> 372 // opt.SetMscStepLimitation(false); >> 373 opt.SetMscStepLimitation(fMinimal); >> 374 478 } 375 } 479 376 >> 377 480 // Optical Processes ///////////////////////// 378 // Optical Processes //////////////////////////////////////////////////////// >> 379 #include "G4Scintillation.hh" >> 380 #include "G4OpAbsorption.hh" >> 381 //#include "G4OpRayleigh.hh" >> 382 #include "G4OpBoundaryProcess.hh" >> 383 481 void DMXPhysicsList::ConstructOp() 384 void DMXPhysicsList::ConstructOp() 482 { 385 { 483 G4OpticalParameters* opParams = G4OpticalPar << 386 // default scintillation process 484 G4Scintillation* theScintProcessDef = new G4 387 G4Scintillation* theScintProcessDef = new G4Scintillation("Scintillation"); 485 opParams->SetScintTrackSecondariesFirst(true << 388 // theScintProcessDef->DumpPhysicsTable(); 486 opParams->SetScintByParticleType(true); << 389 theScintProcessDef->SetTrackSecondariesFirst(true); >> 390 theScintProcessDef->SetScintillationYieldFactor(1.0); // >> 391 theScintProcessDef->SetScintillationExcitationRatio(0.0); // >> 392 theScintProcessDef->SetVerboseLevel(OpVerbLevel); >> 393 >> 394 // scintillation process for alpha: >> 395 G4Scintillation* theScintProcessAlpha = new G4Scintillation("Scintillation"); >> 396 // theScintProcessNuc->DumpPhysicsTable(); >> 397 theScintProcessAlpha->SetTrackSecondariesFirst(true); >> 398 theScintProcessAlpha->SetScintillationYieldFactor(1.1); >> 399 theScintProcessAlpha->SetScintillationExcitationRatio(1.0); >> 400 theScintProcessAlpha->SetVerboseLevel(OpVerbLevel); >> 401 >> 402 // scintillation process for heavy nuclei >> 403 G4Scintillation* theScintProcessNuc = new G4Scintillation("Scintillation"); >> 404 // theScintProcessNuc->DumpPhysicsTable(); >> 405 theScintProcessNuc->SetTrackSecondariesFirst(true); >> 406 theScintProcessNuc->SetScintillationYieldFactor(0.2); >> 407 theScintProcessNuc->SetScintillationExcitationRatio(1.0); >> 408 theScintProcessNuc->SetVerboseLevel(OpVerbLevel); 487 409 488 // optical processes 410 // optical processes 489 G4OpAbsorption* theAbsorptionProcess = new G 411 G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption(); >> 412 // G4OpRayleigh* theRayleighScatteringProcess = new G4OpRayleigh(); 490 G4OpBoundaryProcess* theBoundaryProcess = ne 413 G4OpBoundaryProcess* theBoundaryProcess = new G4OpBoundaryProcess(); >> 414 // theAbsorptionProcess->DumpPhysicsTable(); >> 415 // theRayleighScatteringProcess->DumpPhysicsTable(); >> 416 theAbsorptionProcess->SetVerboseLevel(OpVerbLevel); >> 417 // theRayleighScatteringProcess->SetVerboseLevel(OpVerbLevel); >> 418 theBoundaryProcess->SetVerboseLevel(OpVerbLevel); >> 419 G4OpticalSurfaceModel themodel = unified; >> 420 theBoundaryProcess->SetModel(themodel); 491 421 492 auto particleIterator=GetParticleIterator(); << 422 theParticleIterator->reset(); 493 particleIterator->reset(); << 423 while( (*theParticleIterator)() ) 494 while( (*particleIterator)() ) << 495 { 424 { 496 G4ParticleDefinition* particle = particl << 425 G4ParticleDefinition* particle = theParticleIterator->value(); 497 G4ProcessManager* pmanager = particle->G 426 G4ProcessManager* pmanager = particle->GetProcessManager(); 498 G4String particleName = particle->GetPar 427 G4String particleName = particle->GetParticleName(); 499 if (theScintProcessDef->IsApplicable(*pa 428 if (theScintProcessDef->IsApplicable(*particle)) { 500 pmanager->AddProcess(theScintProcessDe << 429 // if(particle->GetPDGMass() > 5.0*GeV) 501 pmanager->SetProcessOrderingToLast(the << 430 if(particle->GetParticleName() == "GenericIon") { 502 pmanager->SetProcessOrderingToLast(the << 431 pmanager->AddProcess(theScintProcessNuc); // AtRestDiscrete >> 432 pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxAtRest); >> 433 pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxPostStep); >> 434 } >> 435 else if(particle->GetParticleName() == "alpha") { >> 436 pmanager->AddProcess(theScintProcessAlpha); >> 437 pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxAtRest); >> 438 pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxPostStep); >> 439 } >> 440 else { >> 441 pmanager->AddProcess(theScintProcessDef); >> 442 pmanager->SetProcessOrderingToLast(theScintProcessDef,idxAtRest); >> 443 pmanager->SetProcessOrderingToLast(theScintProcessDef,idxPostStep); >> 444 } 503 } 445 } 504 446 505 if (particleName == "opticalphoton") { 447 if (particleName == "opticalphoton") { 506 pmanager->AddDiscreteProcess(theAbsorptionPr 448 pmanager->AddDiscreteProcess(theAbsorptionProcess); >> 449 // pmanager->AddDiscreteProcess(theRayleighScatteringProcess); 507 pmanager->AddDiscreteProcess(theBoundaryProc 450 pmanager->AddDiscreteProcess(theBoundaryProcess); 508 } 451 } 509 } 452 } 510 } 453 } 511 454 >> 455 512 // Hadronic processes //////////////////////// 456 // Hadronic processes //////////////////////////////////////////////////////// 513 457 >> 458 // Elastic processes: >> 459 #include "G4HadronElasticProcess.hh" >> 460 >> 461 // Inelastic processes: >> 462 #include "G4PionPlusInelasticProcess.hh" >> 463 #include "G4PionMinusInelasticProcess.hh" >> 464 #include "G4KaonPlusInelasticProcess.hh" >> 465 #include "G4KaonZeroSInelasticProcess.hh" >> 466 #include "G4KaonZeroLInelasticProcess.hh" >> 467 #include "G4KaonMinusInelasticProcess.hh" >> 468 #include "G4ProtonInelasticProcess.hh" >> 469 #include "G4AntiProtonInelasticProcess.hh" >> 470 #include "G4NeutronInelasticProcess.hh" >> 471 #include "G4AntiNeutronInelasticProcess.hh" >> 472 #include "G4DeuteronInelasticProcess.hh" >> 473 #include "G4TritonInelasticProcess.hh" >> 474 #include "G4AlphaInelasticProcess.hh" >> 475 >> 476 // Low-energy Models: < 20GeV >> 477 #include "G4LElastic.hh" >> 478 #include "G4LEPionPlusInelastic.hh" >> 479 #include "G4LEPionMinusInelastic.hh" >> 480 #include "G4LEKaonPlusInelastic.hh" >> 481 #include "G4LEKaonZeroSInelastic.hh" >> 482 #include "G4LEKaonZeroLInelastic.hh" >> 483 #include "G4LEKaonMinusInelastic.hh" >> 484 #include "G4LEProtonInelastic.hh" >> 485 #include "G4LEAntiProtonInelastic.hh" >> 486 #include "G4LENeutronInelastic.hh" >> 487 #include "G4LEAntiNeutronInelastic.hh" >> 488 #include "G4LEDeuteronInelastic.hh" >> 489 #include "G4LETritonInelastic.hh" >> 490 #include "G4LEAlphaInelastic.hh" >> 491 >> 492 // High-energy Models: >20 GeV >> 493 #include "G4HEPionPlusInelastic.hh" >> 494 #include "G4HEPionMinusInelastic.hh" >> 495 #include "G4HEKaonPlusInelastic.hh" >> 496 #include "G4HEKaonZeroInelastic.hh" >> 497 #include "G4HEKaonZeroInelastic.hh" >> 498 #include "G4HEKaonMinusInelastic.hh" >> 499 #include "G4HEProtonInelastic.hh" >> 500 #include "G4HEAntiProtonInelastic.hh" >> 501 #include "G4HENeutronInelastic.hh" >> 502 #include "G4HEAntiNeutronInelastic.hh" >> 503 >> 504 // Neutron high-precision models: <20 MeV >> 505 #include "G4NeutronHPElastic.hh" >> 506 #include "G4NeutronHPElasticData.hh" >> 507 #include "G4NeutronHPCapture.hh" >> 508 #include "G4NeutronHPCaptureData.hh" >> 509 #include "G4NeutronHPInelastic.hh" >> 510 #include "G4NeutronHPInelasticData.hh" >> 511 #include "G4LCapture.hh" >> 512 >> 513 // Stopping processes >> 514 #include "G4PiMinusAbsorptionAtRest.hh" >> 515 #include "G4KaonMinusAbsorptionAtRest.hh" >> 516 #include "G4AntiProtonAnnihilationAtRest.hh" >> 517 #include "G4AntiNeutronAnnihilationAtRest.hh" >> 518 >> 519 >> 520 // ConstructHad() >> 521 // Makes discrete physics processes for the hadrons, at present limited >> 522 // to those particles with GHEISHA interactions (INTRC > 0). >> 523 // The processes are: Elastic scattering and Inelastic scattering. >> 524 // F.W.Jones 09-JUL-1998 514 void DMXPhysicsList::ConstructHad() 525 void DMXPhysicsList::ConstructHad() 515 { 526 { 516 //Elastic models << 527 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess; 517 G4HadronElastic* elastic_lhep0 = new G4Hadro << 528 G4LElastic* theElasticModel = new G4LElastic; 518 G4ChipsElasticModel* elastic_chip = new G4Ch << 529 theElasticProcess->RegisterMe(theElasticModel); 519 G4ElasticHadrNucleusHE* elastic_he = new G4E << 520 530 521 // Inelastic scattering << 531 theParticleIterator->reset(); 522 const G4double theFTFMin0 = 0.0*GeV; << 532 while ((*theParticleIterator)()) 523 const G4double theFTFMin1 = 3.0*GeV; << 524 const G4double theFTFMax = G4HadronicParamet << 525 const G4double theBERTMin0 = 0.0*GeV; << 526 const G4double theBERTMin1 = 19.0*MeV; << 527 const G4double theBERTMax = 6.0*GeV; << 528 const G4double theHPMin = 0.0*GeV; << 529 const G4double theHPMax = 20.0*MeV; << 530 << 531 G4FTFModel * theStringModel = new G4FTFModel << 532 G4ExcitedStringDecay * theStringDecay = new << 533 theStringModel->SetFragmentationModel( theSt << 534 G4PreCompoundModel * thePreEquilib = new G4P << 535 G4GeneratorPrecompoundInterface * theCascade << 536 << 537 G4TheoFSGenerator * theFTFModel0 = new G4The << 538 theFTFModel0->SetHighEnergyGenerator( theStr << 539 theFTFModel0->SetTransport( theCascade ); << 540 theFTFModel0->SetMinEnergy( theFTFMin0 ); << 541 theFTFModel0->SetMaxEnergy( theFTFMax ); << 542 << 543 G4TheoFSGenerator * theFTFModel1 = new G4The << 544 theFTFModel1->SetHighEnergyGenerator( theStr << 545 theFTFModel1->SetTransport( theCascade ); << 546 theFTFModel1->SetMinEnergy( theFTFMin1 ); << 547 theFTFModel1->SetMaxEnergy( theFTFMax ); << 548 << 549 G4CascadeInterface * theBERTModel0 = new G4C << 550 theBERTModel0->SetMinEnergy( theBERTMin0 ); << 551 theBERTModel0->SetMaxEnergy( theBERTMax ); << 552 << 553 G4CascadeInterface * theBERTModel1 = new G4C << 554 theBERTModel1->SetMinEnergy( theBERTMin1 ); << 555 theBERTModel1->SetMaxEnergy( theBERTMax ); << 556 << 557 G4VCrossSectionDataSet * theAntiNucleonData << 558 G4ComponentGGNuclNuclXsc * ggNuclNuclXsec = << 559 G4VCrossSectionDataSet * theGGNuclNuclData = << 560 G4VCrossSectionDataSet * theGGNNEl = new G4C << 561 G4ComponentGGHadronNucleusXsc * ggHNXsec = n << 562 G4VCrossSectionDataSet * theGGHNEl = new G4C << 563 G4VCrossSectionDataSet * theGGHNInel = new G << 564 << 565 auto particleIterator=GetParticleIterator(); << 566 particleIterator->reset(); << 567 while ((*particleIterator)()) << 568 { 533 { 569 G4ParticleDefinition* particle = particl << 534 G4ParticleDefinition* particle = theParticleIterator->value(); 570 G4ProcessManager* pmanager = particle->G 535 G4ProcessManager* pmanager = particle->GetProcessManager(); 571 G4String particleName = particle->GetPar 536 G4String particleName = particle->GetParticleName(); 572 537 573 if (particleName == "pi+") 538 if (particleName == "pi+") 574 { 539 { 575 // Elastic scattering << 540 pmanager->AddDiscreteProcess(theElasticProcess); 576 G4HadronElasticProcess* theElasticPr << 541 G4PionPlusInelasticProcess* theInelasticProcess = 577 theElasticProcess->AddDataSet( new G << 542 new G4PionPlusInelasticProcess("inelastic"); 578 theElasticProcess->RegisterMe( elast << 543 G4LEPionPlusInelastic* theLEInelasticModel = 579 pmanager->AddDiscreteProcess( theElasticPr << 544 new G4LEPionPlusInelastic; 580 //Inelastic scattering << 545 theInelasticProcess->RegisterMe(theLEInelasticModel); 581 G4HadronInelasticProcess* theInelasticProc << 546 G4HEPionPlusInelastic* theHEInelasticModel = 582 new G4HadronInelasticProcess( "inelastic << 547 new G4HEPionPlusInelastic; 583 theInelasticProcess->AddDataSet( new G4BGG << 548 theInelasticProcess->RegisterMe(theHEInelasticModel); 584 theInelasticProcess->RegisterMe( theFTFMod << 549 pmanager->AddDiscreteProcess(theInelasticProcess); 585 theInelasticProcess->RegisterMe( the << 586 pmanager->AddDiscreteProcess( theInelastic << 587 } 550 } 588 551 589 else if (particleName == "pi-") 552 else if (particleName == "pi-") 590 { 553 { 591 // Elastic scattering << 554 pmanager->AddDiscreteProcess(theElasticProcess); 592 G4HadronElasticProcess* theElasticPr << 555 G4PionMinusInelasticProcess* theInelasticProcess = 593 theElasticProcess->AddDataSet( new G << 556 new G4PionMinusInelasticProcess("inelastic"); 594 theElasticProcess->RegisterMe( elast << 557 G4LEPionMinusInelastic* theLEInelasticModel = 595 pmanager->AddDiscreteProcess( theElasticPr << 558 new G4LEPionMinusInelastic; 596 //Inelastic scattering << 559 theInelasticProcess->RegisterMe(theLEInelasticModel); 597 G4HadronInelasticProcess* theInelasticProc << 560 G4HEPionMinusInelastic* theHEInelasticModel = 598 new G4HadronInelasticProcess( "inelastic << 561 new G4HEPionMinusInelastic; 599 theInelasticProcess->AddDataSet( new G4BGG << 562 theInelasticProcess->RegisterMe(theHEInelasticModel); 600 theInelasticProcess->RegisterMe( theFTFMod << 563 pmanager->AddDiscreteProcess(theInelasticProcess); 601 theInelasticProcess->RegisterMe( the << 564 G4String prcNam; 602 pmanager->AddDiscreteProcess( theInelastic << 565 pmanager->AddRestProcess(new G4PiMinusAbsorptionAtRest, ordDefault); 603 //Absorption << 604 pmanager->AddRestProcess(new G4HadronicAbs << 605 } 566 } >> 567 606 else if (particleName == "kaon+") 568 else if (particleName == "kaon+") 607 { 569 { 608 // Elastic scattering << 570 pmanager->AddDiscreteProcess(theElasticProcess); 609 G4HadronElasticProcess* theElasticPr << 571 G4KaonPlusInelasticProcess* theInelasticProcess = 610 theElasticProcess->AddDataSet( theGGHNEl ) << 572 new G4KaonPlusInelasticProcess("inelastic"); 611 theElasticProcess->RegisterMe( elast << 573 G4LEKaonPlusInelastic* theLEInelasticModel = 612 pmanager->AddDiscreteProcess( theElasticPr << 574 new G4LEKaonPlusInelastic; 613 // Inelastic scattering << 575 theInelasticProcess->RegisterMe(theLEInelasticModel); 614 G4HadronInelasticProcess* theInelasticProc << 576 G4HEKaonPlusInelastic* theHEInelasticModel = 615 new G4HadronInelasticProcess( "inelastic << 577 new G4HEKaonPlusInelastic; 616 theInelasticProcess->AddDataSet( theGGHNIn << 578 theInelasticProcess->RegisterMe(theHEInelasticModel); 617 theInelasticProcess->RegisterMe( theFTFMod << 579 pmanager->AddDiscreteProcess(theInelasticProcess); 618 theInelasticProcess->RegisterMe( the << 580 } 619 pmanager->AddDiscreteProcess( theInelastic << 581 620 } << 621 else if (particleName == "kaon0S") 582 else if (particleName == "kaon0S") 622 { 583 { 623 // Elastic scattering << 584 pmanager->AddDiscreteProcess(theElasticProcess); 624 G4HadronElasticProcess* theElasticPr << 585 G4KaonZeroSInelasticProcess* theInelasticProcess = 625 theElasticProcess->AddDataSet( theGGHNEl ) << 586 new G4KaonZeroSInelasticProcess("inelastic"); 626 theElasticProcess->RegisterMe( elast << 587 G4LEKaonZeroSInelastic* theLEInelasticModel = 627 pmanager->AddDiscreteProcess( theElasticPr << 588 new G4LEKaonZeroSInelastic; 628 // Inelastic scattering << 589 theInelasticProcess->RegisterMe(theLEInelasticModel); 629 G4HadronInelasticProcess* theInelasticProc << 590 G4HEKaonZeroInelastic* theHEInelasticModel = 630 new G4HadronInelasticProcess( "inelastic << 591 new G4HEKaonZeroInelastic; 631 theInelasticProcess->AddDataSet( theGGHNIn << 592 theInelasticProcess->RegisterMe(theHEInelasticModel); 632 theInelasticProcess->RegisterMe( theFTFMod << 593 pmanager->AddDiscreteProcess(theInelasticProcess); 633 theInelasticProcess->RegisterMe( the << 634 pmanager->AddDiscreteProcess( theInelastic << 635 } 594 } 636 595 637 else if (particleName == "kaon0L") 596 else if (particleName == "kaon0L") 638 { 597 { 639 // Elastic scattering << 598 pmanager->AddDiscreteProcess(theElasticProcess); 640 G4HadronElasticProcess* theElasticPr << 599 G4KaonZeroLInelasticProcess* theInelasticProcess = 641 theElasticProcess->AddDataSet( theGGHNEl ) << 600 new G4KaonZeroLInelasticProcess("inelastic"); 642 theElasticProcess->RegisterMe( elast << 601 G4LEKaonZeroLInelastic* theLEInelasticModel = 643 pmanager->AddDiscreteProcess( theElasticPr << 602 new G4LEKaonZeroLInelastic; 644 // Inelastic scattering << 603 theInelasticProcess->RegisterMe(theLEInelasticModel); 645 G4HadronInelasticProcess* theInelasticProc << 604 G4HEKaonZeroInelastic* theHEInelasticModel = 646 new G4HadronInelasticProcess( "inelastic << 605 new G4HEKaonZeroInelastic; 647 theInelasticProcess->AddDataSet( theGGHNIn << 606 theInelasticProcess->RegisterMe(theHEInelasticModel); 648 theInelasticProcess->RegisterMe( theFTFMod << 607 pmanager->AddDiscreteProcess(theInelasticProcess); 649 theInelasticProcess->RegisterMe( the << 650 pmanager->AddDiscreteProcess( theInelastic << 651 } 608 } 652 609 653 else if (particleName == "kaon-") 610 else if (particleName == "kaon-") 654 { 611 { 655 // Elastic scattering << 612 pmanager->AddDiscreteProcess(theElasticProcess); 656 G4HadronElasticProcess* theElasticPr << 613 G4KaonMinusInelasticProcess* theInelasticProcess = 657 theElasticProcess->AddDataSet( theGGHNEl ) << 614 new G4KaonMinusInelasticProcess("inelastic"); 658 theElasticProcess->RegisterMe( elast << 615 G4LEKaonMinusInelastic* theLEInelasticModel = 659 pmanager->AddDiscreteProcess( theElasticPr << 616 new G4LEKaonMinusInelastic; 660 // Inelastic scattering << 617 theInelasticProcess->RegisterMe(theLEInelasticModel); 661 G4HadronInelasticProcess* theInelasticProc << 618 G4HEKaonMinusInelastic* theHEInelasticModel = 662 new G4HadronInelasticProcess( "inelastic << 619 new G4HEKaonMinusInelastic; 663 theInelasticProcess->AddDataSet( the << 620 theInelasticProcess->RegisterMe(theHEInelasticModel); 664 theInelasticProcess->RegisterMe( theFTFMod << 621 pmanager->AddDiscreteProcess(theInelasticProcess); 665 theInelasticProcess->RegisterMe( the << 622 pmanager->AddRestProcess(new G4KaonMinusAbsorptionAtRest, ordDefault); 666 pmanager->AddDiscreteProcess( theInelastic << 667 pmanager->AddRestProcess(new G4HadronicAbs << 668 } 623 } 669 624 670 else if (particleName == "proton") 625 else if (particleName == "proton") 671 { 626 { 672 // Elastic scattering << 627 pmanager->AddDiscreteProcess(theElasticProcess); 673 G4HadronElasticProcess* theElasticPr << 628 G4ProtonInelasticProcess* theInelasticProcess = 674 theElasticProcess->AddDataSet( new G << 629 new G4ProtonInelasticProcess("inelastic"); 675 theElasticProcess->RegisterMe( elast << 630 G4LEProtonInelastic* theLEInelasticModel = new G4LEProtonInelastic; 676 pmanager->AddDiscreteProcess( theElasticPr << 631 theInelasticProcess->RegisterMe(theLEInelasticModel); 677 // Inelastic scattering << 632 G4HEProtonInelastic* theHEInelasticModel = new G4HEProtonInelastic; 678 G4HadronInelasticProcess* theInelasticProc << 633 theInelasticProcess->RegisterMe(theHEInelasticModel); 679 new G4HadronInelasticProcess( "inelastic << 634 pmanager->AddDiscreteProcess(theInelasticProcess); 680 theInelasticProcess->AddDataSet( new G4BGG << 681 theInelasticProcess->RegisterMe( theFTFMod << 682 theInelasticProcess->RegisterMe( the << 683 pmanager->AddDiscreteProcess( theInelastic << 684 } 635 } >> 636 685 else if (particleName == "anti_proton") 637 else if (particleName == "anti_proton") 686 { 638 { 687 // Elastic scattering << 639 pmanager->AddDiscreteProcess(theElasticProcess); 688 const G4double elastic_elimitAntiNuc << 640 G4AntiProtonInelasticProcess* theInelasticProcess = 689 G4AntiNuclElastic* elastic_anuc = ne << 641 new G4AntiProtonInelasticProcess("inelastic"); 690 elastic_anuc->SetMinEnergy( elastic_ << 642 G4LEAntiProtonInelastic* theLEInelasticModel = 691 G4CrossSectionElastic* elastic_anucx << 643 new G4LEAntiProtonInelastic; 692 G4HadronElastic* elastic_lhep2 = new << 644 theInelasticProcess->RegisterMe(theLEInelasticModel); 693 elastic_lhep2->SetMaxEnergy( elastic << 645 G4HEAntiProtonInelastic* theHEInelasticModel = 694 G4HadronElasticProcess* theElasticPr << 646 new G4HEAntiProtonInelastic; 695 theElasticProcess->AddDataSet( elast << 647 theInelasticProcess->RegisterMe(theHEInelasticModel); 696 theElasticProcess->RegisterMe( elast << 648 pmanager->AddDiscreteProcess(theInelasticProcess); 697 theElasticProcess->RegisterMe( elast << 698 pmanager->AddDiscreteProcess( theElasticPr << 699 // Inelastic scattering << 700 G4HadronInelasticProcess* theInelasticProc << 701 new G4HadronInelasticProcess( "inelastic << 702 theInelasticProcess->AddDataSet( theAntiNu << 703 theInelasticProcess->RegisterMe( theFTFMod << 704 pmanager->AddDiscreteProcess( theInelastic << 705 // Absorption << 706 pmanager->AddRestProcess(new G4HadronicAbs << 707 } 649 } >> 650 708 else if (particleName == "neutron") { 651 else if (particleName == "neutron") { 709 // elastic scattering 652 // elastic scattering 710 G4HadronElasticProcess* theElasticProcess = << 653 G4HadronElasticProcess* theNeutronElasticProcess = 711 theElasticProcess->AddDataSet(new G4Ne << 654 new G4HadronElasticProcess; 712 G4HadronElastic* elastic_neutronChipsM << 655 G4LElastic* theElasticModel1 = new G4LElastic; 713 elastic_neutronChipsModel->SetMinEnergy( 19. << 656 G4NeutronHPElastic * theElasticNeutron = new G4NeutronHPElastic; 714 theElasticProcess->RegisterMe( elastic << 657 theNeutronElasticProcess->RegisterMe(theElasticModel1); 715 G4ParticleHPElastic * theElasticNeutronHP = << 658 theElasticModel1->SetMinEnergy(19*MeV); 716 theElasticNeutronHP->SetMinEnergy( the << 659 theNeutronElasticProcess->RegisterMe(theElasticNeutron); 717 theElasticNeutronHP->SetMaxEnergy( the << 660 G4NeutronHPElasticData * theNeutronData = new G4NeutronHPElasticData; 718 theElasticProcess->RegisterMe( theElasticNeu << 661 theNeutronElasticProcess->AddDataSet(theNeutronData); 719 theElasticProcess->AddDataSet( new G4Particl << 662 pmanager->AddDiscreteProcess(theNeutronElasticProcess); 720 pmanager->AddDiscreteProcess( theElasticProc << 663 // inelastic scattering 721 // inelastic scattering << 664 G4NeutronInelasticProcess* theInelasticProcess = 722 G4HadronInelasticProcess* theInelasticProces << 665 new G4NeutronInelasticProcess("inelastic"); 723 new G4HadronInelasticProcess( "inelastic", << 666 G4LENeutronInelastic* theInelasticModel = new G4LENeutronInelastic; 724 theInelasticProcess->AddDataSet( new G4Neutr << 667 theInelasticModel->SetMinEnergy(19*MeV); 725 theInelasticProcess->RegisterMe( theFTFModel << 668 theInelasticProcess->RegisterMe(theInelasticModel); 726 theInelasticProcess->RegisterMe( theBE << 669 G4NeutronHPInelastic * theLENeutronInelasticModel = 727 G4ParticleHPInelastic * theNeutronInelasticH << 670 new G4NeutronHPInelastic; 728 theNeutronInelasticHPModel->SetMinEner << 671 theInelasticProcess->RegisterMe(theLENeutronInelasticModel); 729 theNeutronInelasticHPModel->SetMaxEner << 672 G4NeutronHPInelasticData * theNeutronData1 = 730 theInelasticProcess->RegisterMe( theNeutronI << 673 new G4NeutronHPInelasticData; 731 theInelasticProcess->AddDataSet( new G4Parti << 674 theInelasticProcess->AddDataSet(theNeutronData1); 732 pmanager->AddDiscreteProcess(theInelasticPro 675 pmanager->AddDiscreteProcess(theInelasticProcess); 733 // capture 676 // capture 734 G4NeutronCaptureProcess* theCaptureProcess = << 677 G4HadronCaptureProcess* theCaptureProcess = 735 new G4NeutronCaptureProcess; << 678 new G4HadronCaptureProcess; 736 G4ParticleHPCapture * theLENeutronCaptureMod << 679 G4LCapture* theCaptureModel = new G4LCapture; 737 theLENeutronCaptureModel->SetMinEnergy(theHP << 680 theCaptureModel->SetMinEnergy(19*MeV); 738 theLENeutronCaptureModel->SetMaxEnergy(theHP << 681 theCaptureProcess->RegisterMe(theCaptureModel); >> 682 G4NeutronHPCapture * theLENeutronCaptureModel = new G4NeutronHPCapture; 739 theCaptureProcess->RegisterMe(theLENeutronCa 683 theCaptureProcess->RegisterMe(theLENeutronCaptureModel); 740 theCaptureProcess->AddDataSet( new G4Particl << 684 G4NeutronHPCaptureData * theNeutronData3 = new G4NeutronHPCaptureData; >> 685 theCaptureProcess->AddDataSet(theNeutronData3); 741 pmanager->AddDiscreteProcess(theCaptureProce 686 pmanager->AddDiscreteProcess(theCaptureProcess); >> 687 // G4ProcessManager* pmanager = G4Neutron::Neutron->GetProcessManager(); >> 688 // pmanager->AddProcess(new G4UserSpecialCuts(),-1,-1,1); 742 } 689 } 743 else if (particleName == "anti_neutron") 690 else if (particleName == "anti_neutron") 744 { 691 { 745 // Elastic scattering << 692 pmanager->AddDiscreteProcess(theElasticProcess); 746 G4HadronElasticProcess* theElasticPr << 693 G4AntiNeutronInelasticProcess* theInelasticProcess = 747 theElasticProcess->AddDataSet( theGGHNEl ) << 694 new G4AntiNeutronInelasticProcess("inelastic"); 748 theElasticProcess->RegisterMe( elast << 695 G4LEAntiNeutronInelastic* theLEInelasticModel = 749 pmanager->AddDiscreteProcess( theElasticPr << 696 new G4LEAntiNeutronInelastic; 750 // Inelastic scattering (include ann << 697 theInelasticProcess->RegisterMe(theLEInelasticModel); 751 G4HadronInelasticProcess* theInelasticProc << 698 G4HEAntiNeutronInelastic* theHEInelasticModel = 752 new G4HadronInelasticProcess( "inelastic << 699 new G4HEAntiNeutronInelastic; 753 theInelasticProcess->AddDataSet( theAntiNu << 700 theInelasticProcess->RegisterMe(theHEInelasticModel); 754 theInelasticProcess->RegisterMe( theFTFMod << 701 pmanager->AddDiscreteProcess(theInelasticProcess); 755 pmanager->AddDiscreteProcess( theInelastic << 756 } 702 } >> 703 757 else if (particleName == "deuteron") 704 else if (particleName == "deuteron") 758 { 705 { 759 // Elastic scattering << 706 pmanager->AddDiscreteProcess(theElasticProcess); 760 G4HadronElasticProcess* theElasticPr << 707 G4DeuteronInelasticProcess* theInelasticProcess = 761 theElasticProcess->AddDataSet( theGGNNEl ) << 708 new G4DeuteronInelasticProcess("inelastic"); 762 theElasticProcess->RegisterMe( elast << 709 G4LEDeuteronInelastic* theLEInelasticModel = 763 pmanager->AddDiscreteProcess( theElasticPr << 710 new G4LEDeuteronInelastic; 764 // Inelastic scattering << 711 theInelasticProcess->RegisterMe(theLEInelasticModel); 765 G4HadronInelasticProcess* theInelasticProc << 712 pmanager->AddDiscreteProcess(theInelasticProcess); 766 new G4HadronInelasticProcess( "inelastic << 767 theInelasticProcess->AddDataSet( theGGNucl << 768 theInelasticProcess->RegisterMe( theFTFMod << 769 theInelasticProcess->RegisterMe( the << 770 pmanager->AddDiscreteProcess( theInelastic << 771 } 713 } >> 714 772 else if (particleName == "triton") 715 else if (particleName == "triton") 773 { 716 { 774 // Elastic scattering << 717 pmanager->AddDiscreteProcess(theElasticProcess); 775 G4HadronElasticProcess* theElasticPr << 718 G4TritonInelasticProcess* theInelasticProcess = 776 theElasticProcess->AddDataSet( theGGNNEl ) << 719 new G4TritonInelasticProcess("inelastic"); 777 theElasticProcess->RegisterMe( elast << 720 G4LETritonInelastic* theLEInelasticModel = 778 pmanager->AddDiscreteProcess( theElasticPr << 721 new G4LETritonInelastic; 779 // Inelastic scattering << 722 theInelasticProcess->RegisterMe(theLEInelasticModel); 780 G4HadronInelasticProcess* theInelasticProc << 723 pmanager->AddDiscreteProcess(theInelasticProcess); 781 new G4HadronInelasticProcess( "inelastic << 782 theInelasticProcess->AddDataSet( theGGNucl << 783 theInelasticProcess->RegisterMe( theFTFMod << 784 theInelasticProcess->RegisterMe( the << 785 pmanager->AddDiscreteProcess( theInelastic << 786 } 724 } >> 725 787 else if (particleName == "alpha") 726 else if (particleName == "alpha") 788 { 727 { 789 // Elastic scattering << 728 pmanager->AddDiscreteProcess(theElasticProcess); 790 G4HadronElasticProcess* theElasticPr << 729 G4AlphaInelasticProcess* theInelasticProcess = 791 theElasticProcess->AddDataSet( theGGNNEl ) << 730 new G4AlphaInelasticProcess("inelastic"); 792 theElasticProcess->RegisterMe( elast << 731 G4LEAlphaInelastic* theLEInelasticModel = 793 pmanager->AddDiscreteProcess( theElasticPr << 732 new G4LEAlphaInelastic; 794 // Inelastic scattering << 733 theInelasticProcess->RegisterMe(theLEInelasticModel); 795 G4HadronInelasticProcess* theInelasticProc << 734 pmanager->AddDiscreteProcess(theInelasticProcess); 796 new G4HadronInelasticProcess( "inelastic << 797 theInelasticProcess->AddDataSet( the << 798 theInelasticProcess->RegisterMe( theFTFMod << 799 theInelasticProcess->RegisterMe( the << 800 pmanager->AddDiscreteProcess( theInelastic << 801 } 735 } >> 736 802 } 737 } 803 } 738 } 804 739 >> 740 805 // Decays //////////////////////////////////// 741 // Decays /////////////////////////////////////////////////////////////////// >> 742 #include "G4Decay.hh" >> 743 #include "G4RadioactiveDecay.hh" >> 744 #include "G4IonTable.hh" >> 745 #include "G4Ions.hh" >> 746 806 void DMXPhysicsList::ConstructGeneral() { 747 void DMXPhysicsList::ConstructGeneral() { 807 748 808 // Add Decay Process 749 // Add Decay Process 809 G4Decay* theDecayProcess = new G4Decay(); 750 G4Decay* theDecayProcess = new G4Decay(); 810 auto particleIterator=GetParticleIterator(); << 751 theParticleIterator->reset(); 811 particleIterator->reset(); << 752 while( (*theParticleIterator)() ) 812 while( (*particleIterator)() ) << 813 { 753 { 814 G4ParticleDefinition* particle = particl << 754 G4ParticleDefinition* particle = theParticleIterator->value(); 815 G4ProcessManager* pmanager = particle->G 755 G4ProcessManager* pmanager = particle->GetProcessManager(); 816 756 817 if (theDecayProcess->IsApplicable(*parti 757 if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived()) 818 { 758 { 819 pmanager ->AddProcess(theDecayProcess); 759 pmanager ->AddProcess(theDecayProcess); 820 // set ordering for PostStepDoIt and AtRes 760 // set ordering for PostStepDoIt and AtRestDoIt 821 pmanager ->SetProcessOrdering(theDecayProc 761 pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep); 822 pmanager ->SetProcessOrdering(theDecayProc 762 pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest); 823 } 763 } 824 } 764 } 825 765 826 // Declare radioactive decay to the GenericI 766 // Declare radioactive decay to the GenericIon in the IonTable. 827 G4LossTableManager* man = G4LossTableManager << 767 const G4IonTable *theIonTable = 828 G4VAtomDeexcitation* ad = man->AtomDeexcitat << 768 G4ParticleTable::GetParticleTable()->GetIonTable(); 829 if(!ad) { << 769 G4RadioactiveDecay *theRadioactiveDecay = new G4RadioactiveDecay(); 830 G4EmParameters::Instance()->SetAugerCascad << 831 ad = new G4UAtomicDeexcitation(); << 832 man->SetAtomDeexcitation(ad); << 833 ad->InitialiseAtomicDeexcitation(); << 834 } << 835 770 836 G4PhysicsListHelper::GetPhysicsListHelper()- << 771 for (G4int i=0; i<theIonTable->Entries(); i++) 837 RegisterProcess(new G4RadioactiveDecay(), << 772 { >> 773 G4String particleName = theIonTable->GetParticle(i)->GetParticleName(); >> 774 G4String particleType = theIonTable->GetParticle(i)->GetParticleType(); >> 775 >> 776 if (particleName == "GenericIon") >> 777 { >> 778 G4ProcessManager* pmanager = >> 779 theIonTable->GetParticle(i)->GetProcessManager(); >> 780 pmanager->SetVerboseLevel(VerboseLevel); >> 781 pmanager ->AddProcess(theRadioactiveDecay); >> 782 pmanager ->SetProcessOrdering(theRadioactiveDecay, idxPostStep); >> 783 pmanager ->SetProcessOrdering(theRadioactiveDecay, idxAtRest); >> 784 } >> 785 } 838 } 786 } 839 787 840 // Cuts ////////////////////////////////////// 788 // Cuts ///////////////////////////////////////////////////////////////////// 841 void DMXPhysicsList::SetCuts() 789 void DMXPhysicsList::SetCuts() 842 { 790 { 843 791 844 if (verboseLevel >1) 792 if (verboseLevel >1) 845 G4cout << "DMXPhysicsList::SetCuts:"; 793 G4cout << "DMXPhysicsList::SetCuts:"; 846 794 847 if (verboseLevel>0){ 795 if (verboseLevel>0){ 848 G4cout << "DMXPhysicsList::SetCuts:"; 796 G4cout << "DMXPhysicsList::SetCuts:"; 849 G4cout << "CutLength : " 797 G4cout << "CutLength : " 850 << G4BestUnit(defaultCutValue,"Length") < 798 << G4BestUnit(defaultCutValue,"Length") << G4endl; 851 } 799 } 852 800 853 //special for low energy physics 801 //special for low energy physics 854 G4double lowlimit=250*eV; 802 G4double lowlimit=250*eV; 855 G4ProductionCutsTable::GetProductionCutsTabl 803 G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(lowlimit,100.*GeV); 856 804 857 // set cut values for gamma at first and for 805 // set cut values for gamma at first and for e- second and next for e+, 858 // because some processes for e+/e- need cut 806 // because some processes for e+/e- need cut values for gamma 859 SetCutValue(cutForGamma, "gamma"); 807 SetCutValue(cutForGamma, "gamma"); 860 SetCutValue(cutForElectron, "e-"); 808 SetCutValue(cutForElectron, "e-"); 861 SetCutValue(cutForPositron, "e+"); 809 SetCutValue(cutForPositron, "e+"); 862 810 863 if (verboseLevel>0) DumpCutValuesTable(); 811 if (verboseLevel>0) DumpCutValuesTable(); 864 } 812 } 865 813 866 814