4 #include "G4ParticleDefinition.hh"
5 #include "G4ProcessManager.hh"
7 #include "G4StepLimiter.hh"
8 #include "G4Transportation.hh"
9 #include "G4MonopoleTransportation.hh"
10 #include "CMSG4mplIonisation.hh"
11 #include "G4mplIonisationWithDeltaModel.hh"
12 #include "G4hMultipleScattering.hh"
13 #include "G4hIonisation.hh"
14 #include "G4hhIonisation.hh"
16 #include "CLHEP/Units/GlobalSystemOfUnits.h"
21 G4VPhysicsConstructor(
"Monopole Physics"), fieldBuilder(fB_) {
30 for (HepPDT::ParticleDataTable::const_iterator p=pdt->begin();
31 p != pdt->end(); ++
p,++ii) {
33 std::string particleName = (particle.name()).substr(0,8);
34 if (strcmp(particleName.c_str(),
"Monopole") == 0) {
35 names.push_back(particle.name());
36 masses.push_back((particle.mass())*CLHEP::GeV);
37 elCharges.push_back((
int)(particle.charge()));
40 if (
verbose > 0) G4cout <<
"CMSMonopolePhysics: Monopole[" << ii
41 <<
"] " << particleName <<
" Mass "
42 << particle.mass() <<
" GeV, Magnetic Charge "
43 << magCharge <<
", Electric Charge "
44 << particle.charge() << G4endl;
48 if (
verbose > 0) G4cout <<
"CMSMonopolePhysics has " <<
names.size()
49 <<
" monopole candidates and delta Ray option "
50 << deltaRay << G4endl;
57 for (
unsigned int ii=0; ii<
names.size(); ++ii) {
62 if (
verbose > 0) G4cout <<
"Create G4Monopole " <<
names[ii]
63 <<
" of mass " <<
masses[ii]/CLHEP::GeV
76 G4cout <<
"### CMSMonopolePhysics ConstructProcess()" << G4endl;
77 for (
unsigned int ii=0; ii<
monopoles.size(); ++ii) {
80 G4ProcessManager* pmanager =
new G4ProcessManager(mpl);
81 mpl->SetProcessManager(pmanager);
83 G4String particleName = mpl->GetParticleName();
84 G4double magn = mpl->MagneticCharge();
86 G4cout <<
"### CMSMonopolePhysics instantiates for " << particleName
87 <<
" at " << mpl <<
" Mass " << mpl->GetPDGMass()/CLHEP::GeV
88 <<
" GeV Mag " << magn <<
" Process manager " << pmanager
93 G4double emin = mpl->GetPDGMass()/20000.;
94 if(emin < CLHEP::keV) emin = CLHEP::keV;
95 G4double emax = 10.*CLHEP::TeV;
96 G4int nbin = G4int(std::log10(emin/CLHEP::eV));
97 emin =
std::pow(10.,G4double(nbin))*CLHEP::eV;
99 nbin = G4int(std::log10(emax/emin));
100 if (nbin < 1) nbin = 1;
105 G4cout <<
"### Magnetic charge " << magn <<
" and electric charge "
106 << mpl->GetPDGCharge() <<
"\n # of bins in dE/dx table = "
107 << nbin <<
" in the range " << emin <<
":" << emax << G4endl;
110 pmanager->AddProcess(
new G4Transportation(
verbose), -1, 0, 0);
115 if (mpl->GetPDGCharge() != 0.0) {
117 G4hMultipleScattering* hmsc =
new G4hMultipleScattering();
118 pmanager->AddProcess(hmsc, -1, idx, idx);
122 G4hIonisation* hhioni =
new G4hIonisation();
123 hhioni->SetDEDXBinning(nbin);
124 hhioni->SetMinKinEnergy(emin);
125 hhioni->SetMaxKinEnergy(emax);
126 pmanager->AddProcess(hhioni, -1, idx, idx);
128 G4hhIonisation* hhioni =
new G4hhIonisation();
129 hhioni->SetDEDXBinning(nbin);
130 hhioni->SetMinKinEnergy(emin);
131 hhioni->SetMaxKinEnergy(emax);
132 pmanager->AddProcess(hhioni, -1, idx, idx);
137 CMSG4mplIonisation* mplioni =
new CMSG4mplIonisation(magn);
138 mplioni->SetDEDXBinning(nbin);
139 mplioni->SetMinKinEnergy(emin);
140 mplioni->SetMaxKinEnergy(emax);
142 G4mplIonisationWithDeltaModel*
mod =
new G4mplIonisationWithDeltaModel(magn,
"PAI");
143 mplioni->AddEmModel(0,mod,mod);
145 pmanager->AddProcess(mplioni, -1, idx, idx);
148 pmanager->AddProcess(
new G4StepLimiter(), -1, -1, idx);
149 if (
verbose > 1) pmanager->DumpInfo();
T getUntrackedParameter(std::string const &, T const &) const
std::vector< int > pdgEncodings
HepPDT::ParticleDataTable ParticleDataTable
virtual ~CMSMonopolePhysics()
std::vector< G4Monopole * > monopoles
std::vector< int > elCharges
HepPDT::ParticleData ParticleData
std::vector< std::string > names
std::vector< double > masses
sim::FieldBuilder * fieldBuilder
T mod(const T &a, const T &b)
Power< A, B >::type pow(const A &a, const B &b)
CMSMonopolePhysics(const HepPDT::ParticleDataTable *table_, sim::FieldBuilder *fB_, const edm::ParameterSet &p)