DefaultPhysicsList.cc

Go to the documentation of this file.

#include <G4Version.hh>

#include "DefaultPhysicsList.h"

#include <utl/ErrorLogger.h>

#include <G4UImanager.hh>
#include <G4ParticleDefinition.hh>
#include <G4ParticleTypes.hh>
#include <G4ParticleWithCuts.hh>
#include <G4ParticleTable.hh>
#include <G4ProcessManager.hh>

#if (G4VERSION_NUMBER >= 930)
#  include <G4hMultipleScattering.hh>
#  include <G4eMultipleScattering.hh>
#  include <G4MuMultipleScattering.hh>
#else
#  include "G4MultipleScattering.hh"
#endif
#include <G4Decay.hh>

#include <G4ComptonScattering.hh>
#include <G4GammaConversion.hh>
#include <G4PhotoElectricEffect.hh>

#include <G4eIonisation.hh>
#include <G4eBremsstrahlung.hh>
#include <G4eplusAnnihilation.hh>

#include <G4MuIonisation.hh>
#include <G4MuBremsstrahlung.hh>
#include <G4MuPairProduction.hh>

#include <G4hIonisation.hh>

#include <G4Cerenkov.hh>

#include <G4IonConstructor.hh>
#include <G4BaryonConstructor.hh>
#include <G4MesonConstructor.hh>
#include <G4LeptonConstructor.hh>

// using Geant4 boundary process...
//#include <G4OpBoundaryProcess.hh>
#include <G4OpAbsorption.hh>
#include <G4OpRayleigh.hh>

#if (G4VERSION_NUMBER >= 920)
#  include "tls/DummyProcessAtRest.h"
#endif

using namespace tls;


DefaultPhysicsList::DefaultPhysicsList(const int verbosity) :
  //G4VUserPhysicsList(verbosity),
  fVerbosity(verbosity)
{
  G4VUserPhysicsList::SetVerboseLevel(fVerbosity);
  //defaultCutValue = 0.7*mm;
  fUImanager = G4UImanager::GetUIpointer();
}


DefaultPhysicsList::~DefaultPhysicsList()
{
  INFO("destructing physics list");
}


void
DefaultPhysicsList::ConstructParticle()
{
  // pseudo-particle: useful for testing
  //G4Geantino::GeantinoDefinition();
  // gamma
  G4Gamma::GammaDefinition();
  // optical photon
  G4OpticalPhoton::OpticalPhotonDefinition();

  G4LeptonConstructor::ConstructParticle();
  G4MesonConstructor::ConstructParticle();
  G4BaryonConstructor::ConstructParticle();
}


void
DefaultPhysicsList::ConstructProcess()
{
  this->ConstructOpticalProcesses();
  this->ConstructEMProcesses();
  this->ConstructDecayProcesses();
}


void
DefaultPhysicsList::ConstructDecayProcesses()
{
  auto decayProcess = new G4Decay();

  auto& mumProcMan = *G4MuonMinus::MuonMinusDefinition()->GetProcessManager();
  mumProcMan.AddDiscreteProcess(decayProcess);
  mumProcMan.SetProcessOrdering(decayProcess, idxPostStep);
  mumProcMan.SetProcessOrdering(decayProcess, idxAtRest);

  auto& mupProcMan = *G4MuonPlus::MuonPlusDefinition()->GetProcessManager();
  mupProcMan.AddDiscreteProcess(decayProcess);
  mupProcMan.SetProcessOrdering(decayProcess, idxPostStep);
  mupProcMan.SetProcessOrdering(decayProcess, idxAtRest);

#if (G4VERSION_NUMBER >= 920)
  /* DummyProcessAtRest is only needed if all other AtRestDoIt processes are inactivated
  (Scintillation, Decay, MuonMinusCaptureAtRest) for a given particle.
  It is inactivated by default. Remind that changing the activation while in
  G4State_PreInit or G4State_Init has no effect, i.e. has to happen after the first
  call to G4RunManagerKernel::InitializePhysics() */
  auto dummyProcess = new DummyProcessAtRest("DummyProcessAtRest");
  mumProcMan.AddProcess(dummyProcess);
  mumProcMan.SetProcessOrdering(dummyProcess, idxAtRest);
  mupProcMan.AddProcess(dummyProcess);
  mupProcMan.SetProcessOrdering(dummyProcess, idxAtRest);
#endif
}


void
DefaultPhysicsList::ConstructOpticalProcesses()
{
  auto& procMan = *G4OpticalPhoton::OpticalPhotonDefinition()->GetProcessManager();
  procMan.AddDiscreteProcess(new G4OpAbsorption());
  procMan.AddDiscreteProcess(new G4OpRayleigh());
}


void
DefaultPhysicsList::ConstructEMProcesses()
{
  auto& partIt = *G4ParticleTable::GetParticleTable()->GetIterator();
  for (partIt.reset(); partIt(); ) {
    auto& particle = *partIt.value();
    auto& procMan = *particle.GetProcessManager();
    const G4String& particleName = particle.GetParticleName();

    if (particleName == "gamma") {
      procMan.AddDiscreteProcess(new G4PhotoElectricEffect);
      procMan.AddDiscreteProcess(new G4ComptonScattering);
      procMan.AddDiscreteProcess(new G4GammaConversion);
    } else if (particleName == "e-") {
      procMan.AddProcess(new G4eMultipleScattering, -1, 1, 1);
      procMan.AddProcess(new G4eIonisation,         -1, 2, 2);
      procMan.AddProcess(new G4eBremsstrahlung,     -1, 3, 3);
    } else if (particleName == "e+") {
      procMan.AddProcess(new G4eMultipleScattering, -1, 1, 1);
      procMan.AddProcess(new G4eIonisation,         -1,  2, 2);
      procMan.AddProcess(new G4eBremsstrahlung,     -1,  3, 3);
      procMan.AddProcess(new G4eplusAnnihilation,    0, -1, 4);
    } else if (particleName == "mu+" || particleName == "mu-") {
      procMan.AddProcess(new G4MuMultipleScattering, -1, 1, 1);
      procMan.AddProcess(new G4MuIonisation,       -1, 2, 2);
      procMan.AddProcess(new G4MuBremsstrahlung,   -1, 3, 3);
      procMan.AddProcess(new G4MuPairProduction,   -1, 4, 4);
    } else if (particle.GetPDGCharge() &&
               particle.GetParticleName() != "chargedgeantino") {
      // all others charged particles except geantino
      procMan.AddProcess(new G4hMultipleScattering, -1, 1, 1);
      procMan.AddProcess(new G4hIonisation,         -1, 2, 2);
    }
  }

  // NOTE: At low energy, where dE/dX is large for the muon, it is best
  // to set the following number to something low
  const G4int maxNumPhotons = 3;
  auto cerenkov = new G4Cerenkov("StandardCerenkov");
  cerenkov->SetTrackSecondariesFirst(true);
  cerenkov->SetMaxNumPhotonsPerStep(maxNumPhotons);
  for (partIt.reset(); partIt(); ) {
    auto& particle = *partIt.value();
    auto& procMan = *particle.GetProcessManager();
    if (cerenkov->IsApplicable(particle)) {
      procMan.AddProcess(cerenkov);
      procMan.SetProcessOrdering(cerenkov, idxPostStep);
    }
  }
}


bool
DefaultPhysicsList::ActivateCustomProcesses()
{
  if (false) {
    // not needed as long as decay process is activated
    G4cerr << "activating DummyProcessAtRest for mu+" << G4endl;
    fUImanager->ApplyCommand("/process/activate DummyProcessAtRest mu+");
    G4cerr << "activating DummyProcessAtRest for mu-" << G4endl;
    fUImanager->ApplyCommand("/process/activate DummyProcessAtRest mu-");
  }
  return false;
}


bool
DefaultPhysicsList::InactivateCustomProcesses()
{
  G4cerr << "de-activating DummyProcessAtRest for mu+" << G4endl;
  fUImanager->ApplyCommand("/process/inactivate DummyProcessAtRest mu+");
  G4cerr << "de-activating DummyProcessAtRest for mu-" << G4endl;
  fUImanager->ApplyCommand("/process/inactivate DummyProcessAtRest mu-");
  return true;
}