//-----------------------------------------------------------------------------
// File and Version Information:
// $Id: BgsPhysicsList.cc,v 1.76 2005/01/14 02:04:40 lockman Exp $
//
// Description:
// BgsPhysicsList
//
// Environment:
// Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author List:
// Torre Wenaus, David Williams, Marc Verderi, Bill Lockman, Dennis Wright
//
// Copyright Information:
// Copyright (C) 2001 SLAC
//
// Created:
// Modification history:
//
//-----------------------------------------------------------------------------
#include "BaBar/BaBar.hh"
#include
#include "Bogus/BgsPhysicsList.hh"
#include
#include "Bogus/BgsPDGEncode.hh"
#include "Bogus/BgsControl.hh"
#include "Bogus/BgsLooperDeath.hh"
#include "Bogus/BgsChargedStepDeath.hh"
#include "Bogus/BgsPhysicsRegistrar.hh"
#include "Bogus/BgsGenocide.hh"
#include "Bogus/BgsGentleGenocide.hh"
#include "Bogus/BgsLimitedTransportation.hh"
#include "ErrLogger/ErrLog.hh"
#include "geant4/G4ParticleDefinition.hh"
#include "geant4/G4ParticleWithCuts.hh"
#include "geant4/G4ProcessManager.hh"
#include "geant4/G4ProcessVector.hh"
#include "geant4/G4ParticleTypes.hh"
#include "geant4/G4ParticleTable.hh"
#include "geant4/G4FastSimulationManagerProcess.hh"
#include "geant4/G4ComptonScattering.hh"
#include "geant4/G4GammaConversion.hh"
#include "geant4/G4PhotoElectricEffect.hh"
#include "geant4/G4MultipleScattering.hh"
#include "geant4/G4eIonisation.hh"
#include "geant4/G4eBremsstrahlung.hh"
#include "geant4/G4eplusAnnihilation.hh"
#include "geant4/G4MuIonisation.hh"
#include "geant4/G4MuBremsstrahlung.hh"
#include "geant4/G4MuPairProduction.hh"
#include "geant4/G4hIonisation.hh"
#include "geant4/G4Decay.hh"
#include "geant4/G4Material.hh"
#include "geant4/G4ProductionCuts.hh"
#include "geant4/G4ProductionCutsTable.hh"
#include "geant4/G4MaterialCutsCouple.hh"
//
// Gamma- and electro-nuclear models and processes
//
#include "geant4/G4GammaNuclearReaction.hh"
#include "geant4/G4ElectroNuclearReaction.hh"
#include "geant4/G4TheoFSGenerator.hh"
#include "geant4/G4GeneratorPrecompoundInterface.hh"
#include "geant4/G4QGSModel.hh"
#include "geant4/G4GammaParticipants.hh"
#include "geant4/G4QGSMFragmentation.hh"
#include "geant4/G4ExcitedStringDecay.hh"
#include "geant4/G4PhotoNuclearProcess.hh"
#include "geant4/G4ElectronNuclearProcess.hh"
#include "geant4/G4PositronNuclearProcess.hh"
//
// Hadronic processes
//
#include "geant4/G4HadronElasticProcess.hh"
#include "geant4/G4HadronCaptureProcess.hh"
#include "geant4/G4HadronFissionProcess.hh"
#include "geant4/G4PionPlusInelasticProcess.hh"
#include "geant4/G4PionMinusInelasticProcess.hh"
#include "geant4/G4KaonPlusInelasticProcess.hh"
#include "geant4/G4KaonMinusInelasticProcess.hh"
#include "geant4/G4KaonZeroLInelasticProcess.hh"
#include "geant4/G4KaonZeroSInelasticProcess.hh"
#include "geant4/G4ProtonInelasticProcess.hh"
#include "geant4/G4AntiProtonInelasticProcess.hh"
#include "geant4/G4NeutronInelasticProcess.hh"
#include "geant4/G4AntiNeutronInelasticProcess.hh"
#include "geant4/G4LambdaInelasticProcess.hh"
#include "geant4/G4AntiLambdaInelasticProcess.hh"
#include "geant4/G4DeuteronInelasticProcess.hh"
#include "geant4/G4TritonInelasticProcess.hh"
#include "geant4/G4AlphaInelasticProcess.hh"
#include "geant4/G4ShortLivedConstructor.hh"
//
// Hadronic interaction models
// Low energy (E < 20 GeV) part only
//
#include "geant4/G4LElastic.hh"
// #include "geant4/G4LEPionPlusInelastic.hh"
// #include "geant4/G4LEPionMinusInelastic.hh"
#include "geant4/G4LEKaonPlusInelastic.hh"
#include "geant4/G4LEKaonMinusInelastic.hh"
#include "geant4/G4LEKaonZeroLInelastic.hh"
#include "geant4/G4LEKaonZeroSInelastic.hh"
#include "geant4/G4LEProtonInelastic.hh"
#include "geant4/G4LEAntiProtonInelastic.hh"
#include "geant4/G4LENeutronInelastic.hh"
#include "geant4/G4LEAntiNeutronInelastic.hh"
#include "geant4/G4LELambdaInelastic.hh"
#include "geant4/G4LEAntiLambdaInelastic.hh"
#include "geant4/G4LEDeuteronInelastic.hh"
#include "geant4/G4LETritonInelastic.hh"
#include "geant4/G4LEAlphaInelastic.hh"
#include "geant4/G4NeutronHPElastic.hh"
#include "geant4/G4NeutronHPElasticData.hh"
#include "geant4/G4NeutronHPInelastic.hh"
#include "geant4/G4NeutronHPInelasticData.hh"
#include "geant4/G4NeutronHPFission.hh"
#include "geant4/G4NeutronHPFissionData.hh"
#include "geant4/G4NeutronHPCapture.hh"
#include "geant4/G4NeutronHPCaptureData.hh"
#include "geant4/G4LFission.hh"
#include "geant4/G4LCapture.hh"
#include "geant4/G4CascadeInterface.hh"
//
// Pi+/pi- inelastic cross sections
//
#include "geant4/G4PiNuclearCrossSection.hh"
using std::ostream;
BgsPhysicsList::BgsPhysicsList( BgsControl *theControl,
BgsPhysicsRegistrar *pr)
: G4VUserPhysicsList(),
control(theControl),
physicsRegistrar(pr),
theLooperDeath(0) // looperdeath process
{
SetVerboseLevel(2);
bertini_model = new G4CascadeInterface;
}
BgsPhysicsList::~BgsPhysicsList()
{
delete theLooperDeath;
}
void
BgsPhysicsList::ConstructParticle()
{
// In this method, static member functions should be called
// for all particles which you want to use.
// This ensures that objects of these particle types will be
// created in the program.
ConstructBosons();
ConstructLeptons();
ConstructMesons();
ConstructBaryons();
ConstructIons();
//
// Check to make sure our particles have valid PDG encodings
//
BgsPDGEncode encode(false);
G4ParticleTable *table = G4ParticleTable::GetParticleTable();
G4int nProb = 0;
G4int n = table->entries();
while(n--) {
G4ParticleDefinition *part = table->GetParticle(n);
if (encode.pdt(part) == 0) {
nProb++;
ErrMsg(error) << "The Geant4 particle \""
<< part->GetParticleName()
<< "\" is not recognized by BgsPDGEncode" << endmsg;
}
}
if (nProb > 0) ErrMsg(fatal) << "One or more PDG encoding errors" << endmsg;
// Add short lived particles for high energy models,
// but don't check PDG codes - they are not propagated in Bogus anyway
G4ShortLivedConstructor shortLived;
shortLived.ConstructParticle();
}
void
BgsPhysicsList::ConstructBosons()
{
// pseudo-particles
G4Geantino::GeantinoDefinition();
G4ChargedGeantino::ChargedGeantinoDefinition();
// gamma
G4Gamma::GammaDefinition();
// optical photon
G4OpticalPhoton::OpticalPhotonDefinition();
}
void
BgsPhysicsList::ConstructLeptons()
{
// leptons
G4Electron::ElectronDefinition();
G4Positron::PositronDefinition();
G4MuonPlus::MuonPlusDefinition();
G4MuonMinus::MuonMinusDefinition();
G4TauPlus::TauPlusDefinition();
G4TauMinus::TauMinusDefinition();
G4NeutrinoE::NeutrinoEDefinition();
G4AntiNeutrinoE::AntiNeutrinoEDefinition();
G4NeutrinoMu::NeutrinoMuDefinition();
G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
G4NeutrinoTau::NeutrinoTauDefinition();
G4AntiNeutrinoTau::AntiNeutrinoTauDefinition();
}
void
BgsPhysicsList::ConstructMesons()
{
// mesons
G4PionPlus ::PionPlusDefinition();
G4PionMinus ::PionMinusDefinition();
G4PionZero ::PionZeroDefinition();
G4Eta ::EtaDefinition();
G4EtaPrime ::EtaPrimeDefinition();
// G4RhoZero ::RhoZeroDefinition();
G4KaonPlus ::KaonPlusDefinition();
G4KaonMinus ::KaonMinusDefinition();
G4KaonZero ::KaonZeroDefinition();
G4AntiKaonZero ::AntiKaonZeroDefinition();
G4KaonZeroLong ::KaonZeroLongDefinition();
G4KaonZeroShort::KaonZeroShortDefinition();
}
void
BgsPhysicsList::ConstructBaryons()
{
// baryons
G4Proton ::ProtonDefinition();
G4AntiProton ::AntiProtonDefinition();
G4Neutron ::NeutronDefinition();
G4AntiNeutron ::AntiNeutronDefinition();
G4Lambda ::LambdaDefinition();
G4AntiLambda ::AntiLambdaDefinition();
G4SigmaPlus ::SigmaPlusDefinition();
G4SigmaZero ::SigmaZeroDefinition();
G4SigmaMinus ::SigmaMinusDefinition();
G4AntiSigmaPlus ::AntiSigmaPlusDefinition();
G4AntiSigmaZero ::AntiSigmaZeroDefinition();
G4AntiSigmaMinus::AntiSigmaMinusDefinition();
G4XiZero ::XiZeroDefinition();
G4XiMinus ::XiMinusDefinition();
G4AntiXiZero ::AntiXiZeroDefinition();
G4AntiXiMinus ::AntiXiMinusDefinition();
G4OmegaMinus ::OmegaMinusDefinition();
G4AntiOmegaMinus::AntiOmegaMinusDefinition();
G4Deuteron ::DeuteronDefinition();
G4Triton ::TritonDefinition();
G4Alpha ::AlphaDefinition();
}
void
BgsPhysicsList::ConstructIons()
{
G4GenericIon::GenericIonDefinition();
}
void
BgsPhysicsList::ConstructProcess()
{
if (control->UseBgsTran()) {
AddBgsTransportation(control->GetMaxTrackingStepSize());
ErrMsg(trace) << " +^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^+ " <<
endmsg;
ErrMsg(trace) << " +v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v+ " <<
endmsg;
ErrMsg(trace) << " +--- ---+ " <<
endmsg;
ErrMsg(trace) << " +--- BgsTransportation ---+ " <<
endmsg;
ErrMsg(trace) << " +--- USED !! ---+ " <<
endmsg;
ErrMsg(trace) << " +--- ---+ " <<
endmsg;
ErrMsg(trace) << " +^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^+ " <<
endmsg;
ErrMsg(trace) << " +v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v+ " <<
endmsg;
}
else {
AddTransportation();
ErrMsg(trace) << " +^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^+ " <<
endmsg;
ErrMsg(trace) << " +v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v+ " <<
endmsg;
ErrMsg(trace) << " +--- ---+ " <<
endmsg;
ErrMsg(trace) << " +--- G4Transportation ---+ " <<
endmsg;
ErrMsg(trace) << " +--- USED !! ---+ " <<
endmsg;
ErrMsg(trace) << " +--- ---+ " <<
endmsg;
ErrMsg(trace) << " +^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^+ " <<
endmsg;
ErrMsg(trace) << " +v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v^v+ " <<
endmsg;
}
AddParameterisation();
ConstructEM();
ConstructLeptHad();
ConstructHad();
ConstructGeneral();
ConstructNeutrinoGenocide();
ConstructIonFix();
ConstructNeutronFix();
}
void
BgsPhysicsList::AddBgsTransportation( G4double maxStepSize )
{
static const bool beParanoid = false;
BgsTransportation *theTransportationProcess=
new BgsLimitedTransportation(maxStepSize*mm, beParanoid,
control->GetEnableStepLimitVolume());
// loop over all particles in G4ParticleTable
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
if (!particle->IsShortLived()) {
// Add transportation process for all particles other than "shortlived"
if ( pmanager == 0) {
// Error !! no process manager
ErrMsg(fatal) << "G4VUserPhysicsList::AddTransportation : no process manager" << endmsg;
} else {
// add transportation with ordering = ( -1, "first", "first" )
pmanager ->AddProcess(theTransportationProcess);
pmanager ->SetProcessOrderingToFirst(theTransportationProcess,
idxAlongStep);
pmanager ->SetProcessOrderingToFirst(theTransportationProcess,
idxPostStep);
physicsRegistrar->Register( theTransportationProcess, pmanager,
"transport" );
}
} else {
// shortlived particle case
}
}
}
// **************************** EM Processes ********************************
void
BgsPhysicsList::ConstructEM()
{
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if (particleName == "gamma") {
// gamma
// Construct processes for gamma
AddDiscreteProcess( new G4PhotoElectricEffect(), pmanager, "phot",
GVertex::photoElectricEffect );
AddDiscreteProcess( new G4ComptonScattering(), pmanager, "comp",
GVertex::comptonScattering );
AddDiscreteProcess( new G4GammaConversion(), pmanager, "conv",
GVertex::gammaConversion );
} else if (particleName == "e-") {
//electron
// Construct processes for electron
G4MultipleScattering *ms = new G4MultipleScattering();
ms->SetLateralDisplasmentFlag(false);
AddProcess( ms, -1, 1,1, pmanager, "mult",
GVertex::trackingError );
AddProcess( new G4eIonisation(), -1, 2,2, pmanager, "ioni",
GVertex::ionization );
AddProcess( new G4eBremsstrahlung(), -1,-1,3, pmanager, "brem",
GVertex::bremsstrahlung );
} else if (particleName == "e+") {
//positron
// Construct processes for positron
G4MultipleScattering *ms = new G4MultipleScattering();
ms->SetLateralDisplasmentFlag(false);
AddProcess( ms, -1, 1,1, pmanager, "mult",
GVertex::trackingError );
AddProcess( new G4eIonisation(), -1, 2,2, pmanager, "ioni",
GVertex::ionization );
AddProcess( new G4eBremsstrahlung(), -1,-1,3, pmanager, "brem",
GVertex::bremsstrahlung );
AddProcess( new G4eplusAnnihilation(), 0,-1,4, pmanager, "anni",
GVertex::eplusAnnihilation );
} else if( particleName == "mu+" ||
particleName == "mu-" ) {
//muon
// Construct processes for muon+
G4MultipleScattering *ms = new G4MultipleScattering();
ms->SetLateralDisplasmentFlag(false);
AddProcess( ms, -1, 1,1, pmanager, "mult",
GVertex::trackingError );
AddProcess( new G4MuIonisation(), -1, 2,2, pmanager, "ioni",
GVertex::ionization );
AddProcess( new G4MuBremsstrahlung(), -1,-1,3, pmanager, "mbre",
GVertex::bremsstrahlung );
AddProcess( new G4MuPairProduction(), -1,-1,4, pmanager, "pair",
GVertex::muonPairProduction );
} else if ((!particle->IsShortLived()) &&
(particle->GetPDGCharge() != 0.0) &&
(particle->GetParticleName() != "chargedgeantino")) {
// all others charged particles except geantino
G4int AP=1;
if (particle->GetParticleName() == "GenericIon") {
ostream& o = ErrMsg(warning);
o << "*********************************************************************" << endmsg;
o << "*** Disabling G4MultipleScattering process for particle " << particle->GetParticleName() << endmsg;
o << "*** This process should be re-enabled in Geant4 7.0" << endmsg;
o << "*********************************************************************" << endmsg;
} else {
G4MultipleScattering *ms = new G4MultipleScattering();
ms->SetLateralDisplasmentFlag(false);
AddProcess( ms, -1,AP,AP, pmanager, "mult", GVertex::trackingError );
AP++;
}
AddProcess( new G4hIonisation(), -1,AP,AP, pmanager, "ioni",
GVertex::ionization );
}
}
}
void BgsPhysicsList::AddProcess( G4VProcess *process,
G4int ordAtRestDoIt,
G4int ordAlongStepDoIt,
G4int ordPostStepDoIt,
G4ProcessManager *manager,
const char *category,
GVertex::Cause cause )
{
manager->AddProcess( process, ordAtRestDoIt, ordAlongStepDoIt,
ordPostStepDoIt );
physicsRegistrar->Register( process, manager, category, cause );
}
void BgsPhysicsList::AddDiscreteProcess( G4VProcess *process,
G4ProcessManager *manager,
const char *category,
GVertex::Cause cause )
{
manager->AddDiscreteProcess( process );
physicsRegistrar->Register( process, manager, category, cause );
}
void BgsPhysicsList::AddElasticHadronProcess(G4HadronElasticProcess *process,
G4LElastic *model,
G4ProcessManager *manager,
const char *category,
GVertex::Cause cause )
{
process->RegisterMe(model); // Register interaction model with process
manager->AddDiscreteProcess(process);
physicsRegistrar->Register(process,manager,category,cause);
}
void
BgsPhysicsList::SetStatusEmProcess()
{
}
// ************************** Parameterisation ******************************
void
BgsPhysicsList::AddParameterisation()
{
G4FastSimulationManagerProcess* theFastSimulationManagerProcess =
new G4FastSimulationManagerProcess();
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
pmanager->AddDiscreteProcess( theFastSimulationManagerProcess );
}
}
// ************************** Generic Processes ******************************
void
BgsPhysicsList::ConstructGeneral()
{
// Add Decay Process
G4Decay* theDecayProcess = new G4Decay();
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
if (theDecayProcess->IsApplicable(*particle)) {
pmanager ->AddProcess( theDecayProcess );
// set ordering for PostStepDoIt and AtRestDoIt
pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
physicsRegistrar->
Register( theDecayProcess, pmanager, "dcay", GVertex::decay );
}
}
if (control->GetLooperCut()>0) {
// Set special process to kill loopers
theLooperDeath = new BgsLooperDeath( control->GetLooperCut()*MeV );
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
if (theLooperDeath->IsApplicable(*particle)) {
G4ProcessManager* pmanager = particle->GetProcessManager();
pmanager->AddProcess(theLooperDeath, -1, -1, 5);
physicsRegistrar->
Register( theLooperDeath, pmanager, "loop", GVertex::looperDeath );
}
}
ErrMsg(warning) << "Loopers with pt < " << control->GetLooperCut()
<< " MeV will be killed" << endmsg;
}
if (control->GetMaxNumberOfSteps()>0) {
//
// Set special process to kill runaway particles
// Only needed if dE/dx is turned off!
//
// Do not abuse!
//
theStepDeath = new BgsChargedStepDeath( control->GetMaxNumberOfSteps() );
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
if (theStepDeath->IsApplicable(*particle)) {
G4ProcessManager* pmanager = particle->GetProcessManager();
pmanager->AddProcess(theStepDeath, -1, -1, 5);
physicsRegistrar->
Register( theStepDeath, pmanager, "maxStep", GVertex::runAway );
}
}
ErrMsg(warning)
<< "\n Charged particles will be killed if they take more than "
<< control->GetMaxNumberOfSteps() << " steps.\n"
<< " If you do not understand this message, you should be very concerned.\n"
<< " If this message appears in production, you should be very upset." << endmsg;
}
}
// ************************** Hadronic Processes ******************************
void
BgsPhysicsList::ConstructLeptHad()
{
//
// Gamma-nuclear process
//
// low energy part
G4GammaNuclearReaction* lowEGammaModel = new G4GammaNuclearReaction();
lowEGammaModel->SetMaxEnergy(3.5*GeV);
G4PhotoNuclearProcess* thePhotoNuclearProcess = new G4PhotoNuclearProcess();
thePhotoNuclearProcess->RegisterMe(lowEGammaModel);
// bias the cross section
//
double thePhotoNuclearBias = control->GetPhotoNuclearBias();
if (thePhotoNuclearBias != 1.) {
thePhotoNuclearProcess->BiasCrossSectionByFactor(thePhotoNuclearBias);
// print out a warning if biasing
//
ErrMsg(warning) << "*** Biasing the photo-nuclear process by factor "
<< thePhotoNuclearBias << endmsg;
}
// high energy part
G4TheoFSGenerator* highEGammaModel = new G4TheoFSGenerator();
G4GeneratorPrecompoundInterface* preComModel =
new G4GeneratorPrecompoundInterface();
highEGammaModel->SetTransport(preComModel);
G4QGSModel* theStringModel =
new G4QGSModel;
G4QGSMFragmentation* fragModel = new G4QGSMFragmentation();
G4ExcitedStringDecay* stringDecay =
new G4ExcitedStringDecay(fragModel);
theStringModel->SetFragmentationModel(stringDecay);
highEGammaModel->SetHighEnergyGenerator(theStringModel);
highEGammaModel->SetMinEnergy(3.*GeV);
highEGammaModel->SetMaxEnergy(20.*GeV);
thePhotoNuclearProcess->RegisterMe(highEGammaModel);
G4ProcessManager* gamMan = G4Gamma::Gamma()->GetProcessManager();
AddDiscreteProcess(thePhotoNuclearProcess, gamMan, "gnuc",
GVertex::gammaNuclear );
//
// Electron-nuclear process
//
G4ElectroNuclearReaction* theElectronReaction =
new G4ElectroNuclearReaction();
G4ElectronNuclearProcess* theElectronNuclearProcess =
new G4ElectronNuclearProcess();
theElectronNuclearProcess->RegisterMe(theElectronReaction);
G4ProcessManager* electronMan = G4Electron::Electron()->GetProcessManager();
AddProcess(theElectronNuclearProcess, -1, -1, 4, electronMan, "enuc",
GVertex::electroNuclear );
// bias the cross section
//
G4double theElectroNuclearBias = control->GetElectroNuclearBias();
if (theElectroNuclearBias != 1.) {
theElectronNuclearProcess->BiasCrossSectionByFactor(theElectroNuclearBias);
// print out a warning if biasing
//
ErrMsg(warning) << "*** Biasing the electron-nuclear process by factor "
<< theElectroNuclearBias << endmsg;
}
//
// Positron-nuclear process
//
G4PositronNuclearProcess* thePositronNuclearProcess =
new G4PositronNuclearProcess();
thePositronNuclearProcess->RegisterMe(theElectronReaction);
G4ProcessManager* positronMan = G4Positron::Positron()->GetProcessManager();
AddProcess(thePositronNuclearProcess, -1, -1, 5, positronMan, "enuc",
GVertex::electroNuclear );
// bias the cross section
//
if (theElectroNuclearBias != 1.) {
thePositronNuclearProcess->BiasCrossSectionByFactor(theElectroNuclearBias);
ErrMsg(warning) << "*** Biasing the positron-nuclear process by factor "
<< theElectroNuclearBias << endmsg;
}
}
void
BgsPhysicsList::ConstructHad()
{
// One process handles hadronic elastic processes for all hadrons.
// However hadronic inelastic processes are unique to each hadron.
// For pi+ and pi- only, substitute pi-Nuclear cross sections
G4PiNuclearCrossSection* piNucCS = new G4PiNuclearCrossSection();
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
// *******
// * pi+ *
// *******
if (particleName == "pi+") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4PionPlusInelasticProcess *inel_process = new G4PionPlusInelasticProcess();
inel_process->AddDataSet(piNucCS);
// G4LEPionPlusInelastic *inel_model =
// new G4LEPionPlusInelastic();
inel_process->RegisterMe(bertini_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// *******
// * pi- *
// *******
} else if (particleName == "pi-") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4PionMinusInelasticProcess *inel_process = new G4PionMinusInelasticProcess();
inel_process->AddDataSet(piNucCS);
// G4LEPionMinusInelastic *inel_model =
// new G4LEPionMinusInelastic();
inel_process->RegisterMe(bertini_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// *******
// * K+ *
// *******
} else if (particleName == "kaon+") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4KaonPlusInelasticProcess *inel_process = new G4KaonPlusInelasticProcess();
G4LEKaonPlusInelastic *inel_model = new G4LEKaonPlusInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// *******
// * K- *
// *******
} else if (particleName == "kaon-") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4KaonMinusInelasticProcess *inel_process = new G4KaonMinusInelasticProcess();
G4LEKaonMinusInelastic *inel_model = new G4LEKaonMinusInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// *******
// * K0L *
// *******
} else if (particleName == "kaon0L") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4KaonZeroLInelasticProcess *inel_process = new G4KaonZeroLInelasticProcess();
G4LEKaonZeroLInelastic *inel_model = new G4LEKaonZeroLInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// *******
// * K0S *
// *******
} else if (particleName == "kaon0S") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4KaonZeroSInelasticProcess *inel_process =
new G4KaonZeroSInelasticProcess();
G4LEKaonZeroSInelastic *inel_model = new G4LEKaonZeroSInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// *******
// * p *
// *******
} else if (particleName == "proton") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4ProtonInelasticProcess *inel_process = new G4ProtonInelasticProcess();
// G4LEProtonInelastic *inel_model = new G4LEProtonInelastic();
inel_process->RegisterMe(bertini_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// *********
// * p-bar *
// *********
} else if (particleName == "anti_proton") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4AntiProtonInelasticProcess *inel_process =
new G4AntiProtonInelasticProcess();
G4LEAntiProtonInelastic *inel_model = new G4LEAntiProtonInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// *******
// * n *
// *******
} else if (particleName == "neutron") {
if (control->UseHPNeutrons()) {
G4cout << "High precision neutron models chosen" << G4endl;
putenv("NeutronHPCrossSections=/afs/slac/g/babar/simu/geant4_1/G4NDL3.7");
// Elastic process
G4HadronElasticProcess* el_process = new G4HadronElasticProcess();
// High precision model and data below 20 MeV
G4NeutronHPElastic* hpel_model = new G4NeutronHPElastic();
G4NeutronHPElasticData* el_data = new G4NeutronHPElasticData();
el_process->AddDataSet(el_data);
el_process->RegisterMe(hpel_model);
// LEP model above 20 MeV
G4LElastic* el_model = new G4LElastic();
el_model->SetMinEnergy(19.9*MeV);
el_process->RegisterMe(el_model);
pmanager->AddDiscreteProcess(el_process);
physicsRegistrar->Register(el_process,pmanager,"hade",
GVertex::hadronElastic);
// Inelastic process
G4NeutronInelasticProcess* inel_process =
new G4NeutronInelasticProcess();
// High precision model and data below 20 MeV
G4NeutronHPInelastic* hpinel_model = new G4NeutronHPInelastic();
G4NeutronHPInelasticData* hpinel_data = new G4NeutronHPInelasticData();
inel_process->AddDataSet(hpinel_data);
inel_process->RegisterMe(hpinel_model);
// Bertini model above 20 MeV
G4CascadeInterface* neutron_bertini = new G4CascadeInterface;
neutron_bertini->SetMinEnergy(19.9*MeV);
inel_process->RegisterMe(neutron_bertini);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// Capture process
G4HadronCaptureProcess* cap_process = new G4HadronCaptureProcess();
// High precision model and data below 20 MeV
G4NeutronHPCapture* hpcap_model = new G4NeutronHPCapture();
G4NeutronHPCaptureData* hpcap_data = new G4NeutronHPCaptureData();
cap_process->AddDataSet(hpcap_data);
cap_process->RegisterMe(hpcap_model);
// LEP model above 20 MeV
G4LCapture* cap_model = new G4LCapture();
cap_model->SetMinEnergy(19.9*MeV);
cap_process->RegisterMe(cap_model);
// Fission process
G4HadronFissionProcess* fis_process = new G4HadronFissionProcess();
// High precision model and data below 20 MeV
G4NeutronHPFission* hpfis_model = new G4NeutronHPFission();
G4NeutronHPFissionData* hpfis_data = new G4NeutronHPFissionData();
fis_process->AddDataSet(hpfis_data);
fis_process->RegisterMe(hpfis_model);
// LEP model above 20 MeV
G4LFission* fis_model = new G4LFission();
fis_model->SetMinEnergy(19.9*MeV);
fis_process->RegisterMe(fis_model);
} else {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4NeutronInelasticProcess *inel_process =
new G4NeutronInelasticProcess();
inel_process->RegisterMe(bertini_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
}
// *********
// * n-bar *
// *********
} else if (particleName == "anti_neutron") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4AntiNeutronInelasticProcess *inel_process =
new G4AntiNeutronInelasticProcess();
G4LEAntiNeutronInelastic *inel_model = new G4LEAntiNeutronInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// **********
// * lambda *
// **********
} else if (particleName == "lambda") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4LambdaInelasticProcess *inel_process = new G4LambdaInelasticProcess();
G4LELambdaInelastic *inel_model = new G4LELambdaInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// **************
// * lambda-bar *
// **************
} else if (particleName == "anti_lambda") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4AntiLambdaInelasticProcess *inel_process =
new G4AntiLambdaInelasticProcess();
G4LEAntiLambdaInelastic *inel_model = new G4LEAntiLambdaInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// **************
// * deuteron *
// **************
} else if (particleName == "deuteron") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4DeuteronInelasticProcess *inel_process =
new G4DeuteronInelasticProcess();
G4LEDeuteronInelastic *inel_model = new G4LEDeuteronInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// **************
// * triton *
// **************
} else if (particleName == "triton") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4TritonInelasticProcess *inel_process =
new G4TritonInelasticProcess();
G4LETritonInelastic *inel_model = new G4LETritonInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
// **************
// * alpha *
// **************
} else if (particleName == "alpha") {
// Elastic process
AddElasticHadronProcess(new G4HadronElasticProcess(), new G4LElastic(),
pmanager, "hade", GVertex::hadronElastic);
// Inelastic process
G4AlphaInelasticProcess *inel_process = new G4AlphaInelasticProcess();
G4LEAlphaInelastic *inel_model = new G4LEAlphaInelastic();
inel_process->RegisterMe(inel_model);
pmanager->AddDiscreteProcess(inel_process);
physicsRegistrar->Register(inel_process,pmanager,"hadi",
GVertex::hadronInelastic);
}
} // while
}
//
// ConstructNeutrinoGenocide
//
void BgsPhysicsList::ConstructNeutrinoGenocide()
{
BgsGenocide *genocide = new BgsGenocide();
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if (particleName == "nu_e" ||
particleName == "nu_mu" ||
particleName == "nu_tau" ||
particleName == "anti_nu_e" ||
particleName == "anti_nu_mu" ||
particleName == "anti_nu_tau" ) {
pmanager->AddProcess(genocide, -1, -1, 5);
physicsRegistrar->Register( genocide, pmanager, "neutrinoDeath",
GVertex::neutrino );
}
}
}
//
// ConstructIonFix
//
void BgsPhysicsList::ConstructIonFix()
{
BgsGenocide *genocide = new BgsGentleGenocide(control->GetIonEnergyCut()*MeV,
60);
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if ( particleName == "triton" ||
particleName == "alpha" ||
particleName == "proton" ||
particleName == "deuteron" ) {
pmanager->AddProcess(genocide, -1, -1, 5);
physicsRegistrar->Register( genocide, pmanager, "ionFix",
GVertex::minimumEnergy );
}
}
}
//
// ConstructNeutronFix
//
void BgsPhysicsList::ConstructNeutronFix()
{
BgsGenocide *genocide = new BgsGentleGenocide(control->GetNeutronEnergyCut()*MeV,0);
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if ( particleName == "neutron" ) {
pmanager->AddProcess(genocide, -1, -1, 1);
physicsRegistrar->Register( genocide, pmanager, "neutronFix",
GVertex::minimumEnergy );
}
}
}
// ****************************** Cuts ***************************************
void BgsPhysicsList::SetCuts()
{
// Set default cuts, all volumes
SetDefaultCutValue(control->GetDefaultRangeCut()*mm);
SetCutsWithDefault();
// Enable print out of cuts after tables are built
// This is now done in BgsRunAction
//
// if (verboseLevel > 1) DumpCutValuesTable();
}