G4TankSimulatorOG/G4TankPhysicsList.cc

Go to the documentation of this file.

#include "G4TankPhysicsList.h"

#include <G4ParticleDefinition.hh>
#include <G4ParticleTypes.hh>
#include <G4ParticleWithCuts.hh>
#include <G4ParticleTable.hh>
#include <G4Material.hh>
#include <G4MaterialTable.hh>
#include <G4ProcessManager.hh>
#include <G4ProcessVector.hh>

#warning "GEANT4 version" G4VERSION_NUMBER
#include <G4MuonMinusCapture.hh>

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

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

#include <G4Decay.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 <G4OpAbsorption.hh>
#include <G4OpRayleigh.hh>
//#include <G4OpBoundaryProcess.hh>

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

#include "G4TankOpBoundaryProcess.h"

#include "G4TankFastCerenkov.h"
#include "G4TankSimulator.h"

#include <iostream>

using namespace std;
using namespace G4TankSimulatorOG;


G4TankPhysicsList::G4TankPhysicsList(const bool fastCerenkov) :
  G4VUserPhysicsList(),
  fFastCerenkov(fastCerenkov)
{
}


void
G4TankPhysicsList::ConstructParticle()
{
  ConstructBosons();
  ConstructLeptons();
  ConstructMesons();
  ConstructBaryons();
}


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


void
G4TankPhysicsList::ConstructLeptons()
{
  G4LeptonConstructor::ConstructParticle();
}


void
G4TankPhysicsList::ConstructMesons()
{
  G4MesonConstructor::ConstructParticle();
}


void
G4TankPhysicsList::ConstructBaryons()
{
  G4BaryonConstructor::ConstructParticle();
}


void
G4TankPhysicsList::ConstructProcess()
{
  AddTransportation();
  ConstructGeneral();
  ConstructEM();
  ConstructOp();
}


void
G4TankPhysicsList::ConstructGeneral()
{   
  G4Decay* const decayProcess = new G4Decay();
  
  G4ProcessManager& mumManager =
    *G4MuonMinus::MuonMinusDefinition()->GetProcessManager();  
  mumManager.AddDiscreteProcess(decayProcess);
  mumManager.SetProcessOrdering(decayProcess, idxPostStep);
  mumManager.SetProcessOrdering(decayProcess, idxAtRest);
  
  G4ProcessManager& mupManager =
    *G4MuonPlus::MuonPlusDefinition()->GetProcessManager();  
  mupManager.AddDiscreteProcess(decayProcess);
  mupManager.SetProcessOrdering(decayProcess, idxPostStep);
  mupManager.SetProcessOrdering(decayProcess, idxAtRest);
}


void
G4TankPhysicsList::ConstructEM()
{
  theParticleIterator->reset();
  while ((*theParticleIterator)()) {
    const G4ParticleDefinition& particle = *theParticleIterator->value();
    G4ProcessManager& pManager = *particle.GetProcessManager();
    const G4String& particleName = particle.GetParticleName();
     
    if (particleName == "gamma") {
      // gamma         
      pManager.AddDiscreteProcess(new G4PhotoElectricEffect);
      pManager.AddDiscreteProcess(new G4ComptonScattering);
      pManager.AddDiscreteProcess(new G4GammaConversion);
  
    } else if (particleName == "e-") {
      //electron
#if (G4VERSION_NUMBER>=930)
      pManager.AddProcess(new G4eMultipleScattering, -1, 1, 1);
#else
      pManager.AddProcess(new G4MultipleScattering, -1, 1, 1);
#endif
      pManager.AddProcess(new G4eIonisation,        -1, 2, 2);
      pManager.AddProcess(new G4eBremsstrahlung,    -1, 3, 3);

    } else if (particleName == "e+") {
      //positron
#if (G4VERSION_NUMBER>=930)
      pManager.AddProcess(new G4eMultipleScattering, -1,  1, 1);
#else
      pManager.AddProcess(new G4MultipleScattering, -1,  1, 1);
#endif
      pManager.AddProcess(new G4eIonisation,        -1,  2, 2);
      pManager.AddProcess(new G4eBremsstrahlung,    -1,  3, 3);
      pManager.AddProcess(new G4eplusAnnihilation,   0, -1, 4);

    } else if (particleName == "mu+" || 
               particleName == "mu-") {
      //muon  
#if (G4VERSION_NUMBER>=930)
      pManager.AddProcess(new G4MuMultipleScattering, -1, 1, 1);
#else
      pManager.AddProcess(new G4MultipleScattering, -1, 1, 1);
#endif
      pManager.AddProcess(new G4MuIonisation,       -1, 2, 2);
      pManager.AddProcess(new G4MuBremsstrahlung,   -1, 3, 3);
      pManager.AddProcess(new G4MuPairProduction,   -1, 4, 4);       
      if (particleName == "mu-" &&
          G4TankSimulator::MuCaptureOn()) {
        pManager.AddProcess(new G4MuonMinusCaptureAtRest);
      }

    } else if (particle.GetPDGCharge() &&
               particle.GetParticleName() != "chargedgeantino") {
        // all others charged particles except geantino
#if (G4VERSION_NUMBER>=930)
        pManager.AddProcess(new G4hMultipleScattering, -1, 1, 1);
#else
        pManager.AddProcess(new G4MultipleScattering, -1, 1, 1);
#endif
        pManager.AddProcess(new G4hIonisation,        -1, 2, 2);
    }
  }

  if (!fFastCerenkov) {

    G4Cerenkov* const cerenkov = new G4Cerenkov("Cerenkov");
    // 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;
    cerenkov->SetTrackSecondariesFirst(true);
    cerenkov->SetMaxNumPhotonsPerStep(maxNumPhotons);
 
    theParticleIterator->reset();
    while ((*theParticleIterator)()) {
      G4ParticleDefinition& particle = *theParticleIterator->value();
      G4ProcessManager& pManager = *particle.GetProcessManager();
      if (cerenkov->IsApplicable(particle)) {
        pManager.AddProcess(cerenkov);
        pManager.SetProcessOrdering(cerenkov, idxPostStep);
      }
    }

  } else {

    G4TankFastCerenkov* const cerenkov = new G4TankFastCerenkov("Cerenkov");
    // 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;
    cerenkov->SetTrackSecondariesFirst(true);
    cerenkov->SetMaxNumPhotonsPerStep(maxNumPhotons);

    G4Electron::Electron()->GetProcessManager()->AddContinuousProcess(cerenkov);
  
    G4Positron::Positron()->GetProcessManager()->AddContinuousProcess(cerenkov);

    G4MuonMinus::MuonMinusDefinition()->GetProcessManager()->AddContinuousProcess(cerenkov);

    G4MuonPlus::MuonPlusDefinition()->GetProcessManager()->AddContinuousProcess(cerenkov);

  }
}


void
G4TankPhysicsList::ConstructOp()
{  
  G4OpAbsorption* const absorption = new G4OpAbsorption();
  G4OpRayleigh* const rayleighScattering = new G4OpRayleigh(); 
  //  G4OpBoundaryProcess *boundary = new G4OpBoundaryProcess();
  G4TankOpBoundaryProcess* const boundary = new G4TankOpBoundaryProcess();
  
  // ATTENTION : Isn't this set in the DetectorConstruction as well?
  const G4OpticalSurfaceModel model = unified;
  boundary->SetModel(model);
  
  G4ProcessManager& pManager =
    *G4OpticalPhoton::OpticalPhotonDefinition()->GetProcessManager();
  pManager.AddDiscreteProcess(absorption);
  pManager.AddDiscreteProcess(rayleighScattering);
  pManager.AddDiscreteProcess(boundary);
}


void
G4TankPhysicsList::SetCuts()
{ 
  defaultCutValue = 0.1*mm;
  
  const double cutForGamma    = 30*mm;
  const double cutForElectron = 0.5*mm;
  const double cutForPositron = 1e-5*mm;
  
  // set cut values for gamma at first and for e- second and next for e+,
  // because some processes for e+/e- need cut values for gamma 
  SetCutValue(cutForGamma, "gamma");
  SetCutValue(cutForElectron, "e-");
  SetCutValue(cutForPositron, "e+");

  DumpCutValuesTable();
}