PhysicsLists/src/CMSMonopolePhysics.cc

0001 #include "SimG4Core/PhysicsLists/interface/CMSMonopolePhysics.h"
0002 #include "SimG4Core/PhysicsLists/interface/MonopoleTransportation.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "G4ParticleDefinition.hh"
0006 #include "G4ProcessManager.hh"
0007
0008 #include "G4StepLimiter.hh"
0009 #include "G4hMultipleScattering.hh"
0010 #include "G4hIonisation.hh"
0011 #include "G4mplIonisation.hh"
0012 #include "G4mplIonisationWithDeltaModel.hh"
0013
0014 #include <CLHEP/Units/SystemOfUnits.h>
0015
0016 CMSMonopolePhysics::CMSMonopolePhysics(const HepPDT::ParticleDataTable* pdt, const edm::ParameterSet& p)
0017     : G4VPhysicsConstructor("Monopole Physics") {
0018   verbose = p.getUntrackedParameter<int>("Verbosity", 0);
0019   magCharge = p.getUntrackedParameter<int>("MonopoleCharge", 1);
0020   deltaRay = p.getUntrackedParameter<bool>("MonopoleDeltaRay", true);
0021   multiSc = p.getUntrackedParameter<bool>("MonopoleMultiScatter", false);
0022   transport = p.getUntrackedParameter<bool>("MonopoleTransport", true);
0023   double mass = p.getUntrackedParameter<double>("MonopoleMass", 200);
0024   if (pdt && mass > 0.0) {
0025     int ii = 0;
0026     for (HepPDT::ParticleDataTable::const_iterator p = pdt->begin(); p != pdt->end(); ++p, ++ii) {
0027       HepPDT::ParticleData particle = (p->second);
0028       std::string particleName = (particle.name()).substr(0, 8);
0029       if (strcmp(particleName.c_str(), "Monopole") == 0) {
0030         names.push_back(particle.name());
0031         masses.push_back(mass * CLHEP::GeV);
0032         elCharges.push_back((int)(particle.charge()));
0033         pdgEncodings.push_back(particle.pid());
0034         monopoles.push_back(nullptr);
0035         if (verbose > 0)
0036           G4cout << "CMSMonopolePhysics: Monopole[" << ii << "] " << particleName << " Mass " << particle.mass()
0037                  << " GeV, Magnetic Charge " << magCharge << ", Electric Charge " << particle.charge() << G4endl;
0038       } else if (strcmp(particleName.c_str(), "AntiMono") == 0) {
0039         names.push_back(particle.name());
0040         masses.push_back(mass * CLHEP::GeV);
0041         elCharges.push_back((int)(particle.charge()));
0042         pdgEncodings.push_back(particle.pid());
0043         monopoles.push_back(nullptr);
0044         if (verbose > 0)
0045           G4cout << "CMSMonopolePhysics: Monopole[" << ii << "] " << particleName << " Mass " << particle.mass()
0046                  << " GeV, Magnetic Charge " << magCharge << ", Electric Charge " << particle.charge() << G4endl;
0047       }
0048     }
0049   }
0050   if (verbose > 0)
0051     G4cout << "CMSMonopolePhysics has " << names.size() << " monopole candidates and delta Ray option " << deltaRay
0052            << G4endl;
0053 }
0054
0055 CMSMonopolePhysics::~CMSMonopolePhysics() {}
0056
0057 void CMSMonopolePhysics::ConstructParticle() {
0058   for (unsigned int ii = 0; ii < names.size(); ++ii) {
0059     // monopoles are created once in the master thread
0060     if (!monopoles[ii]) {
0061       G4int mc = (pdgEncodings[ii] >= 0) ? magCharge : -magCharge;
0062       Monopole* mpl = new Monopole(names[ii], pdgEncodings[ii], masses[ii], mc, elCharges[ii]);
0063       monopoles[ii] = mpl;
0064       if (verbose > 0)
0065         G4cout << "Create Monopole " << names[ii] << " of mass " << masses[ii] / CLHEP::GeV << " GeV, magnetic charge "
0066                << mc << ", electric charge " << elCharges[ii] << " and PDG encoding " << pdgEncodings[ii] << " at "
0067                << monopoles[ii] << G4endl;
0068     }
0069   }
0070 }
0071
0072 void CMSMonopolePhysics::ConstructProcess() {
0073   // Add standard EM Processes
0074   if (verbose > 0) {
0075     G4cout << "### CMSMonopolePhysics ConstructProcess()" << G4endl;
0076   }
0077   G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0078
0079   for (unsigned int ii = 0; ii < monopoles.size(); ++ii) {
0080     if (monopoles[ii]) {
0081       Monopole* mpl = monopoles[ii];
0082       G4ProcessManager* pmanager = mpl->GetProcessManager();
0083       if (!pmanager) {
0084         std::ostringstream o;
0085         o << "Monopole without a Process Manager";
0086         throw edm::Exception(edm::errors::Configuration, o.str().c_str());
0087         return;
0088       }
0089
0090       G4double magn = mpl->MagneticCharge();
0091       G4double mass = mpl->GetPDGMass();
0092       if (verbose > 1) {
0093         G4cout << "### CMSMonopolePhysics instantiates for " << mpl->GetParticleName() << " at " << mpl << " Mass "
0094                << mass / CLHEP::GeV << " GeV Mag " << magn << " Process manager " << pmanager << G4endl;
0095       }
0096
0097       if (magn != 0.0) {
0098         G4int idxt(0);
0099         pmanager->RemoveProcess(idxt);
0100         pmanager->AddProcess(new MonopoleTransportation(mpl, verbose), -1, 0, 0);
0101       }
0102
0103       if (mpl->GetPDGCharge() != 0.0) {
0104         if (multiSc) {
0105           G4hMultipleScattering* hmsc = new G4hMultipleScattering();
0106           ph->RegisterProcess(hmsc, mpl);
0107         }
0108         G4hIonisation* hioni = new G4hIonisation();
0109         ph->RegisterProcess(hioni, mpl);
0110       }
0111       if (magn != 0.0) {
0112         G4mplIonisation* mplioni = new G4mplIonisation(magn);
0113         if (!deltaRay) {
0114           G4mplIonisationWithDeltaModel* ion = new G4mplIonisationWithDeltaModel(magn, "PAI");
0115           ion->SetParticle(mpl);
0116           mplioni->AddEmModel(0, ion, ion);
0117         }
0118         ph->RegisterProcess(mplioni, mpl);
0119       }
0120       pmanager->AddDiscreteProcess(new G4StepLimiter());
0121       if (verbose > 1) {
0122         pmanager->DumpInfo();
0123       }
0124     }
0125   }
0126 }