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