application/src/eRositaPhysicsList.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026
0027 #include "eRositaPhysicsList.hh"
0028
0029 #include "globals.hh"
0030 #include "G4ios.hh"
0031 #include "G4ParticleDefinition.hh"
0032 #include "G4ParticleTypes.hh"
0033 #include "G4PhysicsListHelper.hh"
0034 #include "G4ProductionCutsTable.hh"
0035 #include "G4StepLimiter.hh"
0036 #include "G4SystemOfUnits.hh"
0037 #include "G4VPhysicsConstructor.hh"
0038
0039 // Physics List
0040 #include "G4DecayPhysics.hh"
0041 #include "G4EmLivermorePhysics.hh"
0042 #include "G4EmPenelopePhysics.hh"
0043 #include "G4EmStandardPhysics.hh"
0044
0045 // Process
0046 #include "G4ComptonScattering.hh"
0047 #include "G4eBremsstrahlung.hh"
0048 #include "G4eIonisation.hh"
0049 #include "G4eMultipleScattering.hh"
0050 #include "G4eplusAnnihilation.hh"
0051 #include "G4GammaConversion.hh"
0052 #include "G4hImpactIonisation.hh"
0053 #include "G4hIonisation.hh"
0054 #include "G4hMultipleScattering.hh"
0055 #include "G4ionIonisation.hh"
0056 #include "G4PhotoElectricEffect.hh"
0057 #include "G4RayleighScattering.hh"
0058
0059 // Model
0060 #include "G4LivermoreBremsstrahlungModel.hh"
0061 #include "G4LivermoreComptonModel.hh"
0062 #include "G4LivermoreGammaConversionModel.hh"
0063 #include "G4LivermoreIonisationModel.hh"
0064 #include "G4LivermorePhotoElectricModel.hh"
0065 #include "G4LivermoreRayleighModel.hh"
0066 #include "G4PenelopeAnnihilationModel.hh"
0067 #include "G4UniversalFluctuation.hh"
0068
0069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0070
0071 eRositaPhysicsList::eRositaPhysicsList()
0072 {
0073     SetVerboseLevel(1);
0074
0075     constexpr auto DEFAULT_CUT_VALUE{0.001 * mm};
0076     SetDefaultCutValue(DEFAULT_CUT_VALUE);
0077
0078     G4cout << "============================================================================================="
0079               << G4endl
0080               << "Geant4 eRosita example, based on a simplified version of eROSITA simulation."
0081               << G4endl
0082               << "Further details can be found in:"
0083               << G4endl
0084               << " M. G. Pia et al.,"
0085               << G4endl
0086               << "  'PIXE Simulation With Geant4',"
0087               << G4endl
0088               << "  IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp. 3614-3649, 2009"
0089               << G4endl
0090               << " N. Meidinger et al.,"
0091               << G4endl
0092               << "  'Development of the focal plane PNCCD camera system for the X-ray space telescope eROSITA',"
0093               << G4endl
0094               << "  NIM A 624, 321-329, 2010"
0095               << G4endl
0096               << "============================================================================================="
0097               << G4endl;
0098
0099     G4cout << G4endl;
0100 }
0101
0102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0103
0104 eRositaPhysicsList::~eRositaPhysicsList()
0105 {
0106 }
0107
0108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0109
0110 void eRositaPhysicsList::ConstructBosons()
0111 {
0112     // geantino (pseudo-particle)
0113     // G4Geantino::GeantinoDefinition();
0114
0115     // charged geantino (pseudo-particle)
0116     // G4ChargedGeantino::ChargedGeantinoDefinition();
0117
0118     // photon (gamma)
0119     G4Gamma::GammaDefinition();
0120 }
0121
0122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0123
0124 void eRositaPhysicsList::ConstructLeptons()
0125 {
0126     // leptons
0127
0128     // e+ / e-
0129     G4Electron::ElectronDefinition();
0130     G4Positron::PositronDefinition();
0131
0132     // mu+ / mu-
0133     // G4MuonPlus::MuonPlusDefinition();
0134     // G4MuonMinus::MuonMinusDefinition();
0135
0136     // nu_e
0137     // G4NeutrinoE::NeutrinoEDefinition();
0138     // G4AntiNeutrinoE::AntiNeutrinoEDefinition();
0139
0140     // nu_mu
0141     // G4NeutrinoMu::NeutrinoMuDefinition();
0142     // G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
0143 }
0144
0145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0146
0147 void eRositaPhysicsList::ConstructMesons()
0148 {
0149     // light mesons
0150
0151     // pion
0152     // G4PionPlus::PionPlusDefinition();
0153     // G4PionMinus::PionMinusDefinition();
0154     // G4PionZero::PionZeroDefinition();
0155
0156     // eta
0157     // G4Eta::EtaDefinition();
0158     // G4EtaPrime::EtaPrimeDefinition();
0159
0160     // kaon
0161     // G4KaonPlus::KaonPlusDefinition();
0162     // G4KaonMinus::KaonMinusDefinition();
0163     // G4KaonZero::KaonZeroDefinition();
0164     // G4AntiKaonZero::AntiKaonZeroDefinition();
0165     // G4KaonZeroLong::KaonZeroLongDefinition();
0166     // G4KaonZeroShort::KaonZeroShortDefinition();
0167 }
0168
0169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0170
0171 void eRositaPhysicsList::ConstructBaryons()
0172 {
0173     // baryons
0174
0175     // proton
0176     G4Proton::ProtonDefinition();
0177     G4AntiProton::AntiProtonDefinition();
0178
0179     // neutron
0180     // G4Neutron::NeutronDefinition();
0181     // G4AntiNeutron::AntiNeutronDefinition();
0182 }
0183
0184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0185
0186 void eRositaPhysicsList::ConstructParticle()
0187 {
0188     ConstructBosons();
0189     ConstructLeptons();
0190     ConstructMesons();
0191     ConstructBaryons();
0192 }
0193
0194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0195
0196 void eRositaPhysicsList::ConstructEM()
0197 {
0198     auto *helper = G4PhysicsListHelper::GetPhysicsListHelper();
0199
0200     auto particleIterator = GetParticleIterator();
0201     particleIterator->reset();
0202
0203     while ((*particleIterator)()) {
0204         auto *particle = particleIterator->value();
0205         auto particleName = particle->GetParticleName();
0206
0207         if (particleName == "gamma") { // photon
0208             // photoelectric effect
0209             auto *photoelectricEffect = new G4PhotoElectricEffect();
0210             // photoelectricEffect->ActivateAuger(true);
0211             // photoelectricEffect->SetCutForLowEnSecElectrons(0.250 * keV);
0212             // photoelectricEffect->SetCutForLowEnSecPhotons(0.250 * keV);
0213             photoelectricEffect->SetEmModel(new G4LivermorePhotoElectricModel()); // Set the model (Livermore)
0214             helper->RegisterProcess(photoelectricEffect, particle);
0215
0216             // Compton scattering
0217             auto *comptonScattering = new G4ComptonScattering();
0218             comptonScattering->SetEmModel(new G4LivermoreComptonModel()); // Set the model (Livermore)
0219             helper->RegisterProcess(comptonScattering, particle);
0220
0221             // gamma conversion
0222             auto *gammaConversion = new G4GammaConversion();
0223             gammaConversion->SetEmModel(new G4LivermoreGammaConversionModel()); // Set the model (Livermore)
0224             helper->RegisterProcess(gammaConversion, particle);
0225
0226             // Rayleigh scattering
0227             auto *rayleighScattering = new G4RayleighScattering();
0228             rayleighScattering->SetEmModel(new G4LivermoreRayleighModel()); // Set the model (Livermore)
0229             helper->RegisterProcess(rayleighScattering, particle);
0230         } else if (particleName == "e-") { // electron
0231             // multiple scattering
0232             helper->RegisterProcess(new G4eMultipleScattering(), particle);
0233
0234             // ionization
0235             auto *ionization = new G4eIonisation();
0236             ionization->SetEmModel(new G4LivermoreIonisationModel()); // Set the model (Livermore)
0237             ionization->SetFluctModel(new G4UniversalFluctuation());
0238             helper->RegisterProcess(ionization, particle);
0239
0240             // Bremsstrahlung
0241             auto *bremsstrahlung = new G4eBremsstrahlung();
0242             bremsstrahlung->SetEmModel(new G4LivermoreBremsstrahlungModel()); // Set the model (Livermore)
0243             helper->RegisterProcess(bremsstrahlung, particle);
0244         } else if (particleName == "e+") { // positron
0245             // multiple scattering
0246             helper->RegisterProcess(new G4eMultipleScattering(), particle);
0247
0248             // ionization
0249             helper->RegisterProcess(new G4eIonisation(), particle);
0250
0251             // Bremsstrahlung
0252             auto *bremsstrahlung = new G4eBremsstrahlung();
0253             bremsstrahlung->SetEmModel(new G4LivermoreBremsstrahlungModel()); // Set the model (Livermore)
0254             helper->RegisterProcess(bremsstrahlung, particle);
0255
0256             // annihilation
0257             auto *annihilation = new G4eplusAnnihilation();
0258             annihilation->SetEmModel(new G4PenelopeAnnihilationModel()); // Set the model (Penelope)
0259             helper->RegisterProcess(annihilation, particle);
0260         // } else if( particleName == "mu+" || particleName == "mu-") {
0261             // // muon
0262             // helper->RegisterProcess(new G4MuMultipleScattering, particle);
0263             // helper->RegisterProcess(new G4MuIonisation, particle);
0264             // helper->RegisterProcess(new G4MuBremsstrahlung, particle);
0265             // helper->RegisterProcess(new G4MuPairProduction, particle);
0266         } else if (particleName == "proton" || particleName == "pi-" || particleName == "pi+") {
0267             helper->RegisterProcess(new G4hMultipleScattering(), particle);
0268             helper->RegisterProcess(new G4hIonisation(), particle);
0269 /*
0270             // proton
0271             // auto *ionization = new G4hImpactIonisation();
0272             // ionization->SetPixeCrossSectionK("ecpssr");
0273             // ionization->SetPixeCrossSectionL("ecpssr");
0274             // ionization->SetPixeCrossSectionM("ecpssr");
0275             // ionization->SetPixeProjectileMinEnergy(1.* keV);
0276             // ionization->SetPixeProjectileMaxEnergy(200. * MeV);
0277             // ionization->SetCutForSecondaryPhotons(250. * eV);
0278             // ionization->SetCutForAugerElectrons(250. * eV);
0279
0280             auto *ionization = new G4hIonisation();
0281             auto *multipleScattering = new G4hMultipleScattering();
0282
0283             processManager->AddProcess(multipleScattering, -1, 1, 1);
0284             processManager->AddProcess(ionization, -1, 2, 2);
0285 */
0286         } else if (particleName == "alpha" || particleName == "He3" || particleName == "pi-" || particleName == "pi+" || particleName == "GenericIon") {
0287             // pion, alpha, ion (should never occur in the current example)
0288             helper->RegisterProcess(new G4hMultipleScattering, particle);
0289             helper->RegisterProcess(new G4ionIonisation, particle);
0290         } else if ((!particle->IsShortLived()) && (particle->GetPDGCharge() != 0.0) && (particle->GetParticleName() != "chargedgeantino")) {
0291             // every other charged particle, except geantino
0292             helper->RegisterProcess(new G4hMultipleScattering, particle);
0293             helper->RegisterProcess(new G4hIonisation, particle);
0294         }
0295     }
0296 }
0297
0298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0299
0300 void eRositaPhysicsList::ConstructGeneral()
0301 {
0302     auto *helper = G4PhysicsListHelper::GetPhysicsListHelper();
0303
0304     // Add decay process
0305     auto *decay = new G4Decay();
0306
0307     auto particleIterator = GetParticleIterator();
0308     particleIterator->reset();
0309
0310     while ((*particleIterator)()) {
0311         auto *particle = particleIterator->value();
0312
0313         if (decay->IsApplicable(*particle)) {
0314             if (verboseLevel > 1) {
0315                 G4cout << "### Decays for " << particle->GetParticleName() << G4endl;
0316             }
0317             helper->RegisterProcess(decay, particle);
0318 /*
0319             // Set ordering for PostStepDoIt and AtRestDoIt
0320             processManager->SetProcessOrdering(decay, idxPostStep);
0321             processManager->SetProcessOrdering(decay, idxAtRest);
0322 */
0323         }
0324     }
0325 }
0326
0327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0328
0329 void eRositaPhysicsList::ConstructProcess()
0330 {
0331     AddTransportation();
0332     ConstructEM();
0333     ConstructGeneral();
0334     // AddStepMax();
0335 }
0336
0337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0338
0339 void eRositaPhysicsList::SetCuts()
0340 {
0341     // Set the default cut value for all particle types
0342     SetCutsWithDefault();
0343
0344     // Set the secondary production cut lower than 990. eV. Very important for processes at low energies.
0345     constexpr auto ENERGY_LOW_LIMIT{250. * eV};
0346     constexpr auto ENERGY_HIGH_LIMIT{100. * GeV};
0347
0348     G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(ENERGY_LOW_LIMIT, ENERGY_HIGH_LIMIT);
0349
0350     if (verboseLevel > 0) {
0351         DumpCutValuesTable();
0352     }
0353 }
0354
0355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0356
0357 /*
0358 #include "G4StepLimiter.hh"
0359 #include "G4UserSpecialCuts.hh"
0360
0361 void eRositaPhysicsList::AddStepMax()
0362 {
0363     auto *helper = G4PhysicsListHelper::GetPhysicsListHelper();
0364
0365     // Step limitation seen as a process
0366     // auto *stepLimiter = new G4StepLimiter();
0367     // // auto *userCuts = new G4UserSpecialCuts();
0368
0369     particleIterator->reset();
0370
0371     while ((*particleIterator)()){
0372         auto *particle = particleIterator->value();
0373         // auto *processManager = particle->GetProcessManager();
0374
0375         if (particle->GetPDGCharge() != 0.0) {
0376             helper->RegisterProcess(stepLimiter, particle);
0377             // helper->RegisterProcess(userCuts, particle);
0378         }
0379     }
0380 }
0381 */