Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!
Geant4 Version: 11
Operating System: Window 10
Compiler/Version:
CMake Version:
Dear all,
I trying to construct a UserPhysicsList for proton interaction. Here is my PhysicsList class:
#include "PhysicsList.hh"
#include "G4ProcessManager.hh"
#include "G4BosonConstructor.hh"
#include "G4LeptonConstructor.hh"
#include "G4MesonConstructor.hh"
#include "G4BaryonConstructor.hh"
#include "G4IonConstructor.hh"
#include "G4ShortLivedConstructor.hh"
#include "G4Gamma.hh"
#include "G4Proton.hh"
#include "G4PhysicsListHelper.hh"
#include "G4SystemOfUnits.hh"
PhysicsList::PhysicsList()
: G4VUserPhysicsList()
{
}
PhysicsList::~PhysicsList()
{
}
void PhysicsList::ConstructParticle()
{
G4BosonConstructor bosonConst;
bosonConst.ConstructParticle();
G4LeptonConstructor leptonConst;
leptonConst.ConstructParticle();
G4MesonConstructor mesonConst;
mesonConst.ConstructParticle();
G4BaryonConstructor baryonConst;
baryonConst.ConstructParticle();
G4IonConstructor ionConst;
ionConst.ConstructParticle();
G4ShortLivedConstructor shortLivedConst;
shortLivedConst.ConstructParticle();
}
// gamma
#include "G4PhotoElectricEffect.hh"
#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4RayleighScattering.hh"
#include "G4LivermorePhotoElectricModel.hh" // phot models
#include "G4PenelopePhotoElectricModel.hh"
#include "G4KleinNishinaModel.hh" // compt models
#include "G4BetheHeitlerModel.hh" // conv models
#include "G4PenelopeGammaConversionModel.hh"
#include "G4LivermoreRayleighModel.hh" //rayl models
// e- and e+
#include "G4eMultipleScattering.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
// proton
#include "G4hMultipleScattering.hh"
#include "G4hIonisation.hh"
// hadron
#include "G4HadronInelasticProcess.hh"
#include "G4hMultipleScattering.hh"
#include "G4hIonisation.hh"
#include "G4ionIonisation.hh"
#include "G4NuclearStopping.hh"
// process
#include "G4LowEGammaNuclearModel.hh"
#include "G4HadronElasticProcess.hh"
#include "G4HadronInelasticProcess.hh"
// #include "G4AlphaInelasticProcess.hh"
// data
#include "G4VCrossSectionDataSet.hh"
#include "G4PhotoNuclearCrossSection.hh"
#include "G4BGGNucleonElasticXS.hh"
#include "G4ChipsProtonElasticXS.hh"
#include "G4ChipsProtonInelasticXS.hh"
// #include "G4TripathiLightCrossSection.hh"
// model
#include "G4CascadeInterface.hh"
#include "G4IonParametrisedLossModel.hh"
#include "G4ChipsElasticModel.hh"
#include "G4BinaryLightIonReaction.hh"
#include "G4BinaryLightIonReaction.hh"
// #include "G4LEAlphaInelastic.hh"
#include "G4HadronicInteraction.hh"
void PhysicsList::ConstructProcess()
{
G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
auto particleIterator = GetParticleIterator();
particleIterator->reset();
while((*particleIterator)())
{
// get the current particle definition
G4ParticleDefinition* particle = particleIterator->value();
G4String particleName = particle->GetParticleName();
G4Decay* aDecay = new G4Decay();
if (aDecay->IsApplicable(*particle)) ph->RegisterProcess(aDecay, particle);
if (particleName == "gamma")
{
// Em proceses
// 1. photoelectric
G4PhotoElectricEffect* photProcess = new G4PhotoElectricEffect();
photProcess->SetEmModel(new G4LivermorePhotoElectricModel());
ph->RegisterProcess(photProcess, particle);
// 2. compton scattering
G4ComptonScattering* comptProcess = new G4ComptonScattering();
comptProcess->SetEmModel(new G4KleinNishinaModel());
ph->RegisterProcess(comptProcess, particle);
// 3. gamma conversion
G4GammaConversion* convProcess = new G4GammaConversion();
convProcess->SetEmModel(new G4BetheHeitlerModel());
ph->RegisterProcess(convProcess, particle);
// 4. rayleigh scattering
G4RayleighScattering* raylProcess = new G4RayleighScattering();
raylProcess->SetEmModel(new G4LivermoreRayleighModel());
ph->RegisterProcess(raylProcess, particle);
// Hadron proceses
// 5. gamma-nuclear
G4HadronInelasticProcess* photoNuclProcess = new G4HadronInelasticProcess("photonNuclear", G4Gamma::Definition());;
photoNuclProcess->AddDataSet(new G4GammaNuclearXS()); // data
// low gamma-nuclear model from 0 - 200 MeV
G4LowEGammaNuclearModel* theLowGammaNuclear = new G4LowEGammaNuclearModel();
theLowGammaNuclear->SetMaxEnergy(200*MeV);
photoNuclProcess->RegisterMe(theLowGammaNuclear);
G4CascadeInterface* theCascade = new G4CascadeInterface();
theCascade->SetMinEnergy(199*MeV);
theCascade->SetMaxEnergy(10*GeV);
photoNuclProcess->RegisterMe(theCascade);
G4ProcessManager* gammaProcessManager = G4Gamma::Gamma()->GetProcessManager();
gammaProcessManager->AddDiscreteProcess(photoNuclProcess);
}
else if (particleName == "e-")
{
// EM processes
ph->RegisterProcess(new G4eMultipleScattering, particle);
ph->RegisterProcess(new G4eIonisation, particle);
ph->RegisterProcess(new G4eBremsstrahlung, particle);
}
else if (particleName == "e+")
{
ph->RegisterProcess(new G4eMultipleScattering, particle);
ph->RegisterProcess(new G4eIonisation, particle);
ph->RegisterProcess(new G4eBremsstrahlung, particle);
ph->RegisterProcess(new G4eplusAnnihilation, particle);
}
else if (particleName == "proton")
{
// EM processes
ph->RegisterProcess(new G4hMultipleScattering, particle);
ph->RegisterProcess(new G4hIonisation, particle);
// 1. Hadron elastic process
G4HadronElasticProcess * theProtonElasticProcess = new G4HadronElasticProcess("hadElastic");
theProtonElasticProcess->AddDataSet(new G4BGGNucleonElasticXS(G4Proton::Definition()));
theProtonElasticProcess->RegisterMe(new G4ChipsElasticModel());
// 2. Hadron inelastic process
G4HadronInelasticProcess* theProtonInelasticProcess = new G4HadronInelasticProcess("protonInelastic", G4Proton::Definition());
theProtonInelasticProcess->AddDataSet(new G4BGGNucleonInelasticXS (G4Proton::Definition()));
//
// models
G4BinaryCascade* theBinaryCascadeModel = new G4BinaryCascade();
theBinaryCascadeModel->SetMinEnergy(0.*eV);
theBinaryCascadeModel->SetMaxEnergy(1.5*GeV);
//theProtonInelasticProcess->RegisterMe(theBinaryCascadeModel);
G4CascadeInterface* theBertiniCascadeModel = new G4CascadeInterface();
theBertiniCascadeModel->SetMinEnergy(0.*GeV);
theBertiniCascadeModel->SetMaxEnergy(6.*GeV);
theProtonInelasticProcess->RegisterMe(theBertiniCascadeModel);
G4ProcessManager* protonProcessManager = G4Proton::Proton()->GetProcessManager();
protonProcessManager->AddDiscreteProcess(theProtonElasticProcess);
protonProcessManager->AddDiscreteProcess(theProtonInelasticProcess);
}
else if(particleName == "anti_proton")
{
ph->RegisterProcess(new G4hMultipleScattering, particle);
ph->RegisterProcess(new G4hIonisation, particle);
}
else if (particleName == "alpha" || particleName == "He3")
{
ph->RegisterProcess(new G4hMultipleScattering, particle);
ph->RegisterProcess(new G4ionIonisation, particle);
ph->RegisterProcess(new G4NuclearStopping, particle);
// Hadron process
// 1. elastic
G4HadronElasticProcess* alphaElasticProcess = new G4HadronElasticProcess("AlphaElastic");
// cross section
G4VCrossSectionDataSet* theGlauberGribovCrSc = new G4VCrossSectionDataSet("G4GlauberGribovCrossSection");
alphaElasticProcess->AddDataSet(theGlauberGribovCrSc);
// model
G4HadronicInteraction* theChipsElasticModel = new G4ChipsElasticModel();
alphaElasticProcess->RegisterMe(theChipsElasticModel);
G4Alpha::Alpha()->GetProcessManager()->AddDiscreteProcess(alphaElasticProcess);
//G4He3::He3()->GetProcessManager()->AddDiscreteProcess(alphaElasticProcess);
// 2. inelastic
G4HadronInelasticProcess* theAlphaInelasticProcess = new G4HadronInelasticProcess("AlphaInelastic");
//
// cross section dataset
G4VCrossSectionDataSet* theAlphaInelasticDataset = new G4VCrossSectionDataSet("GheishaInelastic");
theAlphaInelasticProcess->AddDataSet(theAlphaInelasticDataset);
//
// model
G4HadronicInteraction* theAlphaInelasticModel = new G4HadronicInteraction("G4LEAlphaInelastic");
theAlphaInelasticProcess->RegisterMe(theAlphaInelasticModel);
ph->RegisterProcess(theAlphaInelasticProcess, particle);
}
else if (particleName == "GenericIon")
{
// EM proceses
//
ph->RegisterProcess(new G4hMultipleScattering, particle);
//
G4ionIonisation* ionIoni = new G4ionIonisation;
ionIoni->SetEmModel(new G4IonParametrisedLossModel());
ph->RegisterProcess(ionIoni, particle);
//
ph->RegisterProcess(new G4NuclearStopping(), particle);
}
else if((!particle->IsShortLived()) &&
(particle->GetPDGCharge() != 0.0) &&
(particle->GetParticleName() != "chargedgeantino"))
{
// all others charged particles except geantino
ph->RegisterProcess(new G4hMultipleScattering(), particle);
ph->RegisterProcess(new G4hIonisation(), particle);
ph->RegisterProcess(new G4NuclearStopping(), particle);
}
}
}
and PhysicsList.hh
#ifndef PhysicsList_h
#define PhysicsList_h 1
#include "G4VUserPhysicsList.hh"
#include "globals.hh"
//particle
#include "G4Gamma.hh"
#include "G4Proton.hh"
#include "G4AntiProton.hh"
#include "G4Alpha.hh"
#include "G4He3.hh"
// ****************************************************************************
// EM process
// ****************************************************************************
// --- gamma ---
#include "G4PhotoElectricEffect.hh"
#include "G4LivermorePhotoElectricModel.hh" // model
#include "G4ComptonScattering.hh"
#include "G4KleinNishinaModel.hh" // model
#include "G4GammaConversion.hh"
#include "G4BetheHeitlerModel.hh" // model
#include "G4RayleighScattering.hh"
#include "G4LivermoreRayleighModel.hh" // model
// --- e- & e+ ---
#include "G4eMultipleScattering.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
// --- proton ---
#include "G4hMultipleScattering.hh"
#include "G4hIonisation.hh"
// --- alpha & He3 ---
#include "G4hMultipleScattering.hh"
#include "G4ionIonisation.hh"
#include "G4NuclearStopping.hh"
/// --- ions ---
#include "G4hMultipleScattering.hh"
// ****************************************************************************
// Hadronic process
// ****************************************************************************
// --- gamma ---
#include "G4HadronInelasticProcess.hh" // "photonNuclear"
#include "G4GammaNuclearXS.hh"
#include "G4PhotoNuclearCrossSection.hh" // datasets
#include "G4LowEGammaNuclearModel.hh" // model
#include "G4CascadeInterface.hh" // model
// --- proton ---
#include "G4HadronElasticProcess.hh" // "hadElastic"
#include "G4BGGNucleonElasticXS.hh" // datasets
#include "G4ChipsElasticModel.hh" // model
// Decay process
#include "G4Decay.hh"
#include "G4HadronicInteraction.hh"
//data
#include "G4VCrossSectionDataSet.hh"
#include "G4BGGNucleonInelasticXS.hh"
//model
#include "G4ChipsElasticModel.hh"
#include "G4BinaryCascade.hh"
#include "G4CascadeInterface.hh"
class G4VPhysicsConstructor;
class PhysicsList: public G4VUserPhysicsList
{
public:
PhysicsList();
virtual ~PhysicsList();
virtual void ConstructParticle();
virtual void ConstructProcess();
};
#endif
- When I run the application, my detector can not record protons.
- In the following figure, electrons (blue) are generated instead of protons. I don’t know what the is reason although my source still emmited protons.
- If I use Reference PhysicsList, I saw protons (red in below image) emitted from the source.
G4String physName = "FTFP_BERT";
G4PhysListFactory factory;
G4VModularPhysicsList* physicsList = factory.GetReferencePhysList(physName);
Can anyone help me to fix problem? Thank you in advance.