CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CMSMonopolePhysics.cc
Go to the documentation of this file.
3 
4 #include "G4ParticleDefinition.hh"
5 #include "G4ProcessManager.hh"
6 
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"
15 
16 #include "CLHEP/Units/GlobalSystemOfUnits.h"
17 
19  sim::FieldBuilder * fB_,
20  const edm::ParameterSet & p) :
21  G4VPhysicsConstructor("Monopole Physics"), fieldBuilder(fB_) {
22 
23  verbose = p.getUntrackedParameter<int>("Verbosity",0);
24  magCharge = p.getUntrackedParameter<int>("MonopoleCharge",1);
25  deltaRay = p.getUntrackedParameter<bool>("MonopoleDeltaRay",true);
26  multiSc = p.getUntrackedParameter<bool>("MonopoleMultiScatter",false);
27  transport = p.getUntrackedParameter<bool>("MonopoleTransport",true);
28  if (pdt) {
29  int ii=0;
30  for (HepPDT::ParticleDataTable::const_iterator p=pdt->begin();
31  p != pdt->end(); ++p,++ii) {
32  HepPDT::ParticleData particle = (p->second);
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()));
38  pdgEncodings.push_back(particle.pid());
39  monopoles.push_back(0);
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;
45  }
46  }
47  }
48  if (verbose > 0) G4cout << "CMSMonopolePhysics has " << names.size()
49  << " monopole candidates and delta Ray option "
50  << deltaRay << G4endl;
51 }
52 
54 
56 
57  for (unsigned int ii=0; ii<names.size(); ++ii) {
58  if (!monopoles[ii]) {
59  G4Monopole* mpl = new G4Monopole(names[ii], pdgEncodings[ii],
60  masses[ii], ((pdgEncodings[ii] > 0 ) ? magCharge : -magCharge), elCharges[ii]);;
61  monopoles[ii] = mpl;
62  if (verbose > 0) G4cout << "Create G4Monopole " << names[ii]
63  << " of mass " << masses[ii]/CLHEP::GeV
64  << " GeV, magnetic charge " << ((pdgEncodings[ii] > 0) ? magCharge : -magCharge)
65  << ", electric charge " << elCharges[ii]
66  << " and PDG encoding " << pdgEncodings[ii]
67  << " at " << monopoles[ii] << G4endl;
68  }
69  }
70 }
71 
73  // Add standard EM Processes
74 
75  if (verbose > 0)
76  G4cout << "### CMSMonopolePhysics ConstructProcess()" << G4endl;
77  for (unsigned int ii=0; ii<monopoles.size(); ++ii) {
78  if (monopoles[ii]) {
79  G4Monopole* mpl = monopoles[ii];
80  G4ProcessManager* pmanager = new G4ProcessManager(mpl);
81  mpl->SetProcessManager(pmanager);
82 
83  G4String particleName = mpl->GetParticleName();
84  G4double magn = mpl->MagneticCharge();
85  if (verbose > 1)
86  G4cout << "### CMSMonopolePhysics instantiates for " << particleName
87  << " at " << mpl << " Mass " << mpl->GetPDGMass()/CLHEP::GeV
88  << " GeV Mag " << magn << " Process manager " << pmanager
89  << G4endl;
90 
91  // defined monopole parameters and binning
92 
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;
98 
99  nbin = G4int(std::log10(emax/emin));
100  if (nbin < 1) nbin = 1;
101  nbin *= 10;
102 
103  G4int idx = 1;
104  if (verbose > 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;
108 
109  if (magn == 0.0 || (!transport)) {
110  pmanager->AddProcess( new G4Transportation(verbose), -1, 0, 0);
111  } else {
112  pmanager->AddProcess( new G4MonopoleTransportation(mpl,fieldBuilder,verbose), -1, 0, 0);
113  }
114 
115  if (mpl->GetPDGCharge() != 0.0) {
116  if (multiSc) {
117  G4hMultipleScattering* hmsc = new G4hMultipleScattering();
118  pmanager->AddProcess(hmsc, -1, idx, idx);
119  ++idx;
120  }
121  if (deltaRay) {
122  G4hIonisation* hhioni = new G4hIonisation();
123  hhioni->SetDEDXBinning(nbin);
124  hhioni->SetMinKinEnergy(emin);
125  hhioni->SetMaxKinEnergy(emax);
126  pmanager->AddProcess(hhioni, -1, idx, idx);
127  } else {
128  G4hhIonisation* hhioni = new G4hhIonisation();
129  hhioni->SetDEDXBinning(nbin);
130  hhioni->SetMinKinEnergy(emin);
131  hhioni->SetMaxKinEnergy(emax);
132  pmanager->AddProcess(hhioni, -1, idx, idx);
133  }
134  ++idx;
135  }
136  if(magn != 0.0) {
137  CMSG4mplIonisation* mplioni = new CMSG4mplIonisation(magn);
138  mplioni->SetDEDXBinning(nbin);
139  mplioni->SetMinKinEnergy(emin);
140  mplioni->SetMaxKinEnergy(emax);
141  if (deltaRay) {
142  G4mplIonisationWithDeltaModel* mod = new G4mplIonisationWithDeltaModel(magn,"PAI");
143  mplioni->AddEmModel(0,mod,mod);
144  }
145  pmanager->AddProcess(mplioni, -1, idx, idx);
146  ++idx;
147  }
148  pmanager->AddProcess( new G4StepLimiter(), -1, -1, idx);
149  if (verbose > 1) pmanager->DumpInfo();
150  }
151  }
152 }
T getUntrackedParameter(std::string const &, T const &) const
std::vector< int > pdgEncodings
HepPDT::ParticleDataTable ParticleDataTable
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)
Definition: ecalDccMap.h:4
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
CMSMonopolePhysics(const HepPDT::ParticleDataTable *table_, sim::FieldBuilder *fB_, const edm::ParameterSet &p)