field04/src/F04PhysicsList.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 //
0027 /// \file field/field04/src/F04PhysicsList.cc
0028 /// \brief Implementation of the F04PhysicsList class
0029 //
0030
0031 #include "F04PhysicsList.hh"
0032
0033 #include "F04PhysicsListMessenger.hh"
0034
0035 #include "G4LossTableManager.hh"
0036 #include "G4OpticalPhysics.hh"
0037 #include "G4ParticleTable.hh"
0038 #include "G4ParticleTypes.hh"
0039 #include "G4ProcessManager.hh"
0040 #include "G4StepLimiterPhysics.hh"
0041
0042 // #include "G4PhysListFactory.hh"
0043 #include "F04StepMax.hh"
0044 #include "FTFP_BERT.hh"
0045 #include "QGSP_BERT.hh"
0046
0047 #include "G4AutoDelete.hh"
0048 #include "G4DecayTable.hh"
0049 #include "G4DecayWithSpin.hh"
0050 #include "G4Electron.hh"
0051 #include "G4Gamma.hh"
0052 #include "G4MuMinusCapturePrecompound.hh"
0053 #include "G4MuonDecayChannelWithSpin.hh"
0054 #include "G4MuonMinusCapture.hh"
0055 #include "G4MuonRadiativeDecayChannelWithSpin.hh"
0056 #include "G4PionDecayMakeSpin.hh"
0057 #include "G4Positron.hh"
0058 #include "G4ProcessTable.hh"
0059 #include "G4SystemOfUnits.hh"
0060
0061 G4ThreadLocal F04StepMax* F04PhysicsList::fStepMaxProcess = nullptr;
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 F04PhysicsList::F04PhysicsList(G4String physName) : G4VModularPhysicsList()
0066 {
0067   G4LossTableManager::Instance();
0068
0069   defaultCutValue = 1. * mm;
0070
0071   fMessenger = new F04PhysicsListMessenger(this);
0072
0073   SetVerboseLevel(1);
0074
0075   //    G4PhysListFactory factory;
0076   G4VModularPhysicsList* phys = nullptr;
0077   if (physName == "QGSP_BERT") {
0078     phys = new QGSP_BERT;
0079   }
0080   else {
0081     phys = new FTFP_BERT;
0082   }
0083
0084   //    if (factory.IsReferencePhysList(physName))
0085   //       phys =factory.GetReferencePhysList(physName);
0086
0087   // Physics List is defined via environment variable PHYSLIST
0088   //    if (!phys) phys = factory.ReferencePhysList();
0089
0090   if (!phys)
0091     G4Exception("F04PhysicsList::F04PhysicsList", "InvalidSetup", FatalException,
0092                 "PhysicsList does not exist");
0093
0094   for (G4int i = 0;; ++i) {
0095     auto elem = const_cast<G4VPhysicsConstructor*>(phys->GetPhysics(i));
0096     if (elem == nullptr) break;
0097     G4cout << "RegisterPhysics: " << elem->GetPhysicsName() << G4endl;
0098     RegisterPhysics(elem);
0099   }
0100
0101   RegisterPhysics(new G4StepLimiterPhysics());
0102   RegisterPhysics(new G4OpticalPhysics());
0103 }
0104
0105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0106
0107 F04PhysicsList::~F04PhysicsList()
0108 {
0109   delete fMessenger;
0110
0111   // delete fStepMaxProcess;
0112 }
0113
0114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0115
0116 void F04PhysicsList::ConstructParticle()
0117 {
0118   G4VModularPhysicsList::ConstructParticle();
0119
0120   G4GenericIon::GenericIonDefinition();
0121
0122   auto muonPlusDecayTable = new G4DecayTable();
0123   muonPlusDecayTable->Insert(new G4MuonDecayChannelWithSpin("mu+", 0.986));
0124   muonPlusDecayTable->Insert(new G4MuonRadiativeDecayChannelWithSpin("mu+", 0.014));
0125   G4MuonPlus::MuonPlusDefinition()->SetDecayTable(muonPlusDecayTable);
0126
0127   auto muonMinusDecayTable = new G4DecayTable();
0128   muonMinusDecayTable->Insert(new G4MuonDecayChannelWithSpin("mu-", 0.986));
0129   muonMinusDecayTable->Insert(new G4MuonRadiativeDecayChannelWithSpin("mu-", 0.014));
0130   G4MuonMinus::MuonMinusDefinition()->SetDecayTable(muonMinusDecayTable);
0131 }
0132
0133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0134
0135 void F04PhysicsList::ConstructProcess()
0136 {
0137   G4VModularPhysicsList::ConstructProcess();
0138
0139   fStepMaxProcess = new F04StepMax();
0140   G4AutoDelete::Register(fStepMaxProcess);
0141
0142   auto decayWithSpin = new G4DecayWithSpin();
0143
0144   G4ProcessTable* processTable = G4ProcessTable::GetProcessTable();
0145
0146   G4VProcess* decay;
0147   decay = processTable->FindProcess("Decay", G4MuonPlus::MuonPlus());
0148
0149   G4ProcessManager* pmanager;
0150   pmanager = G4MuonPlus::MuonPlus()->GetProcessManager();
0151
0152   if (pmanager) {
0153     if (decay) pmanager->RemoveProcess(decay);
0154     pmanager->AddProcess(decayWithSpin);
0155     // set ordering for PostStepDoIt and AtRestDoIt
0156     pmanager->SetProcessOrdering(decayWithSpin, idxPostStep);
0157     pmanager->SetProcessOrdering(decayWithSpin, idxAtRest);
0158   }
0159
0160   decay = processTable->FindProcess("Decay", G4MuonMinus::MuonMinus());
0161
0162   pmanager = G4MuonMinus::MuonMinus()->GetProcessManager();
0163
0164   if (pmanager) {
0165     if (decay) pmanager->RemoveProcess(decay);
0166     pmanager->AddProcess(decayWithSpin);
0167     // set ordering for PostStepDoIt and AtRestDoIt
0168     pmanager->SetProcessOrdering(decayWithSpin, idxPostStep);
0169     pmanager->SetProcessOrdering(decayWithSpin, idxAtRest);
0170   }
0171
0172   G4VProcess* process = processTable->FindProcess("muMinusCaptureAtRest", G4MuonMinus::MuonMinus());
0173
0174   if (pmanager) {
0175     if (process) pmanager->RemoveProcess(process);
0176     process = new G4MuonMinusCapture(new G4MuMinusCapturePrecompound());
0177     pmanager->AddRestProcess(process);
0178   }
0179
0180   auto poldecay = new G4PionDecayMakeSpin();
0181
0182   decay = processTable->FindProcess("Decay", G4PionPlus::PionPlus());
0183
0184   pmanager = G4PionPlus::PionPlus()->GetProcessManager();
0185
0186   if (pmanager) {
0187     if (decay) pmanager->RemoveProcess(decay);
0188     pmanager->AddProcess(poldecay);
0189     // set ordering for PostStepDoIt and AtRestDoIt
0190     pmanager->SetProcessOrdering(poldecay, idxPostStep);
0191     pmanager->SetProcessOrdering(poldecay, idxAtRest);
0192   }
0193
0194   decay = processTable->FindProcess("Decay", G4PionMinus::PionMinus());
0195
0196   pmanager = G4PionMinus::PionMinus()->GetProcessManager();
0197
0198   if (pmanager) {
0199     if (decay) pmanager->RemoveProcess(decay);
0200     pmanager->AddProcess(poldecay);
0201     // set ordering for PostStepDoIt and AtRestDoIt
0202     pmanager->SetProcessOrdering(poldecay, idxPostStep);
0203     pmanager->SetProcessOrdering(poldecay, idxAtRest);
0204   }
0205
0206   AddStepMax();
0207 }
0208
0209 /*
0210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0211
0212 void F04PhysicsList::RemoveFromPhysicsList(const G4String& name)
0213 {
0214     G4bool success = false;
0215     for (G4PhysConstVector::iterator p  = physicsVector->begin();
0216                                      p != physicsVector->end(); ++p) {
0217         G4VPhysicsConstructor* e = (*p);
0218         if (e->GetPhysicsName() == name) {
0219            physicsVector->erase(p);
0220            success = true;
0221            break;
0222         }
0223     }
0224     if (!success) {
0225        std::ostringstream message;
0226        message << "PhysicsList::RemoveFromPhysicsList "<< name << "not found";
0227        G4Exception(message.str().c_str());
0228     }
0229 }
0230
0231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0232
0233 void F04PhysicsList::ClearPhysics()
0234 {
0235     for (G4PhysConstVector::iterator p  = physicsVector->begin();
0236                                      p != physicsVector->end(); ++p) {
0237         delete (*p);
0238     }
0239     physicsVector->clear();
0240 }
0241 */
0242
0243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0244
0245 void F04PhysicsList::SetStepMax(G4double step)
0246 {
0247   fMaxChargedStep = step;
0248   fStepMaxProcess->SetStepMax(fMaxChargedStep);
0249 }
0250
0251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0252
0253 F04StepMax* F04PhysicsList::GetStepMaxProcess()
0254 {
0255   return fStepMaxProcess;
0256 }
0257
0258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0259
0260 void F04PhysicsList::AddStepMax()
0261 {
0262   // Step limitation seen as a process
0263
0264   auto particleIterator = GetParticleIterator();
0265   particleIterator->reset();
0266   while ((*particleIterator)()) {
0267     G4ParticleDefinition* particle = particleIterator->value();
0268     G4ProcessManager* pmanager = particle->GetProcessManager();
0269
0270     if (fStepMaxProcess->IsApplicable(*particle) && !particle->IsShortLived()) {
0271       if (pmanager) pmanager->AddDiscreteProcess(fStepMaxProcess);
0272     }
0273   }
0274 }