mt01physicslist.cc

利用linux环境下蒙特卡罗计算工具geant4实现伽玛射线与探测器乘积能量的模拟过程

//
// ********************************************************************
// * DISCLAIMER                                                       *
// *                                                                  *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in:                      *
// *   http://cern.ch/geant4/license                                  *
// *                                                                  *
// * 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.                                                             *
// *                                                                  *
// * This  code  implementation is the  intellectual property  of the *
// * GEANT4 collaboration.                                            *
// * By copying,  distributing  or modifying the Program (or any work *
// * based  on  the Program)  you indicate  your  acceptance of  this *
// * statement, and all its terms.                                    *
// ********************************************************************
//
//
// --------------------------------------------------------------
//   GEANT 4 - Underground Dark Matter Detector Advanced Example
//
//      For information related to this code contact: Alex Howard
//      e-mail: [e-mail]
// --------------------------------------------------------------
// 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
//
// --------------------------------------------------------------

#include "mt01PhysicsList.hh"

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

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

#include "G4ios.hh"
#include <iomanip>  

#include "G4UserLimits.hh"


// Constructor /////////////////////////////////////////////////////////////
mt01PhysicsList::mt01PhysicsList() : G4VUserPhysicsList() 
{

  defaultCutValue     =10.*keV;
    // 1.0*micrometer; //
  cutForGamma         = defaultCutValue;
  cutForElectron      =10.*keV;
    //1.0*nanometer;
  cutForPositron      = defaultCutValue;

  VerboseLevel = 1;
  OpVerbLevel = 0;

  SetVerboseLevel(VerboseLevel);
}


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


// Construct Particles /////////////////////////////////////////////////////
void mt01PhysicsList::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();
  ConstructMyMesons();
  ConstructMyBaryons();
  ConstructMyIons();
  ConstructMyShortLiveds();

}


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

  //OpticalPhotons
  G4OpticalPhoton::OpticalPhotonDefinition();

}


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

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


// construct Mesons://///////////////////////////////////////////////////
void mt01PhysicsList::ConstructMyMesons()
{
 //  mesons
  G4PionPlus::PionPlusDefinition();
  G4PionMinus::PionMinusDefinition();
  G4PionZero::PionZeroDefinition();
  G4KaonPlus::KaonPlusDefinition();
  G4KaonMinus::KaonMinusDefinition();
  G4Eta::EtaDefinition();
  G4EtaPrime::EtaPrimeDefinition();
  G4KaonZero::KaonZeroDefinition();
  G4AntiKaonZero::AntiKaonZeroDefinition();
  G4KaonZeroLong::KaonZeroLongDefinition();
  G4KaonZeroShort::KaonZeroShortDefinition();
}


// construct Baryons://///////////////////////////////////////////////////
void mt01PhysicsList::ConstructMyBaryons()
{
//  barions
  G4Proton::ProtonDefinition();
  G4AntiProton::AntiProtonDefinition();
  G4Neutron::NeutronDefinition();
  G4AntiNeutron::AntiNeutronDefinition();
}


// construct Ions://///////////////////////////////////////////////////
void mt01PhysicsList::ConstructMyIons()
{
//  Ions
  G4Deuteron::DeuteronDefinition();
  G4Triton::TritonDefinition();
  G4He3::He3Definition();
  G4Alpha::AlphaDefinition();
  G4GenericIon::GenericIonDefinition();
}

// construct Shortliveds://///////////////////////////////////////////////////
void mt01PhysicsList::ConstructMyShortLiveds()
{
  // ShortLiveds
  ;
}




// Construct Processes //////////////////////////////////////////////////////
void mt01PhysicsList::ConstructProcess() 
{

  AddTransportation();

  ConstructEM();

  //ConstructOp();

  ConstructHad();

  ConstructGeneral();

}


// Transportation ///////////////////////////////////////////////////////////
#include "mt01MaxTimeCuts.hh"
#include "mt01MinEkineCuts.hh"

void mt01PhysicsList::AddTransportation() {
  
  G4VUserPhysicsList::AddTransportation();
  
  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
    // time cuts for ONLY neutrons:
    if(particleName == "neutron") 
    pmanager->AddDiscreteProcess(new mt01MaxTimeCuts());
    // Energy cuts to kill charged (embedded in method) particles:
    pmanager->AddDiscreteProcess(new mt01MinEkineCuts());
  }		      
}


// Electromagnetic Processes ////////////////////////////////////////////////
// all charged particles
#include "G4MultipleScattering.hh"

// gamma
#include "G4LowEnergyRayleigh.hh" 
#include "G4LowEnergyPhotoElectric.hh"
#include "G4LowEnergyCompton.hh"  
#include "G4LowEnergyGammaConversion.hh" 


// e-
#include "G4LowEnergyIonisation.hh" 
#include "G4LowEnergyBremsstrahlung.hh" 

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


// alpha and GenericIon and deuterons, triton, He3:
#include "G4hLowEnergyIonisation.hh"
#include "G4EnergyLossTables.hh"
// hLowEnergyIonisation uses Ziegler 1988 as the default


//muon:
#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
#include "G4MuonMinusCaptureAtRest.hh"

//OTHERS:
//#include "G4hIonisation.hh" // standard hadron ionisation

void mt01PhysicsList::ConstructEM() {

// processes:

  G4LowEnergyPhotoElectric* lowePhot = new G4LowEnergyPhotoElectric();
  G4LowEnergyIonisation* loweIon  = new G4LowEnergyIonisation();
  G4LowEnergyBremsstrahlung* loweBrem = new G4LowEnergyBremsstrahlung();

  // note LowEIon uses proton as basis for its data-base, therefore
  // cannot specify different LowEnergyIonisation models for different
  // particles, but can change model globally for Ion, Alpha and Proton.


  //fluorescence apply specific cut for fluorescence from photons, electrons
  //and bremsstrahlung photons:
  G4double fluorcut = 250*eV;
  lowePhot->SetCutForLowEnSecPhotons(fluorcut);
  loweIon->SetCutForLowEnSecPhotons(fluorcut);
  loweBrem->SetCutForLowEnSecPhotons(fluorcut);
  
  // setting tables explicitly for electronic stopping power
  //  ahadronLowEIon->SetElectronicStoppingPowerModel
  //  (G4GenericIon::GenericIonDefinition(), "ICRU_R49p") ;
  //  ahadronLowEIon->SetElectronicStoppingPowerModel
  //  (G4Proton::ProtonDefinition(), "ICRU_R49p") ;

  // Switch off the Barkas and Bloch corrections
  //  ahadronLowEIon->SetBarkasOff();


  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
    G4String particleType = particle->GetParticleType();
    G4double charge = particle->GetPDGCharge();
    
    if (particleName == "gamma") 
      {
	//gamma
	pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh());
	pmanager->AddDiscreteProcess(lowePhot);
	pmanager->AddDiscreteProcess(new G4LowEnergyCompton());
	pmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion());
      } 
    else if (particleName == "e-") 
      {
	//electron
	// process ordering: AddProcess(name, at rest, along step, post step)
	// -1 = not implemented, then ordering
        G4MultipleScattering* aMultipleScattering = new G4MultipleScattering();
	pmanager->AddProcess(aMultipleScattering,     -1, 1, 1);
	pmanager->AddProcess(loweIon,                 -1, 2, 2);
	pmanager->AddProcess(loweBrem,                -1,-1, 3);
      } 
    else if (particleName == "e+") 
      {
	//positron
        G4MultipleScattering* aMultipleScattering = new G4MultipleScattering();
	pmanager->AddProcess(aMultipleScattering,     -1, 1, 1);
	pmanager->AddProcess(new G4eIonisation(),     -1, 2, 2);
	pmanager->AddProcess(new G4eBremsstrahlung(), -1,-1, 3);
	pmanager->AddProcess(new G4eplusAnnihilation(),0,-1, 4);      
      } 
    else if( particleName == "mu+" || 
	     particleName == "mu-"    ) 
      {
	//muon  
        G4MultipleScattering* aMultipleScattering = new G4MultipleScattering();
	pmanager->AddProcess(aMultipleScattering,           -1, 1, 1);
	pmanager->AddProcess(new G4MuIonisation(),          -1, 2, 2);
	pmanager->AddProcess(new G4MuBremsstrahlung(),      -1,-1, 3);
	pmanager->AddProcess(new G4MuPairProduction(),      -1,-1, 4);
	if( particleName == "mu-" )
	  pmanager->AddProcess(new G4MuonMinusCaptureAtRest(), 0,-1,-1);
      } 
    else if (particleName == "proton"     ||
	     particleName == "alpha"      ||
	     particleName == "deuteron"   ||
	     particleName == "triton"     ||
	     particleName == "He3"        ||
	     particleName == "GenericIon" || 
	     (particleType == "nucleus" && charge != 0)) 
      {
	// OBJECT may be dynamically created as either a GenericIon or nucleus
	// G4Nucleus exists and therefore has particle type nucleus
	// genericIon:
        G4MultipleScattering* aMultipleScattering = new G4MultipleScattering();
        G4hLowEnergyIonisation* ahadronLowEIon = new G4hLowEnergyIonisation();
	pmanager->AddProcess(aMultipleScattering,-1,1,1);
	pmanager->AddProcess(ahadronLowEIon,-1,2,2); 
        // ahadronLowEIon->SetNuclearStoppingOff() ;
	//        ahadronLowEIon->SetNuclearStoppingPowerModel("ICRU_R49") ;
	//        ahadronLowEIon->SetNuclearStoppingOn() ;
  
        //fluorescence switch off for hadrons (for now) PIXE:
        ahadronLowEIon->SetFluorescence(false);
      } 
    else if ((!particle->IsShortLived()) &&
	     (charge != 0.0) && 
	     (particle->GetParticleName() != "chargedgeantino")) 
      {
	//all others charged particles except geantino
        G4MultipleScattering* aMultipleScattering = new G4MultipleScattering();
        G4hLowEnergyIonisation* ahadronLowEIon = new G4hLowEnergyIonisation();
	pmanager->AddProcess(aMultipleScattering,-1,1,1);
	pmanager->AddProcess(ahadronLowEIon,       -1,2,2);      
	//      pmanager->AddProcess(new G4hIonisation(),       -1,2,2);      
      }
    
  }
}


// Optical Processes ////////////////////////////////////////////////////////
#include "G4Scintillation.hh"
#include "G4OpAbsorption.hh"
//#include "G4OpRayleigh.hh"
#include "G4OpBoundaryProcess.hh"

void mt01PhysicsList::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();