Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-01 03:54:05

0001 #include "SimG4Core/HelpfulWatchers/interface/MonopoleSteppingAction.h"
0002 
0003 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0004 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0005 #include "SimG4Core/Physics/interface/Monopole.h"
0006 
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 
0011 #include "G4Event.hh"
0012 #include "G4ParticleTable.hh"
0013 #include "G4Run.hh"
0014 #include "G4Track.hh"
0015 
0016 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0017 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0018 
0019 MonopoleSteppingAction::MonopoleSteppingAction(edm::ParameterSet const &p) : actOnTrack(false), bZ(0) {
0020   mode = p.getUntrackedParameter<bool>("ChangeFromFirstStep", true);
0021   edm::LogVerbatim("SimG4CoreWatcher") << "MonopoleSeppingAction set mode for"
0022                                        << " start at first step to " << mode;
0023 }
0024 
0025 MonopoleSteppingAction::~MonopoleSteppingAction() {}
0026 
0027 void MonopoleSteppingAction::registerConsumes(edm::ConsumesCollector cc) {
0028   tok_bFieldH_ = cc.esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>();
0029   edm::LogVerbatim("SimG4CoreWatcher") << "MonopoleSteppingAction::Initialize ESGetToken for MagneticField";
0030 }
0031 
0032 void MonopoleSteppingAction::beginRun(edm::EventSetup const &es) {
0033   const MagneticField *bField = &es.getData(tok_bFieldH_);
0034   const GlobalPoint p(0, 0, 0);
0035   bZ = (bField->inTesla(p)).z();
0036   edm::LogVerbatim("SimG4CoreWatcher") << "Magnetic Field (X): " << (bField->inTesla(p)).x()
0037                                        << " Y: " << (bField->inTesla(p)).y() << " Z: " << bZ;
0038 }
0039 
0040 void MonopoleSteppingAction::update(const BeginOfRun *) {
0041   G4ParticleTable *partTable = G4ParticleTable::GetParticleTable();
0042   magCharge = CLHEP::e_SI / CLHEP::fine_structure_const * 0.5;
0043   for (int ii = 0; ii < partTable->size(); ii++) {
0044     G4ParticleDefinition *particle = partTable->GetParticle(ii);
0045     std::string particleName = (particle->GetParticleName()).substr(0, 8);
0046     if (strcmp(particleName.c_str(), "Monopole") == 0) {
0047       magCharge = CLHEP::e_SI * ((Monopole *)(particle))->MagneticCharge();
0048       pdgCode.push_back(particle->GetPDGEncoding());
0049     }
0050   }
0051   edm::LogVerbatim("SimG4CoreWatcher") << "MonopoleSeppingAction Finds " << pdgCode.size() << " candidates";
0052   for (unsigned int ii = 0; ii < pdgCode.size(); ++ii) {
0053     edm::LogVerbatim("SimG4CoreWatcher") << "PDG Code[" << ii << "] = " << pdgCode[ii];
0054   }
0055   cMevToJ = CLHEP::e_SI / CLHEP::eV;
0056   cMeVToKgMByS = CLHEP::e_SI * CLHEP::meter / (CLHEP::eV * CLHEP::c_light * CLHEP::second);
0057   cInMByS = CLHEP::c_light * CLHEP::second / CLHEP::meter;
0058   edm::LogVerbatim("SimG4CoreWatcher") << "MonopoleSeppingAction Constants:"
0059                                        << " MevToJoules  = " << cMevToJ << ", MevToKgm/s  = " << cMeVToKgMByS
0060                                        << ", c in m/s    = " << cInMByS << ", mag. charge = " << magCharge;
0061 }
0062 
0063 void MonopoleSteppingAction::update(const BeginOfTrack *trk) {
0064   actOnTrack = false;
0065   if (!pdgCode.empty()) {
0066     const G4Track *aTrack = (*trk)();
0067     int code = aTrack->GetDefinition()->GetPDGEncoding();
0068     if (std::count(pdgCode.begin(), pdgCode.end(), code) > 0) {
0069       actOnTrack = true;
0070       eStart = aTrack->GetTotalEnergy();
0071       pxStart = aTrack->GetMomentum().x();
0072       pyStart = aTrack->GetMomentum().y();
0073       pzStart = aTrack->GetMomentum().z();
0074       dirxStart = aTrack->GetMomentumDirection().x();
0075       diryStart = aTrack->GetMomentumDirection().y();
0076       dirzStart = aTrack->GetMomentumDirection().z();
0077       LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction Track " << code << " Flag " << actOnTrack
0078                                    << " (px,py,pz,E) = (" << pxStart / GeV << ", " << pyStart / GeV << ", "
0079                                    << pzStart / GeV << ", " << eStart / GeV << ")";
0080     }
0081   }
0082 }
0083 
0084 void MonopoleSteppingAction::update(const G4Step *aStep) {
0085   if (actOnTrack) {
0086     double eT, pT, pZ, tStep;
0087     G4Track *aTrack = aStep->GetTrack();
0088     G4ThreeVector initialPosition(0, 0, 0);
0089     if (mode) {
0090       tStep = aTrack->GetGlobalTime();
0091       eT = eStart * std::sqrt(dirxStart * dirxStart + diryStart * diryStart);
0092       pT = std::sqrt(pxStart * pxStart + pyStart * pyStart);
0093       pZ = pzStart;
0094     } else {
0095       const G4ThreeVector &dirStep = aTrack->GetMomentumDirection();
0096       double lStep = aTrack->GetStepLength();
0097       double xStep = aTrack->GetPosition().x() - lStep * dirStep.x();
0098       double yStep = aTrack->GetPosition().y() - lStep * dirStep.y();
0099       double zStep = aTrack->GetPosition().z() - lStep * dirStep.z();
0100       double vStep = aTrack->GetVelocity();
0101       initialPosition = G4ThreeVector(xStep, yStep, zStep);
0102       tStep = (vStep > 0. ? (lStep / vStep) : 0.);
0103       double eStep = aTrack->GetTotalEnergy();
0104       eT = eStep * dirStep.perp();
0105       pT = aTrack->GetMomentum().perp();
0106       pZ = aTrack->GetMomentum().z();
0107     }
0108     LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: tStep " << tStep << " eT " << eT << " pT " << pT << " pZ "
0109                                  << pZ;
0110     double eT1 = eT * cMevToJ;
0111     double pT1 = pT * cMeVToKgMByS;
0112     double pZ1 = pZ * cMeVToKgMByS;
0113     double ts1 = tStep / CLHEP::second;
0114     double eT2 = (eT1 > 0. ? (1. / eT1) : 0.);
0115     double fac0 = magCharge * bZ * cInMByS;
0116     double fac2 = pZ1 * cInMByS * eT2;
0117     double fac1 = fac0 * cInMByS * ts1 * eT2 + fac2;
0118     double fact1 = (eT1 / fac0) * (std::sqrt(1 + fac1 * fac1) - std::sqrt(1 + fac2 * fac2));
0119     double fact2 = (pT1 / (magCharge * bZ) * (asinh(fac1) - asinh(fac2)));
0120     LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: Factor eT = " << eT << " " << eT1 << " " << eT2
0121                                  << ", pT = " << pT << " " << pT1 << ", pZ = " << pZ << " " << pZ1
0122                                  << ", time = " << tStep << " " << ts1 << ", Factors " << fac0 << " " << fac1 << " "
0123                                  << fac2 << ", Final ones " << fact1 << " " << fact2;
0124     G4ThreeVector ez(0, 0, 1), et(std::sqrt(0.5), std::sqrt(0.5), 0);
0125     G4ThreeVector displacement = (fact1 * ez + fact2 * et) * CLHEP::m;
0126 
0127     LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: Initial " << initialPosition << " Displacement "
0128                                  << displacement << " Final " << (initialPosition + displacement);
0129     aTrack->SetPosition(initialPosition + displacement);
0130   }
0131 }