CMSMonopolePhysics.cc

Go to the documentation of this file.

#include "SimG4Core/PhysicsLists/interface/CMSMonopolePhysics.h"
#include "SimG4Core/PhysicsLists/interface/MonopoleTransportation.h"
#include "SimG4Core/PhysicsLists/interface/CMSmplIonisation.h"
#include "SimG4Core/MagneticField/interface/ChordFinderSetter.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "G4ParticleDefinition.hh"
#include "G4ProcessManager.hh"

#include "G4StepLimiter.hh"
#include "G4hMultipleScattering.hh"
#include "G4hIonisation.hh"
#include "G4hhIonisation.hh"
#include "G4mplIonisationModel.hh"

#include "CLHEP/Units/GlobalSystemOfUnits.h"

CMSMonopolePhysics::CMSMonopolePhysics(const HepPDT::ParticleDataTable * pdt,
                       sim::ChordFinderSetter * cfs, 
                       const edm::ParameterSet & p) :
  G4VPhysicsConstructor("Monopole Physics"), chordFinderSetter(cfs) {
  
  verbose   = p.getUntrackedParameter<int>("Verbosity",0);
  magCharge = p.getUntrackedParameter<int>("MonopoleCharge",1);
  deltaRay  = p.getUntrackedParameter<bool>("MonopoleDeltaRay",true);
  multiSc   = p.getUntrackedParameter<bool>("MonopoleMultiScatter",false);
  transport = p.getUntrackedParameter<bool>("MonopoleTransport",true);
  double mass = p.getUntrackedParameter<double>("MonopoleMass",200);
  if (pdt && mass > 0.0) {
    int ii=0;
    for (HepPDT::ParticleDataTable::const_iterator p=pdt->begin(); 
     p != pdt->end(); ++p,++ii) {
      HepPDT::ParticleData particle = (p->second);
      std::string particleName = (particle.name()).substr(0,8);
      if (strcmp(particleName.c_str(),"Monopole") == 0) {
    names.push_back(particle.name());
    masses.push_back(mass*CLHEP::GeV);
    elCharges.push_back((int)(particle.charge()));
    pdgEncodings.push_back(particle.pid());
    monopoles.push_back(nullptr);
    if (verbose > 0) G4cout << "CMSMonopolePhysics: Monopole[" << ii
                << "] " << particleName << " Mass "
                << particle.mass() << " GeV, Magnetic Charge "
                << magCharge << ", Electric Charge "
                << particle.charge() << G4endl;
      } else if(strcmp(particleName.c_str(),"AntiMono") == 0) {
    names.push_back(particle.name());
    masses.push_back(mass*CLHEP::GeV);
    elCharges.push_back((int)(particle.charge()));
    pdgEncodings.push_back(particle.pid());
    monopoles.push_back(nullptr);
    if (verbose > 0) G4cout << "CMSMonopolePhysics: Monopole[" << ii
                << "] " << particleName << " Mass "
                << particle.mass() << " GeV, Magnetic Charge "
                << magCharge << ", Electric Charge "
                << particle.charge() << G4endl; 
      }
    }
  }
  if (verbose > 0) G4cout << "CMSMonopolePhysics has " << names.size()
              << " monopole candidates and delta Ray option " 
              << deltaRay << G4endl;
}

CMSMonopolePhysics::~CMSMonopolePhysics() {}

void CMSMonopolePhysics::ConstructParticle() {
  
  for (unsigned int ii=0; ii<names.size(); ++ii) {
    // monopoles are created once in the master thread
    if (!monopoles[ii]) {
      G4int mc = (pdgEncodings[ii] >= 0 ) ? magCharge : -magCharge;
      Monopole* mpl = new Monopole(names[ii], pdgEncodings[ii], masses[ii],
                                   mc, elCharges[ii]);
      monopoles[ii] = mpl;
      if (verbose > 0) G4cout << "Create Monopole " << names[ii] 
                  << " of mass " << masses[ii]/CLHEP::GeV
                  << " GeV, magnetic charge " << mc 
                  << ", electric charge " << elCharges[ii]
                  << " and PDG encoding " << pdgEncodings[ii]
                  << " at " << monopoles[ii] << G4endl;
    }
  }
}

void CMSMonopolePhysics::ConstructProcess() {
  // Add standard EM Processes
  if (verbose > 0) {
    G4cout << "### CMSMonopolePhysics ConstructProcess()" << G4endl;
  }
  G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();

  for (unsigned int ii=0; ii<monopoles.size(); ++ii) {
    if (monopoles[ii]) {
      Monopole* mpl = monopoles[ii];
      G4ProcessManager *pmanager = mpl->GetProcessManager();
      if(!pmanager) {
        std::ostringstream o;
        o << "Monopole without a Process Manager";
        G4Exception("CMSMonopolePhysics::ConstructProcess()","",
                    FatalException,o.str().c_str());
      }

      G4double magn = mpl->MagneticCharge();
      G4double mass = mpl->GetPDGMass();
      if (verbose > 1) {
    G4cout << "### CMSMonopolePhysics instantiates for " 
           << mpl->GetParticleName()
           << " at " << mpl << " Mass " << mass/CLHEP::GeV 
           << " GeV Mag " << magn  << " Process manager " << pmanager 
           << G4endl;
      }
  
      if (magn != 0.0) {
    pmanager->RemoveProcess(0);
    pmanager->AddProcess(new MonopoleTransportation(mpl,chordFinderSetter,verbose),
                                                        -1, 0, 0);
      }

      if (mpl->GetPDGCharge() != 0.0) {
    if (multiSc) {
      G4hMultipleScattering* hmsc = new G4hMultipleScattering();
      ph->RegisterProcess(hmsc, mpl);
    }
    G4hIonisation* hioni = new G4hIonisation();
    ph->RegisterProcess(hioni, mpl);
      }
      if(magn != 0.0) {
    CMSmplIonisation* mplioni = new CMSmplIonisation(magn);
        if(!deltaRay) {
          G4mplIonisationModel* ion = new G4mplIonisationModel(magn,"PAI");
          ion->SetParticle(mpl);
          mplioni->AddEmModel(0,ion,ion);
    }
    ph->RegisterProcess(mplioni, mpl);
      }
      pmanager->AddDiscreteProcess(new G4StepLimiter());
      if (verbose > 1) { pmanager->DumpInfo(); }
    }
  }
}