gorad/src/GRPhysicsList.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //  Gorad (Geant4 Open-source Radiation Analysis and Design)
0027 //
0028 //  Author : Makoto Asai (SLAC National Accelerator Laboratory)
0029 //
0030 //  Development of Gorad is funded by NASA Johnson Space Center (JSC)
0031 //  under the contract NNJ15HK11B.
0032 //
0033 // ********************************************************************
0034 //
0035 // GRPhysicsList.cc
0036 //   Gorad Physics List
0037 //
0038 // History
0039 //   September 8th, 2020 : first implementation
0040 //
0041 // ********************************************************************
0042
0043 #include "GRPhysicsList.hh"
0044 #include "GRPhysicsListMessenger.hh"
0045 #include "G4SystemOfUnits.hh"
0046 #include "G4PhysListFactory.hh"
0047
0048 GRPhysicsList::GRPhysicsList()
0049 : PLName("FTFP_BERT"), physList(nullptr),
0050   EM_opt("Op_0"), Had_opt("FTFP_BERT"),
0051   addHP(false), addRDM(false), addRMC(false), addOptical(false),
0052   stepLimit_opt(-1)
0053 {
0054   factory = nullptr;
0055   messenger = new GRPhysicsListMessenger(this);
0056
0057   globalCuts[0] = 0.7*mm; // e-
0058   globalCuts[1] = 0.7*mm; // e+
0059   globalCuts[2] = 0.7*mm; // gamma
0060   globalCuts[3] = 0.7*mm; // proton
0061 }
0062
0063 GRPhysicsList::~GRPhysicsList()
0064 {
0065   delete factory;
0066   if(physList) delete physList;
0067   delete messenger;
0068 }
0069
0070 void GRPhysicsList::ConstructParticle()
0071 {
0072   if(!physList) GeneratePL();
0073   physList->ConstructParticle();
0074 }
0075
0076 void GRPhysicsList::ConstructProcess()
0077 {
0078   if(!physList) GeneratePL();
0079   physList->ConstructProcess();
0080 }
0081
0082 #include "G4Region.hh"
0083 #include "G4ProductionCuts.hh"
0084
0085 void GRPhysicsList::SetCuts()
0086 {
0087   if(!physList) GeneratePL();
0088   physList->SetCutValue(globalCuts[2],"gamma"); // gamma should be defined first!
0089   physList->SetCutValue(globalCuts[0],"e-");
0090   physList->SetCutValue(globalCuts[1],"e+");
0091   physList->SetCutValue(globalCuts[3],"proton");
0092 }
0093
0094 void GRPhysicsList::SetGlobalCuts(G4double val)
0095 {
0096   for(G4int i=0; i<4; i++)
0097   { SetGlobalCut(i,val); }
0098   if(physList) SetCuts();
0099 }
0100
0101 void GRPhysicsList::SetGlobalCut(G4int i, G4double val)
0102 {
0103   globalCuts[i] = val;
0104   if(physList) SetCuts();
0105 }
0106
0107 void GRPhysicsList::GeneratePLName()
0108 {
0109   G4String plname = Had_opt;
0110   if(addHP && Had_opt != "Shielding") plname += "_HP";
0111
0112   G4String EMopt = "";
0113   if(EM_opt=="Op_1") EMopt = "_EMV";
0114   else if(EM_opt=="Op_3") EMopt = "_EMY";
0115   else if(EM_opt=="Op_4") EMopt = "_EMZ";
0116   else if(EM_opt=="LIV")  EMopt = "_LIV";
0117   else if(EM_opt=="LIV_Pol") G4cout << "EM option <LIV_Pol> is under development." << G4endl;
0118   plname += EMopt;
0119
0120   auto valid = factory->IsReferencePhysList(plname);
0121   if(valid)
0122   { PLName = plname; }
0123   else
0124   {
0125     G4ExceptionDescription ed;
0126     ed << "Physics List <" << plname << "> is not a valid reference physics list.";
0127     G4Exception("GRPhysicsList::GeneratePLName()","GRPHYS0001",
0128                 FatalException,ed);
0129   }
0130 }
0131
0132 #include "G4RadioactiveDecayPhysics.hh"
0133 #include "G4OpticalPhysics.hh"
0134 #include "G4StepLimiterPhysics.hh"
0135 #include "G4ParallelWorldPhysics.hh"
0136 #include "G4GenericBiasingPhysics.hh"
0137 #include "GRParallelWorldPhysics.hh"
0138 #include "G4ProcessTable.hh"
0139 #include "G4EmParameters.hh"
0140 #include "G4HadronicParameters.hh"
0141
0142 void GRPhysicsList::GeneratePL()
0143 {
0144   if(physList) return;
0145
0146   factory = new G4PhysListFactory();
0147   GeneratePLName();
0148   physList = factory->GetReferencePhysList(PLName);
0149   G4cout << "Creating " << PLName << " physics list ################ " << applyGeomImpBias << G4endl;
0150
0151   if(addRDM && Had_opt != "Shielding")
0152   { physList->RegisterPhysics(new G4RadioactiveDecayPhysics());
0153     G4cout << "Adding G4RadioactiveDecayPhysics ################ " << G4endl; }
0154
0155   if(addRMC)
0156   { G4cout << "Reverse Monte Calro option is under development." << G4endl; }
0157
0158   if(stepLimit_opt>=0)
0159   { physList->RegisterPhysics(new G4StepLimiterPhysics());
0160     G4cout << "Adding G4StepLimiterPhysics ################ " << G4endl; }
0161
0162   if(addOptical) // Optical processes should come last!
0163   { physList->RegisterPhysics(new G4OpticalPhysics());
0164     G4cout << "Adding G4OpticalPhysics ################ " << G4endl; }
0165
0166   if(applyGeomImpBias) // Geometry Importance Biasing with parallel world
0167   {
0168     physList->RegisterPhysics(new GRParallelWorldPhysics("GeomBias",false));
0169     G4cout << "Adding G4GenericBiasingPhysics for GeomBias ################ " << G4endl;
0170   }
0171
0172   G4int verbose = G4ProcessTable::GetProcessTable()->GetVerboseLevel();
0173   physList->SetVerboseLevel(verbose);
0174   G4EmParameters::Instance()->SetVerbose(verbose);
0175   G4HadronicParameters::Instance()->SetVerboseLevel(verbose);
0176 }
0177
0178 #include "G4RegionStore.hh"
0179
0180 G4Region* GRPhysicsList::FindRegion(const G4String& reg) const
0181 {
0182   auto store = G4RegionStore::GetInstance();
0183   return store->GetRegion(reg);
0184 }
0185
0186 G4Region* GRPhysicsList::SetLocalCut(const G4String& reg,G4int i,G4double val)
0187 {
0188   auto regPtr = FindRegion(reg);
0189   if(!regPtr) return regPtr;
0190
0191   auto cuts = regPtr->GetProductionCuts();
0192   if(!cuts)
0193   {
0194     cuts = new G4ProductionCuts();
0195     regPtr->SetProductionCuts(cuts);
0196   }
0197
0198   cuts->SetProductionCut(val,i);
0199   return regPtr;
0200 }
0201
0202 G4double GRPhysicsList::GetLocalCut(const G4String& reg,G4int i) const
0203 {
0204   auto regPtr = FindRegion(reg);
0205   G4double val = -1.0;
0206   if(regPtr)
0207   {
0208     auto cuts = regPtr->GetProductionCuts();
0209     if(cuts) val = cuts->GetProductionCut(i);
0210   }
0211   return val;
0212 }
0213
0214 #include "G4UserLimits.hh"
0215
0216 G4Region* GRPhysicsList::SetLocalStepLimit(const G4String& reg,G4double val)
0217 {
0218   auto regPtr = FindRegion(reg);
0219   if(!regPtr) return regPtr;
0220
0221   auto uLim = regPtr->GetUserLimits();
0222   if(!uLim)
0223   {
0224     uLim = new G4UserLimits(val);
0225     regPtr->SetUserLimits(uLim);
0226   }
0227   else
0228   { uLim->SetMaxAllowedStep(val); }
0229   return regPtr;
0230 }
0231
0232 #include "G4Track.hh"
0233 G4double GRPhysicsList::GetLocalStepLimit(const G4String& reg) const
0234 {
0235   static G4Track dummyTrack;
0236   auto regPtr = FindRegion(reg);
0237   G4double val = -1.0;
0238   if(regPtr)
0239   {
0240     auto uLim = regPtr->GetUserLimits();
0241     if(uLim) val = uLim->GetMaxAllowedStep(dummyTrack);
0242   }
0243   return val;
0244 }
0245
0246 void GRPhysicsList::SetGlobalStepLimit(G4double val)
0247 { SetLocalStepLimit("DefaultRegionForTheWorld",val); }
0248
0249 G4double GRPhysicsList::GetGlobalStepLimit() const
0250 { return GetLocalStepLimit("DefaultRegionForTheWorld"); }
0251
0252