splitting/src/PhysicsList.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 /// \file PhysicsList.cc
0028 /// \brief Implementation of the PhysicsList class.
0029
0030 #include "PhysicsList.hh"
0031
0032 #include "G4SystemOfUnits.hh"
0033 #include "G4UIcmdWithAnInteger.hh"
0034
0035 // Geant4-DNA MODELS
0036
0037 #include "G4DNAAttachment.hh"
0038 #include "G4DNABornExcitationModel.hh"
0039 #include "G4DNABornIonisationModel.hh"
0040 #include "G4DNAChampionElasticModel.hh"
0041 #include "G4DNAChargeDecrease.hh"
0042 #include "G4DNAChargeIncrease.hh"
0043 #include "G4DNADingfelderChargeDecreaseModel.hh"
0044 #include "G4DNADingfelderChargeIncreaseModel.hh"
0045 #include "G4DNAElastic.hh"
0046 #include "G4DNAElectronSolvation.hh"
0047 #include "G4DNAExcitation.hh"
0048 #include "G4DNAIonisation.hh"
0049 #include "G4DNAMeltonAttachmentModel.hh"
0050 #include "G4DNAMillerGreenExcitationModel.hh"
0051 #include "G4DNARuddIonisationModel.hh"
0052 #include "G4DNASancheExcitationModel.hh"
0053 #include "G4DNAScreenedRutherfordElasticModel.hh"
0054 #include "G4DNAVibExcitation.hh"
0055
0056 //
0057
0058 #include "SplitProcess.hh"
0059
0060 #include "G4BetheBlochIonGasModel.hh"
0061 #include "G4BraggIonGasModel.hh"
0062 #include "G4DummyModel.hh"
0063 #include "G4ElectronCapture.hh"
0064 #include "G4EmConfigurator.hh"
0065 #include "G4IonFluctuations.hh"
0066 #include "G4LossTableManager.hh"
0067 #include "G4MollerBhabhaModel.hh"
0068 #include "G4UniversalFluctuation.hh"
0069 #include "G4UrbanMscModel.hh"
0070 #include "G4VEmModel.hh"
0071 #include "G4eIonisation.hh"
0072 #include "G4eMultipleScattering.hh"
0073 #include "G4hIonisation.hh"
0074 #include "G4hMultipleScattering.hh"
0075 #include "G4ionIonisation.hh"
0076
0077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0078
0079 PhysicsList::PhysicsList() : G4VUserPhysicsList(), G4UImessenger()
0080
0081 {
0082   fpNumberOfSplit = new G4UIcmdWithAnInteger("/vrt/numberOfSplit", this);
0083   fpNumberOfSplit->SetGuidance("Number of split for particle splitting");
0084   fpNumberOfSplit->SetParameterName("NumberOfSplit", false);
0085   fpNumberOfSplit->SetDefaultValue(1);
0086
0087   defaultCutValue = 1 * micrometer;
0088   fCutForGamma = defaultCutValue;
0089   fCutForElectron = defaultCutValue;
0090   fCutForPositron = defaultCutValue;
0091   fCutForProton = defaultCutValue;
0092
0093   SetVerboseLevel(1);
0094
0095   //  fNumberOfSplit = 1;
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0099
0100 PhysicsList::~PhysicsList() {}
0101
0102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0103 void PhysicsList::SetNewValue(G4UIcommand* command, G4String newValue)
0104 {
0105   if (command == fpNumberOfSplit) fNumberOfSplit = fpNumberOfSplit->GetNewIntValue(newValue);
0106 }
0107
0108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0109
0110 void PhysicsList::ConstructParticle()
0111 {
0112   ConstructBosons();
0113   ConstructLeptons();
0114   ConstructBarions();
0115 }
0116
0117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0118
0119 void PhysicsList::ConstructBosons()
0120 {
0121   // gamma
0122   G4Gamma::GammaDefinition();
0123 }
0124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0125
0126 void PhysicsList::ConstructLeptons()
0127 {
0128   // leptons
0129   G4Electron::ElectronDefinition();
0130   G4Positron::PositronDefinition();
0131 }
0132
0133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0134
0135 // DNA
0136 #include "G4DNAGenericIonsManager.hh"
0137 // ENDDNA
0138
0139 void PhysicsList::ConstructBarions()
0140 {
0141   //  baryons
0142   G4Proton::ProtonDefinition();
0143   G4GenericIon::GenericIonDefinition();
0144
0145   // Geant4 DNA new particles
0146   G4DNAGenericIonsManager* genericIonsManager;
0147   genericIonsManager = G4DNAGenericIonsManager::Instance();
0148   genericIonsManager->GetIon("alpha++");
0149   genericIonsManager->GetIon("alpha+");
0150   genericIonsManager->GetIon("helium");
0151   genericIonsManager->GetIon("hydrogen");
0152 }
0153
0154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0155
0156 void PhysicsList::ConstructProcess()
0157 {
0158   AddTransportation();
0159   ConstructEM();
0160   ConstructGeneral();
0161 }
0162
0163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0164
0165 void PhysicsList::ConstructEM()
0166 {
0167   auto particleIterator = GetParticleIterator();
0168   particleIterator->reset();
0169   G4bool splitProcessActive = false;
0170   // Only if particle split is larger than 1, VRT is set
0171   if (fNumberOfSplit > 1) splitProcessActive = true;
0172
0173   while ((*particleIterator)()) {
0174     G4ParticleDefinition* particle = particleIterator->value();
0175     G4ProcessManager* pmanager = particle->GetProcessManager();
0176     G4String particleName = particle->GetParticleName();
0177
0178     // *********************************
0179     // 1) Processes for the World region
0180     // *********************************
0181
0182     if (particleName == "e-") {
0183       // STANDARD msc is active in the world
0184       G4eMultipleScattering* msc = new G4eMultipleScattering();
0185       msc->AddEmModel(1, new G4UrbanMscModel());
0186       pmanager->AddProcess(msc, -1, 1, -1);
0187
0188       // STANDARD ionisation is active in the world
0189       G4eIonisation* eion = new G4eIonisation();
0190       eion->SetEmModel(new G4MollerBhabhaModel());
0191       pmanager->AddProcess(eion, -1, 2, 2);
0192
0193       // DNA elastic is not active in the world
0194       G4DNAElastic* theDNAElasticProcess = new G4DNAElastic("e-_G4DNAElastic");
0195       theDNAElasticProcess->SetEmModel(new G4DummyModel());
0196       pmanager->AddDiscreteProcess(theDNAElasticProcess);
0197
0198       // DNA excitation is not active in the world
0199       G4DNAExcitation* dnaex = new G4DNAExcitation("e-_G4DNAExcitation");
0200       dnaex->SetEmModel(new G4DummyModel());
0201       pmanager->AddDiscreteProcess(dnaex);
0202
0203       // DNA ionisation is not active in the world
0204       G4DNAIonisation* dnaioni = new G4DNAIonisation("e-_G4DNAIonisation");
0205       dnaioni->SetEmModel(new G4DummyModel());
0206
0207       // Variance reduction is set here
0208       if (splitProcessActive) {
0209         G4String splitRegion = "Target";
0210         G4cout << "- Split for e-e created in e-_G4DNAIonisation process activated"
0211                << "with split number " << fNumberOfSplit << " in region " << splitRegion << "-- "
0212                << G4endl;
0213         SplitProcess* splitProcess = new SplitProcess(splitRegion, fNumberOfSplit);
0214         splitProcess->RegisterProcess(dnaioni);
0215
0216         pmanager->AddDiscreteProcess(splitProcess);
0217       }
0218       else {
0219         pmanager->AddDiscreteProcess(dnaioni);
0220       }
0221
0222       // DNA attachment is not active in the world
0223       G4DNAAttachment* dnaatt = new G4DNAAttachment("e-_G4DNAAttachment");
0224       dnaatt->SetEmModel(new G4DummyModel());
0225       pmanager->AddDiscreteProcess(dnaatt);
0226
0227       // DNA vib. excitation is not active in the world
0228       G4DNAVibExcitation* dnavib = new G4DNAVibExcitation("e-_G4DNAVibExcitation");
0229       dnavib->SetEmModel(new G4DummyModel());
0230       pmanager->AddDiscreteProcess(dnavib);
0231
0232       // THE FOLLOWING PROCESS WILL KILL ALL ELECTRONS BELOW A
0233       // SELECTED ENERY THRESHOLD
0234       // Capture of low-energy e-
0235       G4ElectronCapture* ecap = new G4ElectronCapture("Target", 7.4 * eV);
0236       pmanager->AddDiscreteProcess(ecap);
0237       // 7.4 eV is compatible with the validity range of the Champion's model
0238     }
0239     else if (particleName == "proton") {
0240       // STANDARD msc is active in the world
0241       G4hMultipleScattering* msc = new G4hMultipleScattering();
0242       msc->AddEmModel(1, new G4UrbanMscModel());
0243       pmanager->AddProcess(msc, -1, 1, -1);
0244
0245       // STANDARD ionisation is active in the world
0246       G4hIonisation* hion = new G4hIonisation();
0247       hion->SetEmModel(new G4BraggIonGasModel());
0248       hion->SetEmModel(new G4BetheBlochIonGasModel());
0249       pmanager->AddProcess(hion, -1, 2, 2);
0250
0251       // DNA excitation is not active in the world
0252       G4DNAExcitation* dnaex = new G4DNAExcitation("proton_G4DNAExcitation");
0253       dnaex->SetEmModel(new G4DummyModel());
0254       dnaex->SetEmModel(new G4DummyModel());
0255       pmanager->AddDiscreteProcess(dnaex);
0256
0257       // DNA ionisation is not active in the world
0258       G4DNAIonisation* dnaioni = new G4DNAIonisation("proton_G4DNAIonisation");
0259       dnaioni->SetEmModel(new G4DummyModel());
0260       dnaioni->SetEmModel(new G4DummyModel());
0261       pmanager->AddDiscreteProcess(dnaioni);
0262
0263       // DNA charge decrease is ACTIVE in the world since
0264       // no corresponding STANDARD process exist
0265       pmanager->AddDiscreteProcess(new G4DNAChargeDecrease("proton_G4DNAChargeDecrease"));
0266     }
0267     else if (particleName == "hydrogen") {
0268       // DNA processes are ACTIVE in the world since
0269       // no corresponding STANDARD processes exist
0270       pmanager->AddDiscreteProcess(new G4DNAIonisation("hydrogen_G4DNAIonisation"));
0271       pmanager->AddDiscreteProcess(new G4DNAExcitation("hydrogen_G4DNAExcitation"));
0272       pmanager->AddDiscreteProcess(new G4DNAChargeIncrease("hydrogen_G4DNAChargeIncrease"));
0273     }
0274     else if (particleName == "GenericIon") {
0275       // THIS IS NEEDED FOR STANDARD ALPHA G4ionIonisation PROCESS
0276
0277       // STANDARD msc is active in the world
0278       G4hMultipleScattering* msc = new G4hMultipleScattering();
0279       msc->AddEmModel(1, new G4UrbanMscModel());
0280       pmanager->AddProcess(msc, -1, 1, -1);
0281
0282       // STANDARD ionisation is active in the world
0283       G4ionIonisation* hion = new G4ionIonisation();
0284       hion->SetEmModel(new G4BraggIonGasModel());
0285       hion->SetEmModel(new G4BetheBlochIonGasModel());
0286       pmanager->AddProcess(hion, -1, 2, 2);
0287     }
0288     else if (particleName == "alpha") {
0289       // STANDARD msc is active in the world
0290       G4hMultipleScattering* msc = new G4hMultipleScattering();
0291       msc->AddEmModel(1, new G4UrbanMscModel());
0292       pmanager->AddProcess(msc, -1, 1, -1);
0293
0294       // STANDARD ionisation is active in the world
0295       G4ionIonisation* hion = new G4ionIonisation();
0296       hion->SetEmModel(new G4BraggIonGasModel());
0297       hion->SetEmModel(new G4BetheBlochIonGasModel());
0298       pmanager->AddProcess(hion, -1, 2, 2);
0299
0300       // DNA excitation is not active in the world
0301       G4DNAExcitation* dnaex = new G4DNAExcitation("alpha_G4DNAExcitation");
0302       dnaex->SetEmModel(new G4DummyModel());
0303       pmanager->AddDiscreteProcess(dnaex);
0304
0305       // DNA ionisation is not active in the world
0306       G4DNAIonisation* dnaioni = new G4DNAIonisation("alpha_G4DNAIonisation");
0307       dnaioni->SetEmModel(new G4DummyModel());
0308       pmanager->AddDiscreteProcess(dnaioni);
0309
0310       // DNA charge decrease is ACTIVE in the world since no
0311       // corresponding STANDARD process exist
0312       pmanager->AddDiscreteProcess(new G4DNAChargeDecrease("alpha_G4DNAChargeDecrease"));
0313     }
0314     else if (particleName == "alpha+") {
0315       // DNA processes are ACTIVE in the world since no
0316       // corresponding STANDARD processes exist
0317       pmanager->AddDiscreteProcess(new G4DNAExcitation("alpha+_G4DNAExcitation"));
0318       pmanager->AddDiscreteProcess(new G4DNAIonisation("alpha+_G4DNAIonisation"));
0319       pmanager->AddDiscreteProcess(new G4DNAChargeDecrease("alpha+_G4DNAChargeDecrease"));
0320       pmanager->AddDiscreteProcess(new G4DNAChargeIncrease("alpha+_G4DNAChargeIncrease"));
0321     }
0322     else if (particleName == "helium") {
0323       // DNA processes are ACTIVE in the world since no
0324       // corresponding STANDARD processes exist
0325       pmanager->AddDiscreteProcess(new G4DNAExcitation("helium_G4DNAExcitation"));
0326       pmanager->AddDiscreteProcess(new G4DNAIonisation("helium_G4DNAIonisation"));
0327       pmanager->AddDiscreteProcess(new G4DNAChargeIncrease("helium_G4DNAChargeIncrease"));
0328     }
0329   }
0330
0331   // **************************************
0332   // 2) Define processes for Target region
0333   // **************************************
0334
0335   // STANDARD EM processes should be inactivated when
0336   // corresponding DNA processes are used
0337   // - STANDARD EM e- processes are inactivated below 1 MeV
0338   // - STANDARD EM proton & alpha processes are inactivated below
0339   //   standEnergyLimit
0340   G4double standEnergyLimit = 9.9 * MeV;
0341   //
0342
0343   G4double massFactor = 1.0079 / 4.0026;
0344   G4EmConfigurator* em_config = G4LossTableManager::Instance()->EmConfigurator();
0345
0346   G4VEmModel* mod;
0347
0348   // *** e-
0349
0350   // ---> STANDARD EM processes are inactivated below 1 MeV
0351
0352   mod = new G4UrbanMscModel();
0353   mod->SetActivationLowEnergyLimit(1. * MeV);
0354   em_config->SetExtraEmModel("e-", "msc", mod, "Target");
0355
0356   mod = new G4MollerBhabhaModel();
0357   mod->SetActivationLowEnergyLimit(1 * MeV);
0358   em_config->SetExtraEmModel("e-", "eIoni", mod, "Target", 0.0, 100 * TeV,
0359                              new G4UniversalFluctuation());
0360
0361   // ---> DNA processes activated
0362
0363   mod = new G4DNAChampionElasticModel();
0364   em_config->SetExtraEmModel("e-", "e-_G4DNAElastic", mod, "Target", 7.4 * eV, 1. * MeV);
0365
0366   mod = new G4DNABornIonisationModel();
0367   em_config->SetExtraEmModel("e-", "e-_G4DNAIonisation", mod, "Target", 11. * eV, 1. * MeV);
0368   // Note: valid from 11 eV to 0.999.. MeV then switch to std models at
0369   // higher energies ; same for other models
0370
0371   mod = new G4DNABornExcitationModel();
0372   em_config->SetExtraEmModel("e-", "e-_G4DNAExcitation", mod, "Target", 9. * eV, 1. * MeV);
0373
0374   mod = new G4DNAMeltonAttachmentModel();
0375   em_config->SetExtraEmModel("e-", "e-_G4DNAAttachment", mod, "Target", 4. * eV, 13. * eV);
0376
0377   mod = new G4DNASancheExcitationModel();
0378   em_config->SetExtraEmModel("e-", "e-_G4DNAVibExcitation", mod, "Target", 2. * eV, 100. * eV);
0379
0380   // *** proton
0381
0382   // ---> STANDARD EM processes inactivated below standEnergyLimit
0383
0384   // STANDARD msc is still active
0385   // Inactivate following STANDARD processes
0386
0387   mod = new G4BraggIonGasModel();
0388   mod->SetActivationLowEnergyLimit(standEnergyLimit);
0389   em_config->SetExtraEmModel("proton", "hIoni", mod, "Target", 0.0, 2 * MeV,
0390                              new G4IonFluctuations());
0391
0392   mod = new G4BetheBlochIonGasModel();
0393   mod->SetActivationLowEnergyLimit(standEnergyLimit);
0394   em_config->SetExtraEmModel("proton", "hIoni", mod, "Target", 2 * MeV, 100 * TeV,
0395                              new G4UniversalFluctuation());
0396
0397   // ---> DNA processes activated
0398
0399   mod = new G4DNARuddIonisationModel();
0400   em_config->SetExtraEmModel("proton", "proton_G4DNAIonisation", mod, "Target", 0.0, 0.5 * MeV);
0401
0402   mod = new G4DNABornIonisationModel();
0403   em_config->SetExtraEmModel("proton", "proton_G4DNAIonisation", mod, "Target", 0.5 * MeV,
0404                              10 * MeV);
0405
0406   mod = new G4DNAMillerGreenExcitationModel();
0407   em_config->SetExtraEmModel("proton", "proton_G4DNAExcitation", mod, "Target", 10 * eV, 0.5 * MeV);
0408
0409   mod = new G4DNABornExcitationModel();
0410   em_config->SetExtraEmModel("proton", "proton_G4DNAExcitation", mod, "Target", 0.5 * MeV,
0411                              10 * MeV);
0412
0413   // *** alpha
0414
0415   // ---> STANDARD EM processes inactivated below standEnergyLimit
0416
0417   // STANDARD msc is still active
0418   // Inactivate following STANDARD processes
0419
0420   mod = new G4BraggIonGasModel();
0421   mod->SetActivationLowEnergyLimit(standEnergyLimit);
0422   em_config->SetExtraEmModel("alpha", "ionIoni", mod, "Target", 0.0, 2 * MeV / massFactor,
0423                              new G4IonFluctuations());
0424
0425   mod = new G4BetheBlochIonGasModel();
0426   mod->SetActivationLowEnergyLimit(standEnergyLimit);
0427   em_config->SetExtraEmModel("alpha", "ionIoni", mod, "Target", 2 * MeV / massFactor, 100 * TeV,
0428                              new G4UniversalFluctuation());
0429
0430   // ---> DNA processes activated
0431
0432   mod = new G4DNARuddIonisationModel();
0433   em_config->SetExtraEmModel("alpha", "alpha_G4DNAIonisation", mod, "Target", 0.0, 10 * MeV);
0434
0435   mod = new G4DNAMillerGreenExcitationModel();
0436   em_config->SetExtraEmModel("alpha", "alpha_G4DNAExcitation", mod, "Target", 1 * keV, 10 * MeV);
0437 }
0438
0439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0440
0441 void PhysicsList::ConstructGeneral() {}
0442
0443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0444
0445 void PhysicsList::SetCuts()
0446 {
0447   if (verboseLevel > 0) {
0448     G4cout << "PhysicsList::SetCuts:";
0449     G4cout << "CutLength : " << G4BestUnit(defaultCutValue, "Length") << G4endl;
0450   }
0451
0452   // set cut values for gamma at first and for e- second and next for e+,
0453   // because some processes for e+/e- need cut values for gamma
0454   SetCutValue(fCutForGamma, "gamma");
0455   SetCutValue(fCutForElectron, "e-");
0456   SetCutValue(fCutForPositron, "e+");
0457   SetCutValue(fCutForProton, "proton");
0458
0459   if (verboseLevel > 0) {
0460     DumpCutValuesTable();
0461   }
0462 }
0463
0464 void PhysicsList::SetNumberOfSplit(G4int nb)
0465 {
0466   fNumberOfSplit = nb;
0467 }
0468
0469 G4int PhysicsList::GetNumberOfSplit()
0470 {
0471   return fNumberOfSplit;
0472 }