CustomPhysics/src/CustomPhysicsList.cc

0001 #include <memory>
0002
0003 #include "SimG4Core/CustomPhysics/interface/CustomPhysicsList.h"
0004 #include "SimG4Core/CustomPhysics/interface/CustomParticleFactory.h"
0005 #include "SimG4Core/CustomPhysics/interface/CustomParticle.h"
0006 #include "SimG4Core/CustomPhysics/interface/DummyChargeFlipProcess.h"
0007 #include "SimG4Core/CustomPhysics/interface/G4ProcessHelper.h"
0008 #include "SimG4Core/CustomPhysics/interface/CustomPDGParser.h"
0009
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/FileInPath.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include "G4hMultipleScattering.hh"
0015 #include "G4hIonisation.hh"
0016 #include "G4ProcessManager.hh"
0017 #include "G4HadronicProcess.hh"
0018
0019 #include "SimG4Core/CustomPhysics/interface/FullModelHadronicProcess.h"
0020 #include "SimG4Core/CustomPhysics/interface/CMSDarkPairProductionProcess.h"
0021 #include "SimG4Core/CustomPhysics/interface/CMSQGSPSIMPBuilder.h"
0022 #include "SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h"
0023
0024 #include "SimG4Core/CustomPhysics/interface/CMSSQLoopProcess.h"
0025 #include "SimG4Core/CustomPhysics/interface/CMSSQLoopProcessDiscr.h"
0026 #include "SimG4Core/CustomPhysics/interface/CMSSQNeutronAnnih.h"
0027 #include "SimG4Core/CustomPhysics/interface/CMSSQInelasticCrossSection.h"
0028
0029 using namespace CLHEP;
0030
0031 G4ThreadLocal std::unique_ptr<G4ProcessHelper> CustomPhysicsList::myHelper;
0032
0033 CustomPhysicsList::CustomPhysicsList(const std::string& name, const edm::ParameterSet& p, bool apinew)
0034     : G4VPhysicsConstructor(name) {
0035   myConfig = p;
0036   if (apinew) {
0037     dfactor = p.getParameter<double>("DarkMPFactor");
0038     fHadronicInteraction = p.getParameter<bool>("RhadronPhysics");
0039   } else {
0040     // this is left for backward compatibility
0041     dfactor = p.getParameter<double>("dark_factor");
0042     fHadronicInteraction = p.getParameter<bool>("rhadronPhysics");
0043   }
0044   edm::FileInPath fp = p.getParameter<edm::FileInPath>("particlesDef");
0045   particleDefFilePath = fp.fullPath();
0046   fParticleFactory = std::make_unique<CustomParticleFactory>();
0047   myHelper.reset(nullptr);
0048
0049   edm::LogVerbatim("SimG4CoreCustomPhysics") << "CustomPhysicsList: Path for custom particle definition file: \n"
0050                                              << particleDefFilePath << "\n"
0051                                              << "      dark_factor= " << dfactor;
0052 }
0053
0054 CustomPhysicsList::~CustomPhysicsList() {}
0055
0056 void CustomPhysicsList::ConstructParticle() {
0057   edm::LogVerbatim("SimG4CoreCustomPhysics") << "===== CustomPhysicsList::ConstructParticle ";
0058   fParticleFactory.get()->loadCustomParticles(particleDefFilePath);
0059 }
0060
0061 void CustomPhysicsList::ConstructProcess() {
0062   edm::LogVerbatim("SimG4CoreCustomPhysics") << "CustomPhysicsList: adding CustomPhysics processes "
0063                                              << "for the list of particles";
0064
0065   G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0066
0067   for (auto particle : fParticleFactory.get()->getCustomParticles()) {
0068     if (particle->GetParticleType() == "simp") {
0069       G4ProcessManager* pmanager = particle->GetProcessManager();
0070       if (pmanager) {
0071         CMSSIMPInelasticProcess* simpInelPr = new CMSSIMPInelasticProcess();
0072         CMSQGSPSIMPBuilder* theQGSPSIMPB = new CMSQGSPSIMPBuilder();
0073         theQGSPSIMPB->Build(simpInelPr);
0074         pmanager->AddDiscreteProcess(simpInelPr);
0075       } else
0076         edm::LogVerbatim("CustomPhysics") << "   No pmanager";
0077     }
0078
0079     CustomParticle* cp = dynamic_cast<CustomParticle*>(particle);
0080     if (cp) {
0081       G4ProcessManager* pmanager = particle->GetProcessManager();
0082       edm::LogVerbatim("SimG4CoreCustomPhysics")
0083           << "CustomPhysicsList: " << particle->GetParticleName() << "  PDGcode= " << particle->GetPDGEncoding()
0084           << "  Mass= " << particle->GetPDGMass() / GeV << " GeV.";
0085       if (pmanager) {
0086         if (particle->GetPDGCharge() != 0.0) {
0087           ph->RegisterProcess(new G4hMultipleScattering, particle);
0088           ph->RegisterProcess(new G4hIonisation, particle);
0089         }
0090         if (cp->GetCloud() && fHadronicInteraction &&
0091             (CustomPDGParser::s_isgluinoHadron(particle->GetPDGEncoding()) ||
0092              (CustomPDGParser::s_isstopHadron(particle->GetPDGEncoding())) ||
0093              (CustomPDGParser::s_issbottomHadron(particle->GetPDGEncoding())))) {
0094           edm::LogVerbatim("SimG4CoreCustomPhysics")
0095               << "CustomPhysicsList: " << particle->GetParticleName()
0096               << " CloudMass= " << cp->GetCloud()->GetPDGMass() / GeV
0097               << " GeV; SpectatorMass= " << cp->GetSpectator()->GetPDGMass() / GeV << " GeV.";
0098
0099           if (!myHelper.get()) {
0100             myHelper = std::make_unique<G4ProcessHelper>(myConfig, fParticleFactory.get());
0101           }
0102           pmanager->AddDiscreteProcess(new FullModelHadronicProcess(myHelper.get()));
0103         }
0104         if (particle->GetParticleType() == "darkpho") {
0105           CMSDarkPairProductionProcess* darkGamma = new CMSDarkPairProductionProcess(dfactor);
0106           pmanager->AddDiscreteProcess(darkGamma);
0107         }
0108         if (particle->GetParticleName() == "anti_sexaq") {
0109           // here the different sexaquark interactions get defined
0110           G4HadronicProcess* sqInelPr = new G4HadronicProcess();
0111           CMSSQNeutronAnnih* sqModel = new CMSSQNeutronAnnih(particle->GetPDGMass() / GeV);
0112           sqInelPr->RegisterMe(sqModel);
0113           CMSSQInelasticCrossSection* sqInelXS = new CMSSQInelasticCrossSection(particle->GetPDGMass() / GeV);
0114           sqInelPr->AddDataSet(sqInelXS);
0115           pmanager->AddDiscreteProcess(sqInelPr);
0116           // add also the looping needed to simulate flat interaction probability
0117           CMSSQLoopProcess* sqLoopPr = new CMSSQLoopProcess();
0118           pmanager->AddContinuousProcess(sqLoopPr);
0119           CMSSQLoopProcessDiscr* sqLoopPrDiscr = new CMSSQLoopProcessDiscr(particle->GetPDGMass() / GeV);
0120           pmanager->AddDiscreteProcess(sqLoopPrDiscr);
0121         } else if (particle->GetParticleName() == "sexaq") {
0122           edm::LogVerbatim("CustomPhysics") << "   No pmanager implemented for sexaq, only for anti_sexaq";
0123         }
0124       }
0125     }
0126   }
0127 }