///////////////////////////////////////////////////////////////////////// // G4Tutorial: // // PhysicsList.hh // // Activating a list of physics constructors, which create physics // processes relevant for the simulation, and setting production thresholds. // ///////////////////////////////////////////////////////////////////////// #include "PhysicsList.hh" #include "PhysicsListMessenger.hh" #include "Particles.hh" #include "PhysicsDecay.hh" #include "PhysicsEMElectronEEDL.hh" #include "PhysicsEMElectronPenelope.hh" #include "PhysicsEMElectronStandard.hh" #include "PhysicsEMPositronStandard.hh" #include "PhysicsEMPhotonEPDL.hh" #include "PhysicsEMPhotonPenelope.hh" #include "PhysicsEMPhotonStandard.hh" #include "PhysicsEMHadronIonLowEnergy.hh" #include "PhysicsHEHadronIonLElastic.hh" #include "PhysicsHIPionBinary.hh" #include "PhysicsHIIonBinary.hh" #include "PhysicsHIProtonNeutron.hh" #include "G4VPhysicsConstructor.hh" #include "G4RegionStore.hh" #include "G4Region.hh" #include "G4ProductionCuts.hh" G4VPhysicsConstructor* PhysicsList::decay = 0; G4VPhysicsConstructor* PhysicsList::emElectron = 0; G4VPhysicsConstructor* PhysicsList::emPositron = 0; G4VPhysicsConstructor* PhysicsList::emPhoton = 0; G4VPhysicsConstructor* PhysicsList::emHadronIon = 0; G4VPhysicsConstructor* PhysicsList::heHadronIon = 0; G4VPhysicsConstructor* PhysicsList::hiPion = 0; G4VPhysicsConstructor* PhysicsList::hiIon = 0; G4VPhysicsConstructor* PhysicsList::hiProtonNeutron = 0; PhysicsList::PhysicsList(G4String detRegName) : G4VModularPhysicsList(), detectorRegionName(detRegName) { // Initialization of value for lower production threshold (the data member // defaultCutValue is defined in the parent class) defaultCutValue = 0.001 * mm; // Instantiation of messenger messenger = new PhysicsListMessenger(this); // Registering particles RegisterPhysics(new Particles()); } PhysicsList::~PhysicsList() { // Deleting instance of messenger delete messenger; } void PhysicsList::RegisterPhysConstructor(const G4String& constrName) { // In this example, physics constructors can be registered for several // categories: // A) Electromagnetic (EM) physics for electrons // B) Electromagnetic (EM) physics for positrons // C) Electromagnetic (EM) physics for photons // D) Electromagnetic (EM) physics for charged hadrons and ions // E) Hadronic elastic (HE) physics for hadrons and ions // F) Hadronic inelastic (HI) physics for pions // G) Hadronic inelastic (HI) physics for ions // H) Hadronic inelastic (HI) physics for protons and neutrons // I) Decay for unstable particles // // For each category, one or more physics constructors were implemented. // Different physics constructors for a given category contain alternative // physics models: // A) PhysicsEMElectronEEDL, PhysicsEMElectronPenelope, // PhysicsEMElectronStandard // B) PhysicsEMPhotonEPDL, PhysicsEMPhotonPenelope, PhysicsEMPhotonStandard // C) PhysicsEMPositronStandard // D) PhysicsEMHadronIonLowEnergy // E) PhysicsHEHadronIonLElastic // F) PhysicsHIPionBinary // G) PhysicsHIIonBinary // H) PhysicsHIProtonNeutron // I) PhysicsDecay // (Have a look at the these classes for more details) // // NOTE: The user can only select ONE physics constructor of each category. // // The selected physics constructors for each category are assigned to the // data members emElectron, emPositron and emPhoton, respectively, and they // are registered using the RegisterPhysics member function: // ******************** // *** EM: Electron *** // ******************** if(emElectron == 0) { // The condition is only true if no electron physics // constructor was yet instantiated // If the physics constructor name is matching one of the strings, // the corresponding physics constructor is instantiated and assigned // to the data member emElectron if(constrName == "EM-Electron-EEDL") emElectron = new PhysicsEMElectronEEDL(); if(constrName == "EM-Electron-Penelope") emElectron = new PhysicsEMElectronPenelope(); if(constrName == "EM-Electron-Standard") emElectron = new PhysicsEMElectronStandard(); // The physics constructor is finally registered (The "if" statement is // required for cases where the constructor name was not matching any // option) if(emElectron) RegisterPhysics(emElectron); } // ******************** // *** EM: Positron *** // ******************** if(emPositron == 0) { // The condition is only true if no positron physics // constructor was yet instantiated // If the physics constructor name is matching the string, // the corresponding physics constructor is instantiated and assigned // to the data member emPositron if(constrName == "EM-Positron-Standard") emPositron = new PhysicsEMPositronStandard(); // The physics constructor is finally registered (The "if" statement is // required for cases where the constructor name was not matching any // option) if(emPositron) RegisterPhysics(emPositron); } // ****************** // *** EM: Photon *** // ****************** if(emPhoton == 0) { // The condition is only true if no photon physics // constructor was yet instantiated // If the physics constructor name is matching one of the three strings, // the corresponding physics constructor is instantiated and assigned // to the data member emPhoton if(constrName == "EM-Photon-EPDL") emPhoton = new PhysicsEMPhotonEPDL(); if(constrName == "EM-Photon-Penelope") emPhoton = new PhysicsEMPhotonPenelope(); if(constrName == "EM-Photon-Standard") emPhoton = new PhysicsEMPhotonStandard(); // The physics constructor is finally registered (The "if" statement is // required for cases where the constructor name was not matching any // option) if(emPhoton) RegisterPhysics(emPhoton); } // *********************** // *** EM: Hadron, Ion *** // *********************** if(emHadronIon == 0) { // The condition is only true if no hadron-ion // physics constructor for EM interactions // was yet instantiated // If the physics constructor name is matching the string, // the corresponding physics constructor is instantiated and assigned // to the data member emHadronIon if(constrName == "EM-HadronIon-LowEnergy") emHadronIon = new PhysicsEMHadronIonLowEnergy(); // The physics constructor is finally registered (The "if" statement is // required for cases where the constructor name was not matching any // option) if(emHadronIon) RegisterPhysics(emHadronIon); } // ************************************* // *** Hadronic elastic: Hadron, Ion *** // ************************************* if(heHadronIon == 0) { // The condition is only true if no hadron-ion // physics constructor for elastic hadronic // interactions was yet instantiated // If the physics constructor name is matching the string, // the corresponding physics constructor is instantiated and assigned // to the data member emHadronIon if(constrName == "HE-HadronIon-LElastic") heHadronIon = new PhysicsHEHadronIonLElastic(); // The physics constructor is finally registered (The "if" statement is // required for cases where the constructor name was not matching any // option) if(heHadronIon) RegisterPhysics(heHadronIon); } // ******************************** // *** Hadronic inelastic: Pion *** // ******************************** if(hiPion == 0) { // The condition is only true if no pion // physics constructor for inelastic hadronic // interactions was yet instantiated // If the physics constructor name is matching the string, // the corresponding physics constructor is instantiated and assigned // to the data member emPion if(constrName == "HI-Pion-Binary") hiPion = new PhysicsHIPionBinary(); // The physics constructor is finally registered (The "if" statement is // required for cases where the constructor name was not matching any // option) if(hiPion) RegisterPhysics(hiPion); } // ******************************** // *** Hadronic inelastic: Ion *** // ******************************** if(hiIon == 0) { // The condition is only true if no ion // physics constructor for inelastic hadronic // interactions was yet instantiated // If the physics constructor name is matching the string, // the corresponding physics constructor is instantiated and assigned // to the data member emIon if(constrName == "HI-Ion-Binary") hiIon = new PhysicsHIIonBinary(); // The physics constructor is finally registered (The "if" statement is // required for cases where the constructor name was not matching any // option) if(hiIon) RegisterPhysics(hiIon); } // ******************************************* // *** Hadronic inelastic: Proton, Neutron *** // ******************************************* if(hiProtonNeutron == 0) { // The condition is only true if no proton-neutron // physics constructor for inelastic hadronic // interactions was yet instantiated // If the physics constructor name is matching the string, // the corresponding physics constructor is instantiated and assigned // to the data member emProtonNeutron if(constrName == "HI-ProtonNeutron-inelastic") hiProtonNeutron = new PhysicsHIProtonNeutron(); // The physics constructor is finally registered (The "if" statement is // required for cases where the constructor name was not matching any // option) if(hiProtonNeutron) RegisterPhysics(hiProtonNeutron); } // ********************************* // *** Decay: Unstable particles *** // ********************************* if(decay == 0) { // The condition is only true if the particle decay // physics constructor was not yet instantiated // If the physics constructor name is matching the string, // the corresponding physics constructor is instantiated and assigned // to the data member decay if(constrName == "Decay") decay = new PhysicsDecay(); // The physics constructor is finally registered (The "if" statement is // required for cases where the constructor name was not matching any // option) if(decay) RegisterPhysics(decay); } } void PhysicsList::SetCuts() { // In order to allow the use of production cuts as low as 250 eV (NOTE: // 250 eV is the smallest possible cut for the Livermore library based EM // models for electrons and photons) the limits of the physics table must // be set accordingly (default of the lower value is 990 eV): G4double lowerProdLimit = 250.0 * eV; G4double upperProdLimit = 1.0 * GeV; G4ProductionCutsTable::GetProductionCutsTable() -> SetEnergyRange(lowerProdLimit, upperProdLimit); // The default cut value, given by the data member defaultCutValue, is set // for all particles (electrons, positrons, photons) for the default region, // which is the world region: SetCutsWithDefault(); // To define production cuts for the detector, a G4ProductionCuts object is // first instantiated and the lower cut value is set: G4ProductionCuts* prodCutsDetector = new G4ProductionCuts; prodCutsDetector -> SetProductionCut(defaultCutValue); // The next step is to retrieve the G4Region instance, which represents the // detector, from the region store (by using the detector region name, which // was stored in a data member of the current class): G4Region* detectorRegion = G4RegionStore::GetInstance() -> GetRegion(detectorRegionName); // Then the G4ProductionCuts instance is assigned to the detector region: // (check that detectorRegion is a valid pointer before, otherwise the // program will crash) if (detectorRegion) detectorRegion -> SetProductionCuts(prodCutsDetector); // The production cuts are finally printed: DumpCutValuesTable(); } void PhysicsList::SetProdThreshold(double cut) { // The default cut value is set (NOTE: defaultCutValue is a data member // defined in the parent class) if(cut > 0.0 * mm) { defaultCutValue = cut; } }