DMXPhysicsList.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DMXPhysicsList.cc

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// --------------------------------------------------------------
//   GEANT 4 - Underground Dark Matter Detector Advanced Example
//
//      For information related to this code contact: Alex Howard
//      e-mail: alexander.howard@cern.ch
// --------------------------------------------------------------
// Comments
//
//                  Underground Advanced
//               by A. Howard and H. Araujo 
//                    (27th November 2001)
//
// PhysicsList program
//
// Modified:
//
// 14-02-03 Fix bugs in msc and hIon instanciation + cut per region
//
// 05-02-05 AH - changes to G4Decay - added is not short lived protection
//          and redefined particles to allow non-static creation
//          i.e. changed construction to G4MesonConstructor, G4BaryonConstructor
// 
// 23-10-09 LP - migrated EM physics from the LowEnergy processes (not supported) to 
//          the new G4Livermore model implementation. Results unchanged.
//
// --------------------------------------------------------------

#include <iomanip>  

#include "DMXPhysicsList.hh"

#include "globals.hh"
#include "G4SystemOfUnits.hh"
#include "G4ProcessManager.hh"
#include "G4ProcessVector.hh"

#include "G4ParticleDefinition.hh"
#include "G4ParticleWithCuts.hh"
#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"

#include "G4ios.hh"
#include "G4UserLimits.hh"

#include "G4MesonConstructor.hh"
#include "G4BaryonConstructor.hh"
#include "G4IonConstructor.hh"
#include "G4ShortLivedConstructor.hh"

#include "DMXMaxTimeCuts.hh"
#include "DMXMinEkineCuts.hh"
#include "G4StepLimiter.hh"

// gamma
#include "G4PhotoElectricEffect.hh"
#include "G4LivermorePhotoElectricModel.hh"

#include "G4ComptonScattering.hh"
#include "G4LivermoreComptonModel.hh"

#include "G4GammaConversion.hh"
#include "G4BetheHeitler5DModel.hh"

#include "G4RayleighScattering.hh" 
#include "G4LivermoreRayleighModel.hh"

// e-
#include "G4eMultipleScattering.hh"

#include "G4eIonisation.hh"
#include "G4LivermoreIonisationModel.hh"

#include "G4eBremsstrahlung.hh"
#include "G4UniversalFluctuation.hh"

// e+
#include "G4eIonisation.hh" 
#include "G4eBremsstrahlung.hh" 
#include "G4eplusAnnihilation.hh"

// alpha and GenericIon and deuterons, triton, He3:
//muon:
#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
#include "G4MuMultipleScattering.hh"
#include "G4MuonMinusCapture.hh"

//OTHERS:
#include "G4hIonisation.hh" 
#include "G4hMultipleScattering.hh"
#include "G4hBremsstrahlung.hh"
#include "G4ionIonisation.hh"
#include "G4IonParametrisedLossModel.hh"

//em process options to allow msc step-limitation to be switched off
#include "G4EmParameters.hh"
#include "G4VAtomDeexcitation.hh"
#include "G4UAtomicDeexcitation.hh"
#include "G4LossTableManager.hh"

#include "G4Scintillation.hh"
#include "G4OpAbsorption.hh"
#include "G4OpBoundaryProcess.hh"

// Elastic processes:
#include "G4HadronElasticProcess.hh"
#include "G4ChipsElasticModel.hh"
#include "G4ElasticHadrNucleusHE.hh"

// Inelastic processes:
#include "G4PionPlusInelasticProcess.hh"
#include "G4PionMinusInelasticProcess.hh"
#include "G4KaonPlusInelasticProcess.hh"
#include "G4KaonZeroSInelasticProcess.hh"
#include "G4KaonZeroLInelasticProcess.hh"
#include "G4KaonMinusInelasticProcess.hh"
#include "G4ProtonInelasticProcess.hh"
#include "G4AntiProtonInelasticProcess.hh"
#include "G4NeutronInelasticProcess.hh"
#include "G4AntiNeutronInelasticProcess.hh"
#include "G4DeuteronInelasticProcess.hh"
#include "G4TritonInelasticProcess.hh"
#include "G4AlphaInelasticProcess.hh"

// High energy FTFP model and Bertini cascade
#include "G4FTFModel.hh"
#include "G4LundStringFragmentation.hh"
#include "G4ExcitedStringDecay.hh"
#include "G4PreCompoundModel.hh"
#include "G4GeneratorPrecompoundInterface.hh"
#include "G4TheoFSGenerator.hh"
#include "G4CascadeInterface.hh"

// Cross sections
#include "G4VCrossSectionDataSet.hh"
#include "G4CrossSectionDataSetRegistry.hh"

#include "G4CrossSectionElastic.hh"
#include "G4CrossSectionInelastic.hh"
#include "G4BGGPionElasticXS.hh"
#include "G4BGGPionInelasticXS.hh"
#include "G4AntiNuclElastic.hh"

#include "G4CrossSectionInelastic.hh"
#include "G4CrossSectionPairGG.hh"
#include "G4BGGNucleonInelasticXS.hh"
#include "G4BGGNucleonElasticXS.hh"
#include "G4NeutronInelasticXS.hh"
#include "G4NeutronElasticXS.hh"
#include "G4ComponentAntiNuclNuclearXS.hh"
#include "G4ComponentGGNuclNuclXsc.hh"
#include "G4ComponentGGHadronNucleusXsc.hh"

#include "G4HadronElastic.hh"
#include "G4HadronCaptureProcess.hh"

// Neutron high-precision models: <20 MeV
#include "G4ParticleHPElastic.hh"
#include "G4ParticleHPElasticData.hh"
#include "G4ParticleHPCapture.hh"
#include "G4ParticleHPCaptureData.hh"
#include "G4ParticleHPInelastic.hh"
#include "G4ParticleHPInelasticData.hh"

// Stopping processes
#include "G4PiMinusAbsorptionBertini.hh"
#include "G4KaonMinusAbsorptionBertini.hh"
#include "G4AntiProtonAbsorptionFritiof.hh"
#include "G4HadronicParameters.hh"

#include "G4Decay.hh"
#include "G4RadioactiveDecayBase.hh"
#include "G4PhysicsListHelper.hh"
#include "G4NuclideTable.hh"
#include "G4NuclearLevelData.hh"

// Constructor /////////////////////////////////////////////////////////////
DMXPhysicsList::DMXPhysicsList() : G4VUserPhysicsList() 
{

  defaultCutValue     = 1.0*micrometer; //
  cutForGamma         = defaultCutValue;
  cutForElectron      = 1.0*nanometer;
  cutForPositron      = defaultCutValue;

  VerboseLevel = 1;
  OpVerbLevel = 0;

  //set a finer grid of the physic tables in order to improve precision
  //former LowEnergy models have 200 bins up to 100 GeV
  G4EmParameters* param = G4EmParameters::Instance();
  param->SetMaxEnergy(100*GeV);
  param->SetNumberOfBinsPerDecade(20);
  param->SetMscStepLimitType(fMinimal);
  param->SetFluo(true);
  param->SetPixe(true);
  param->SetAuger(true);

  G4EmParameters::Instance()->AddPhysics("World","G4RadioactiveDecay");
  G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
  deex->SetStoreICLevelData(true);
  deex->SetMaxLifeTime(G4NuclideTable::GetInstance()->GetThresholdOfHalfLife()
                       /std::log(2.));
  SetVerboseLevel(VerboseLevel);
}

// Destructor //////////////////////////////////////////////////////////////
DMXPhysicsList::~DMXPhysicsList() 
{;}

// Construct Particles /////////////////////////////////////////////////////
void DMXPhysicsList::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. 

  ConstructMyBosons();
  ConstructMyLeptons();
  ConstructMyHadrons();
  ConstructMyShortLiveds();
}

// construct Bosons://///////////////////////////////////////////////////
void DMXPhysicsList::ConstructMyBosons()
{
  // pseudo-particles
  G4Geantino::GeantinoDefinition();
  G4ChargedGeantino::ChargedGeantinoDefinition();
  
  // gamma
  G4Gamma::GammaDefinition();

  //OpticalPhotons
  G4OpticalPhoton::OpticalPhotonDefinition();

}

// construct Leptons://///////////////////////////////////////////////////
void DMXPhysicsList::ConstructMyLeptons()
{
  // leptons
  G4Electron::ElectronDefinition();
  G4Positron::PositronDefinition();
  G4MuonPlus::MuonPlusDefinition();
  G4MuonMinus::MuonMinusDefinition();

  G4NeutrinoE::NeutrinoEDefinition();
  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
  G4NeutrinoMu::NeutrinoMuDefinition();
  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
}

// construct Hadrons://///////////////////////////////////////////////////
void DMXPhysicsList::ConstructMyHadrons()
{
 //  mesons
  G4MesonConstructor mConstructor;
  mConstructor.ConstructParticle();

 //  baryons
  G4BaryonConstructor bConstructor;
  bConstructor.ConstructParticle();

 //  ions
  G4IonConstructor iConstructor;
  iConstructor.ConstructParticle();
}

// construct Shortliveds://///////////////////////////////////////////////////
void DMXPhysicsList::ConstructMyShortLiveds()
{
  G4ShortLivedConstructor slConstructor;
  slConstructor.ConstructParticle();
}

// Construct Processes //////////////////////////////////////////////////////
void DMXPhysicsList::ConstructProcess() 
{
  AddTransportation();

  ConstructEM();

  ConstructOp();

  ConstructHad();

  ConstructGeneral();
}

// Transportation ///////////////////////////////////////////////////////////
void DMXPhysicsList::AddTransportation() {
  
  G4VUserPhysicsList::AddTransportation();
  
  auto particleIterator=GetParticleIterator();
  particleIterator->reset();
  while( (*particleIterator)() ){
    G4ParticleDefinition* particle = particleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
    // time cuts for ONLY neutrons:
    if(particleName == "neutron") 
      pmanager->AddDiscreteProcess(new DMXMaxTimeCuts());
    // Energy cuts to kill charged (embedded in method) particles:
    pmanager->AddDiscreteProcess(new DMXMinEkineCuts());

    // Step limit applied to all particles:
    pmanager->AddDiscreteProcess(new G4StepLimiter);
  }         
}

// Electromagnetic Processes ////////////////////////////////////////////////
// all charged particles
void DMXPhysicsList::ConstructEM() {
  
  G4LossTableManager* man = G4LossTableManager::Instance();
  man->SetAtomDeexcitation(new G4UAtomicDeexcitation());

  auto particleIterator=GetParticleIterator();
  particleIterator->reset();
  while( (*particleIterator)() ){
    G4ParticleDefinition* particle = particleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
    G4String particleType = particle->GetParticleType();
    G4double charge = particle->GetPDGCharge();
    
    if (particleName == "gamma") 
      {
  //gamma
  G4RayleighScattering* theRayleigh = new G4RayleighScattering();
  pmanager->AddDiscreteProcess(theRayleigh);

  G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
  thePhotoElectricEffect->SetEmModel(new G4LivermorePhotoElectricModel());
  pmanager->AddDiscreteProcess(thePhotoElectricEffect);
  
  G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
  theComptonScattering->SetEmModel(new G4LivermoreComptonModel());
  pmanager->AddDiscreteProcess(theComptonScattering);
  
  G4GammaConversion* theGammaConversion = new G4GammaConversion();
  theGammaConversion->SetEmModel(new G4BetheHeitler5DModel());
  pmanager->AddDiscreteProcess(theGammaConversion);

      } 
    else if (particleName == "e-") 
      {
  //electron
  // process ordering: AddProcess(name, at rest, along step, post step)
  // Multiple scattering
  G4eMultipleScattering* msc = new G4eMultipleScattering();
  msc->SetStepLimitType(fUseDistanceToBoundary);
  pmanager->AddProcess(msc,-1, 1, -1);

  // Ionisation
  G4eIonisation* eIonisation = new G4eIonisation();
        G4VEmModel* theIoniLiv = new G4LivermoreIonisationModel();
        theIoniLiv->SetHighEnergyLimit(0.1*MeV); 
        eIonisation->AddEmModel(0, theIoniLiv, new G4UniversalFluctuation() );
  eIonisation->SetStepFunction(0.2, 100*um); //improved precision in tracking  
  pmanager->AddProcess(eIonisation,-1, 2, 2);
  
  // Bremsstrahlung
  G4eBremsstrahlung* eBremsstrahlung = new G4eBremsstrahlung();
  pmanager->AddProcess(eBremsstrahlung, -1,-3, 3);
      } 
    else if (particleName == "e+") 
      {
  //positron  
  G4eMultipleScattering* msc = new G4eMultipleScattering();
  msc->SetStepLimitType(fUseDistanceToBoundary);
  pmanager->AddProcess(msc,-1, 1, 1);
  
  // Ionisation
  G4eIonisation* eIonisation = new G4eIonisation();
  eIonisation->SetStepFunction(0.2, 100*um); //     
  pmanager->AddProcess(eIonisation,                 -1, 2, 2);

  //Bremsstrahlung (use default, no low-energy available)
  pmanager->AddProcess(new G4eBremsstrahlung(), -1,-1, 3);

  //Annihilation
  pmanager->AddProcess(new G4eplusAnnihilation(),0,-1, 4);      
      } 
    else if( particleName == "mu+" || 
       particleName == "mu-"    ) 
      {
  //muon  
  pmanager->AddProcess(new G4MuMultipleScattering,    -1, 1,-1);
  pmanager->AddProcess(new G4MuIonisation(),          -1, 2, 1);
  pmanager->AddProcess(new G4MuBremsstrahlung(),      -1,-1, 2);
  pmanager->AddProcess(new G4MuPairProduction(),      -1,-1, 3);
  if( particleName == "mu-" )
    pmanager->AddProcess(new G4MuonMinusCapture(), 0,-1,-1);
      } 
    else if (particleName == "proton" || 
       particleName == "pi+" || 
       particleName == "pi-")
      {
  //multiple scattering
  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, -1);
      
  //ionisation
  G4hIonisation* hIonisation = new G4hIonisation();
  hIonisation->SetStepFunction(0.2, 50*um);
  pmanager->AddProcess(hIonisation,               -1, 2, 1);      
  
  //bremmstrahlung
  pmanager->AddProcess(new G4hBremsstrahlung,     -1,-3, 2);
      }
    else if(particleName == "alpha"      ||
       particleName == "deuteron"   ||
       particleName == "triton"     ||
       particleName == "He3")
      {
  //multiple scattering
  pmanager->AddProcess(new G4hMultipleScattering,-1,1,-1);
  
  //ionisation
  G4ionIonisation* ionIoni = new G4ionIonisation();
  ionIoni->SetStepFunction(0.1, 20*um);
  pmanager->AddProcess(ionIoni,                  -1, 2,-2);
      }
    else if (particleName == "GenericIon")
      {
  // OBJECT may be dynamically created as either a GenericIon or nucleus
  // G4Nucleus exists and therefore has particle type nucleus
  // genericIon:
  
  //multiple scattering
  pmanager->AddProcess(new G4hMultipleScattering,-1,1,-1);

  //ionisation
  G4ionIonisation* ionIoni = new G4ionIonisation();
  ionIoni->SetEmModel(new G4IonParametrisedLossModel());
  ionIoni->SetStepFunction(0.1, 20*um);
  pmanager->AddProcess(ionIoni,                   -1, 2, 1);
      } 

    else if ((!particle->IsShortLived()) &&
       (charge != 0.0) && 
       (particle->GetParticleName() != "chargedgeantino")) 
      {
  //all others charged particles except geantino
        G4hMultipleScattering* aMultipleScattering = new G4hMultipleScattering();
        G4hIonisation* ahadronIon = new G4hIonisation();
  
  //multiple scattering
  pmanager->AddProcess(aMultipleScattering,-1,1,-1);

  //ionisation
  pmanager->AddProcess(ahadronIon,       -1,2,1);      
      }
  }
}

// Optical Processes ////////////////////////////////////////////////////////
void DMXPhysicsList::ConstructOp() 
{
  // default scintillation process
  G4Scintillation* theScintProcessDef = new G4Scintillation("Scintillation");
  // theScintProcessDef->DumpPhysicsTable();
  theScintProcessDef->SetTrackSecondariesFirst(true);
  theScintProcessDef->SetScintillationYieldFactor(1.0); //
  theScintProcessDef->SetScintillationExcitationRatio(0.0); //
  theScintProcessDef->SetVerboseLevel(OpVerbLevel);

  // scintillation process for alpha:
  G4Scintillation* theScintProcessAlpha = new G4Scintillation("Scintillation");
  // theScintProcessNuc->DumpPhysicsTable();
  theScintProcessAlpha->SetTrackSecondariesFirst(true);
  theScintProcessAlpha->SetScintillationYieldFactor(1.1);
  theScintProcessAlpha->SetScintillationExcitationRatio(1.0);
  theScintProcessAlpha->SetVerboseLevel(OpVerbLevel);

  // scintillation process for heavy nuclei
  G4Scintillation* theScintProcessNuc = new G4Scintillation("Scintillation");
  // theScintProcessNuc->DumpPhysicsTable();
  theScintProcessNuc->SetTrackSecondariesFirst(true);
  theScintProcessNuc->SetScintillationYieldFactor(0.2);
  theScintProcessNuc->SetScintillationExcitationRatio(1.0);
  theScintProcessNuc->SetVerboseLevel(OpVerbLevel);

  // optical processes
  G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption();
  //  G4OpRayleigh* theRayleighScatteringProcess = new G4OpRayleigh();
  G4OpBoundaryProcess* theBoundaryProcess = new G4OpBoundaryProcess();
  //  theAbsorptionProcess->DumpPhysicsTable();
  //  theRayleighScatteringProcess->DumpPhysicsTable();
  theAbsorptionProcess->SetVerboseLevel(OpVerbLevel);
  // theRayleighScatteringProcess->SetVerboseLevel(OpVerbLevel);
  theBoundaryProcess->SetVerboseLevel(OpVerbLevel);

  auto particleIterator=GetParticleIterator();
  particleIterator->reset();
  while( (*particleIterator)() )
    {
      G4ParticleDefinition* particle = particleIterator->value();
      G4ProcessManager* pmanager = particle->GetProcessManager();
      G4String particleName = particle->GetParticleName();
      if (theScintProcessDef->IsApplicable(*particle)) {
  //      if(particle->GetPDGMass() > 5.0*GeV) 
  if(particle->GetParticleName() == "GenericIon") {
    pmanager->AddProcess(theScintProcessNuc); // AtRestDiscrete
    pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxAtRest);
    pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxPostStep);
  }   
  else if(particle->GetParticleName() == "alpha") {
    pmanager->AddProcess(theScintProcessAlpha);
    pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxAtRest);
    pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxPostStep);
  }
  else {
    pmanager->AddProcess(theScintProcessDef);
    pmanager->SetProcessOrderingToLast(theScintProcessDef,idxAtRest);
    pmanager->SetProcessOrderingToLast(theScintProcessDef,idxPostStep);
  }   
      }
      
      if (particleName == "opticalphoton") {
  pmanager->AddDiscreteProcess(theAbsorptionProcess);
  //  pmanager->AddDiscreteProcess(theRayleighScatteringProcess);
  pmanager->AddDiscreteProcess(theBoundaryProcess);
      }
    }
}

// Hadronic processes ////////////////////////////////////////////////////////

void DMXPhysicsList::ConstructHad() 
{
  //Elastic models
  G4HadronElastic* elastic_lhep0 = new G4HadronElastic();
  G4ChipsElasticModel* elastic_chip = new G4ChipsElasticModel();
  G4ElasticHadrNucleusHE* elastic_he = new G4ElasticHadrNucleusHE(); 
  
  // Inelastic scattering
  const G4double theFTFMin0 =    0.0*GeV;
  const G4double theFTFMin1 =    3.0*GeV;
  const G4double theFTFMax = G4HadronicParameters::Instance()->GetMaxEnergy();
  const G4double theBERTMin0 =   0.0*GeV;
  const G4double theBERTMin1 =  19.0*MeV;
  const G4double theBERTMax =    6.0*GeV;
  const G4double theHPMin =      0.0*GeV;
  const G4double theHPMax =     20.0*MeV;

  G4FTFModel * theStringModel = new G4FTFModel;
  G4ExcitedStringDecay * theStringDecay = new G4ExcitedStringDecay( new G4LundStringFragmentation );
  theStringModel->SetFragmentationModel( theStringDecay );
  G4PreCompoundModel * thePreEquilib = new G4PreCompoundModel( new G4ExcitationHandler );
  G4GeneratorPrecompoundInterface * theCascade = new G4GeneratorPrecompoundInterface( thePreEquilib );

  G4TheoFSGenerator * theFTFModel0 = new G4TheoFSGenerator( "FTFP" );
  theFTFModel0->SetHighEnergyGenerator( theStringModel );
  theFTFModel0->SetTransport( theCascade );
  theFTFModel0->SetMinEnergy( theFTFMin0 );
  theFTFModel0->SetMaxEnergy( theFTFMax );

  G4TheoFSGenerator * theFTFModel1 = new G4TheoFSGenerator( "FTFP" );
  theFTFModel1->SetHighEnergyGenerator( theStringModel );
  theFTFModel1->SetTransport( theCascade );
  theFTFModel1->SetMinEnergy( theFTFMin1 );
  theFTFModel1->SetMaxEnergy( theFTFMax );

  G4CascadeInterface * theBERTModel0 = new G4CascadeInterface;
  theBERTModel0->SetMinEnergy( theBERTMin0 );
  theBERTModel0->SetMaxEnergy( theBERTMax );

  G4CascadeInterface * theBERTModel1 = new G4CascadeInterface;
  theBERTModel1->SetMinEnergy( theBERTMin1 );
  theBERTModel1->SetMaxEnergy( theBERTMax );

  G4VCrossSectionDataSet * theAntiNucleonData = new G4CrossSectionInelastic( new G4ComponentAntiNuclNuclearXS );
  G4ComponentGGNuclNuclXsc * ggNuclNuclXsec = new G4ComponentGGNuclNuclXsc();
  G4VCrossSectionDataSet * theGGNuclNuclData = new G4CrossSectionInelastic(ggNuclNuclXsec);
  G4VCrossSectionDataSet * theGGNNEl = new G4CrossSectionElastic(ggNuclNuclXsec);
  G4ComponentGGHadronNucleusXsc * ggHNXsec = new G4ComponentGGHadronNucleusXsc();
  G4VCrossSectionDataSet * theGGHNEl = new G4CrossSectionElastic(ggHNXsec);
  G4VCrossSectionDataSet * theGGHNInel = new G4CrossSectionInelastic(ggHNXsec);

  auto particleIterator=GetParticleIterator();
  particleIterator->reset();
  while ((*particleIterator)())
    {
      G4ParticleDefinition* particle = particleIterator->value();
      G4ProcessManager* pmanager = particle->GetProcessManager();
      G4String particleName = particle->GetParticleName();

      if (particleName == "pi+") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
          theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
          theElasticProcess->RegisterMe( elastic_he );
    pmanager->AddDiscreteProcess( theElasticProcess );
    //Inelastic scattering
    G4PionPlusInelasticProcess* theInelasticProcess = 
      new G4PionPlusInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );
  } 

      else if (particleName == "pi-") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
          theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
          theElasticProcess->RegisterMe( elastic_he );
    pmanager->AddDiscreteProcess( theElasticProcess );
    //Inelastic scattering
    G4PionMinusInelasticProcess* theInelasticProcess = 
      new G4PionMinusInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( new G4BGGPionInelasticXS( particle ) );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );    
    //Absorption
    pmanager->AddRestProcess(new G4PiMinusAbsorptionBertini, ordDefault);
  }
      else if (particleName == "kaon+") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
    theElasticProcess->AddDataSet( theGGHNEl );
          theElasticProcess->RegisterMe( elastic_lhep0 );
    pmanager->AddDiscreteProcess( theElasticProcess );
          // Inelastic scattering 
    G4KaonPlusInelasticProcess* theInelasticProcess = 
      new G4KaonPlusInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( theGGHNInel );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );
  }      
      else if (particleName == "kaon0S") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
    theElasticProcess->AddDataSet( theGGHNEl );
          theElasticProcess->RegisterMe( elastic_lhep0 );
    pmanager->AddDiscreteProcess( theElasticProcess );
          // Inelastic scattering  
    G4KaonZeroSInelasticProcess* theInelasticProcess = 
      new G4KaonZeroSInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( theGGHNInel );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );    
  }

      else if (particleName == "kaon0L") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
    theElasticProcess->AddDataSet( theGGHNEl );
          theElasticProcess->RegisterMe( elastic_lhep0 );
    pmanager->AddDiscreteProcess( theElasticProcess );
    // Inelastic scattering
    G4KaonZeroLInelasticProcess* theInelasticProcess = 
      new G4KaonZeroLInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( theGGHNInel );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 ); 
    pmanager->AddDiscreteProcess( theInelasticProcess );    
  }

      else if (particleName == "kaon-") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
    theElasticProcess->AddDataSet( theGGHNEl );
          theElasticProcess->RegisterMe( elastic_lhep0 );
    pmanager->AddDiscreteProcess( theElasticProcess );
          // Inelastic scattering
    G4KaonMinusInelasticProcess* theInelasticProcess = 
      new G4KaonMinusInelasticProcess("inelastic"); 
          theInelasticProcess->AddDataSet( theGGHNInel );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );
    pmanager->AddRestProcess(new G4KaonMinusAbsorptionBertini, ordDefault);
  }

      else if (particleName == "proton") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
          theElasticProcess->AddDataSet( new G4BGGNucleonElasticXS( G4Proton::Proton() ) );
          theElasticProcess->RegisterMe( elastic_chip );
    pmanager->AddDiscreteProcess( theElasticProcess );
    // Inelastic scattering
    G4ProtonInelasticProcess* theInelasticProcess = 
      new G4ProtonInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( new G4BGGNucleonInelasticXS( G4Proton::Proton() ) );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );
  }
      else if (particleName == "anti_proton") 
  {
    // Elastic scattering
          const G4double elastic_elimitAntiNuc = 100.0*MeV;
          G4AntiNuclElastic* elastic_anuc = new G4AntiNuclElastic();
          elastic_anuc->SetMinEnergy( elastic_elimitAntiNuc );
          G4CrossSectionElastic* elastic_anucxs = new G4CrossSectionElastic( elastic_anuc->GetComponentCrossSection() );
          G4HadronElastic* elastic_lhep2 = new G4HadronElastic();
          elastic_lhep2->SetMaxEnergy( elastic_elimitAntiNuc );
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
          theElasticProcess->AddDataSet( elastic_anucxs );
          theElasticProcess->RegisterMe( elastic_lhep2 );
          theElasticProcess->RegisterMe( elastic_anuc );
    pmanager->AddDiscreteProcess( theElasticProcess );
    // Inelastic scattering
    G4AntiProtonInelasticProcess* theInelasticProcess = 
      new G4AntiProtonInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( theAntiNucleonData );
    theInelasticProcess->RegisterMe( theFTFModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );
    // Absorption
    pmanager->AddRestProcess(new G4AntiProtonAbsorptionFritiof, ordDefault);
  }
      else if (particleName == "neutron") {
  // elastic scattering
  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
        theElasticProcess->AddDataSet(new G4NeutronElasticXS());
        G4HadronElastic* elastic_neutronChipsModel = new G4ChipsElasticModel();
  elastic_neutronChipsModel->SetMinEnergy( 19.0*MeV );
        theElasticProcess->RegisterMe( elastic_neutronChipsModel );
  G4ParticleHPElastic * theElasticNeutronHP = new G4ParticleHPElastic;
        theElasticNeutronHP->SetMinEnergy( theHPMin );
        theElasticNeutronHP->SetMaxEnergy( theHPMax );
  theElasticProcess->RegisterMe( theElasticNeutronHP );
  theElasticProcess->AddDataSet( new G4ParticleHPElasticData );
  pmanager->AddDiscreteProcess( theElasticProcess );
  // inelastic scattering   
  G4NeutronInelasticProcess* theInelasticProcess =
    new G4NeutronInelasticProcess("inelastic");
  theInelasticProcess->AddDataSet( new G4NeutronInelasticXS() );
  theInelasticProcess->RegisterMe( theFTFModel1 );
        theInelasticProcess->RegisterMe( theBERTModel1 );
  G4ParticleHPInelastic * theNeutronInelasticHPModel = new G4ParticleHPInelastic;
        theNeutronInelasticHPModel->SetMinEnergy( theHPMin );
        theNeutronInelasticHPModel->SetMaxEnergy( theHPMax );
  theInelasticProcess->RegisterMe( theNeutronInelasticHPModel );
  theInelasticProcess->AddDataSet( new G4ParticleHPInelasticData );
  pmanager->AddDiscreteProcess(theInelasticProcess);
  // capture
  G4HadronCaptureProcess* theCaptureProcess =
    new G4HadronCaptureProcess;
  G4ParticleHPCapture * theLENeutronCaptureModel = new G4ParticleHPCapture;
  theLENeutronCaptureModel->SetMinEnergy(theHPMin);
  theLENeutronCaptureModel->SetMaxEnergy(theHPMax);
  theCaptureProcess->RegisterMe(theLENeutronCaptureModel);
  theCaptureProcess->AddDataSet( new G4ParticleHPCaptureData);
  pmanager->AddDiscreteProcess(theCaptureProcess);
      }
      else if (particleName == "anti_neutron") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
    theElasticProcess->AddDataSet( theGGHNEl );
          theElasticProcess->RegisterMe( elastic_lhep0 );
    pmanager->AddDiscreteProcess( theElasticProcess );
          // Inelastic scattering (include annihilation on-fly)
    G4AntiNeutronInelasticProcess* theInelasticProcess = 
      new G4AntiNeutronInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( theAntiNucleonData );
    theInelasticProcess->RegisterMe( theFTFModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );    
  }
      else if (particleName == "deuteron") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
    theElasticProcess->AddDataSet( theGGNNEl );
          theElasticProcess->RegisterMe( elastic_lhep0 );
    pmanager->AddDiscreteProcess( theElasticProcess );
          // Inelastic scattering
    G4DeuteronInelasticProcess* theInelasticProcess = 
      new G4DeuteronInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( theGGNuclNuclData );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );
  }
      else if (particleName == "triton") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
    theElasticProcess->AddDataSet( theGGNNEl );
          theElasticProcess->RegisterMe( elastic_lhep0 );
    pmanager->AddDiscreteProcess( theElasticProcess );
          // Inelastic scattering
    G4TritonInelasticProcess* theInelasticProcess = 
      new G4TritonInelasticProcess("inelastic");
    theInelasticProcess->AddDataSet( theGGNuclNuclData );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );
  }
      else if (particleName == "alpha") 
  {
    // Elastic scattering
          G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
    theElasticProcess->AddDataSet( theGGNNEl );
          theElasticProcess->RegisterMe( elastic_lhep0 );
    pmanager->AddDiscreteProcess( theElasticProcess );
          // Inelastic scattering
    G4AlphaInelasticProcess* theInelasticProcess = 
      new G4AlphaInelasticProcess("inelastic");  
          theInelasticProcess->AddDataSet( theGGNuclNuclData );
    theInelasticProcess->RegisterMe( theFTFModel1 );
          theInelasticProcess->RegisterMe( theBERTModel0 );
    pmanager->AddDiscreteProcess( theInelasticProcess );
  }
    }
}

// Decays ///////////////////////////////////////////////////////////////////
void DMXPhysicsList::ConstructGeneral() {

  // Add Decay Process
  G4Decay* theDecayProcess = new G4Decay();
  auto particleIterator=GetParticleIterator();
  particleIterator->reset();
  while( (*particleIterator)() )
    {
      G4ParticleDefinition* particle = particleIterator->value();
      G4ProcessManager* pmanager = particle->GetProcessManager();
      
      if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived()) 
  { 
    pmanager ->AddProcess(theDecayProcess);
    // set ordering for PostStepDoIt and AtRestDoIt
    pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
    pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
  }
    }

  // Declare radioactive decay to the GenericIon in the IonTable.
  G4LossTableManager* man = G4LossTableManager::Instance();
  G4VAtomDeexcitation* ad = man->AtomDeexcitation();
  if(!ad) {
    G4EmParameters::Instance()->SetAugerCascade(true);
    ad = new G4UAtomicDeexcitation();
    man->SetAtomDeexcitation(ad);
    ad->InitialiseAtomicDeexcitation();
  }

  G4PhysicsListHelper::GetPhysicsListHelper()->
    RegisterProcess(new G4RadioactiveDecayBase(), G4GenericIon::GenericIon());
}

// Cuts /////////////////////////////////////////////////////////////////////
void DMXPhysicsList::SetCuts() 
{
  
  if (verboseLevel >1)
    G4cout << "DMXPhysicsList::SetCuts:";
  
  if (verboseLevel>0){
    G4cout << "DMXPhysicsList::SetCuts:";
    G4cout << "CutLength : " 
     << G4BestUnit(defaultCutValue,"Length") << G4endl;
  }

  //special for low energy physics
  G4double lowlimit=250*eV;  
  G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(lowlimit,100.*GeV);

  // 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+");
  
  if (verboseLevel>0) DumpCutValuesTable();
}