/* This file is part of MAUS: http://micewww.pp.rl.ac.uk/projects/maus * * MAUS is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * MAUS is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with MAUS. If not, see . * */ #include #include #include #include #include "json/json.h" #include "Geant4/globals.hh" #include "Geant4/G4StepLimiter.hh" #include "Geant4/G4UImanager.hh" #include "Geant4/G4ProcessTable.hh" #include "Geant4/G4ProcessVector.hh" #include "Geant4/G4PhysListFactory.hh" #include "Interface/Squeak.hh" #include "src/common_cpp/Utils/Globals.hh" #include "src/common_cpp/Utils/JsonWrapper.hh" #include "src/common_cpp/Simulation/MAUSPhysicsList.hh" namespace MAUS { // Note muMsc and CoulombScat make an error message from G4.9.2 - but are needed // for G4.9.5 so we leave them in... I could do some logic checking for G4 // version if it is a problem (Rogers) const std::string MAUSPhysicsList::_scatNames[] = {"muBrems", "hBrems", "eBrem", "muPairProd", "hPairProd", "ElectroNuclear", "msc", "muMsc", "CoulombScat"}; const std::string MAUSPhysicsList::_eLossNames[] = {"muBrems", "hBrems", "eBrem", "muPairProd", "hPairProd", "muIoni", "hIoni", "eIoni"}; const int MAUSPhysicsList::_nScatNames = 9; const int MAUSPhysicsList::_nELossNames = 8; MAUSPhysicsList::MAUSPhysicsList(G4VModularPhysicsList* physList) : G4VUserPhysicsList() { _list = physList; } MAUSPhysicsList::~MAUSPhysicsList() {} MAUSPhysicsList* MAUSPhysicsList::GetMAUSPhysicsList() { Json::Value& dc = *Globals::GetInstance()->GetConfigurationCards(); std::string physModel = JsonWrapper::GetProperty (dc, "physics_model", JsonWrapper::stringValue).asString(); G4PhysListFactory fact; if (fact.IsReferencePhysList(physModel)) { MAUSPhysicsList* mpl = new MAUSPhysicsList(fact.GetReferencePhysList(physModel)); try { mpl->Setup(); return mpl; } catch(Squeal squee) { delete mpl; throw squee; } } else { throw(Squeal(Squeal::recoverable, "Failed to recognise physics list model "+physModel, "MAUSPhysicsList::GetMAUSPhysicsList()") ); } } void MAUSPhysicsList::BeginOfReferenceParticleAction() { Squeak::mout(Squeak::debug) << "MAUSPhysicsList::BeginOfReferenceParticleAction() with de model " << _refDEModel << std::endl; SetStochastics(no_scat, _refDEModel, no_hadronic, false); } void MAUSPhysicsList:: BeginOfRunAction() { Squeak::mout(Squeak::debug) << "MAUSPhysicsList::BeginOfRunAction() with" << "\n de model " << _dEModel << "\n msc model " << _msModel << "\n hadronic model " << _hadronicModel << "\n particle decay " << _partDecay << "\n pi 1/2 life " << _piHalfLife << "\n mu 1/2 life " << _muHalfLife << std::endl; SetStochastics(_msModel, _dEModel, _hadronicModel, _partDecay); SetHalfLife(_piHalfLife, _muHalfLife); } void MAUSPhysicsList::SetStochastics(scat scatteringModel, eloss energyLossModel, hadronic hadronicModel, bool decayModel) { SetEnergyLoss(energyLossModel); SetScattering(scatteringModel); SetHadronic(hadronicModel); SetDecay(decayModel); } void MAUSPhysicsList::SetDecay(bool decay) { if (!decay) UIApplyCommand("/process/inactivate Decay"); else UIApplyCommand("/process/activate Decay"); } void MAUSPhysicsList::SetEnergyLoss(eloss eLossModel) { double cutDouble = _productionThreshold; std::stringstream cutStream; std::string elossActive = "activate"; std::string flucActive = "true"; switch (eLossModel) { case energyStraggling: elossActive = "activate"; flucActive = "true"; cutDouble = _productionThreshold; break; case dedx: elossActive = "activate"; flucActive = "false"; cutDouble = 1e12*mm; break; case no_eloss: elossActive = "inactivate"; flucActive = "false"; cutDouble = 1e12*mm; break; } cutStream << cutDouble; std::vector uiCommand; uiCommand.push_back("/run/setCut "+cutStream.str()); uiCommand.push_back("/process/eLoss/fluct "+flucActive); for (int i = 0; i < _nELossNames; i++) uiCommand.push_back("/process/"+elossActive+" "+_eLossNames[i]); for (size_t i = 0; i < uiCommand.size(); i++) { UIApplyCommand(uiCommand[i]); } } void MAUSPhysicsList::SetScattering(scat scatteringModel) { std::string activation; switch (scatteringModel) { case mcs: activation = "/process/activate "; break; case no_scat: activation = "/process/inactivate "; break; } for (int i = 0; i < _nScatNames; i++) { UIApplyCommand(activation+_scatNames[i]); } } void MAUSPhysicsList::SetHadronic(hadronic hadronicModel) { std::string activation; switch (hadronicModel) { case all_hadronic: activation = "/process/activate "; break; case no_hadronic: activation = "/process/inactivate "; break; } G4ProcessVector* pvec = G4ProcessTable::GetProcessTable()-> FindProcesses(fHadronic); for (int i = 0; i < pvec->size(); i++) { std::string pname = (*pvec)[i]->GetProcessName(); UIApplyCommand(activation+pname); } } void MAUSPhysicsList::ConstructParticle() { _list->ConstructParticle(); } void MAUSPhysicsList::ConstructProcess() { _list->ConstructProcess(); SetSpecialProcesses(); } void MAUSPhysicsList::SetSpecialProcesses() { theParticleIterator->reset(); // from G4VUserPhysicsList while ( (*theParticleIterator)() ) { G4ProcessManager* pmanager = theParticleIterator->value()-> GetProcessManager(); pmanager->AddProcess(new G4StepLimiter, -1, -1, 2); } } void MAUSPhysicsList::SetHalfLife(double pionHalfLife, double muonHalfLife) { SetParticleHalfLife("pi+", pionHalfLife); SetParticleHalfLife("pi-", pionHalfLife); SetParticleHalfLife("mu+", muonHalfLife); SetParticleHalfLife("mu-", muonHalfLife); } void MAUSPhysicsList::SetParticleHalfLife(std::string particleName, double halfLife) { if (halfLife <= 0.) return; std::stringstream ss_in; ss_in << halfLife; UIApplyCommand("/particle/select "+particleName); UIApplyCommand("/particle/property/lifetime "+ss_in.str()+" ns"); UIApplyCommand("/particle/property/dump"); } void MAUSPhysicsList::Setup() { Json::Value& dc = *Globals::GetInstance()->GetConfigurationCards(); std::string refPhysicsModel = JsonWrapper::GetProperty(dc, "reference_physics_processes", JsonWrapper::stringValue).asString(); std::string physicsModel = JsonWrapper::GetProperty(dc, "physics_processes", JsonWrapper::stringValue).asString(); if (refPhysicsModel == "none") _refDEModel = no_eloss; if (refPhysicsModel == "mean_energy_loss") _refDEModel = dedx; if (physicsModel == "none") { _msModel = no_scat; _dEModel = no_eloss; _hadronicModel = no_hadronic; } if (physicsModel == "mean_energy_loss") { _msModel = no_scat; _dEModel = dedx; _hadronicModel = no_hadronic; } if (physicsModel == "standard") { _msModel = mcs; _dEModel = energyStraggling; _hadronicModel = all_hadronic; } _partDecay = JsonWrapper::GetProperty(dc, "particle_decay", JsonWrapper::booleanValue).asBool(); _piHalfLife = JsonWrapper::GetProperty(dc, "charged_pion_half_life", JsonWrapper::realValue).asDouble(); _muHalfLife = JsonWrapper::GetProperty(dc, "muon_half_life", JsonWrapper::realValue).asDouble(); _productionThreshold = JsonWrapper::GetProperty(dc, "production_threshold", JsonWrapper::realValue).asDouble(); } void MAUSPhysicsList::UIApplyCommand(std::string command) { Squeak::mout(Squeak::debug) << "Apply G4UI command: " << command << std::endl; G4UImanager::GetUIpointer()->ApplyCommand(command); } }