File indexing completed on 2024-05-10 02:21:26
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/SystemOfUnits.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 using CLHEP::GeV;
0078 LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction Track " << code << " Flag " << actOnTrack
0079 << " (px,py,pz,E) = (" << pxStart / GeV << ", " << pyStart / GeV << ", "
0080 << pzStart / GeV << ", " << eStart / GeV << ")";
0081 }
0082 }
0083 }
0084
0085 void MonopoleSteppingAction::update(const G4Step *aStep) {
0086 if (actOnTrack) {
0087 double eT, pT, pZ, tStep;
0088 G4Track *aTrack = aStep->GetTrack();
0089 G4ThreeVector initialPosition(0, 0, 0);
0090 if (mode) {
0091 tStep = aTrack->GetGlobalTime();
0092 eT = eStart * std::sqrt(dirxStart * dirxStart + diryStart * diryStart);
0093 pT = std::sqrt(pxStart * pxStart + pyStart * pyStart);
0094 pZ = pzStart;
0095 } else {
0096 const G4ThreeVector &dirStep = aTrack->GetMomentumDirection();
0097 double lStep = aTrack->GetStepLength();
0098 double xStep = aTrack->GetPosition().x() - lStep * dirStep.x();
0099 double yStep = aTrack->GetPosition().y() - lStep * dirStep.y();
0100 double zStep = aTrack->GetPosition().z() - lStep * dirStep.z();
0101 double vStep = aTrack->GetVelocity();
0102 initialPosition = G4ThreeVector(xStep, yStep, zStep);
0103 tStep = (vStep > 0. ? (lStep / vStep) : 0.);
0104 double eStep = aTrack->GetTotalEnergy();
0105 eT = eStep * dirStep.perp();
0106 pT = aTrack->GetMomentum().perp();
0107 pZ = aTrack->GetMomentum().z();
0108 }
0109 LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: tStep " << tStep << " eT " << eT << " pT " << pT << " pZ "
0110 << pZ;
0111 double eT1 = eT * cMevToJ;
0112 double pT1 = pT * cMeVToKgMByS;
0113 double pZ1 = pZ * cMeVToKgMByS;
0114 double ts1 = tStep / CLHEP::second;
0115 double eT2 = (eT1 > 0. ? (1. / eT1) : 0.);
0116 double fac0 = magCharge * bZ * cInMByS;
0117 double fac2 = pZ1 * cInMByS * eT2;
0118 double fac1 = fac0 * cInMByS * ts1 * eT2 + fac2;
0119 double fact1 = (eT1 / fac0) * (std::sqrt(1 + fac1 * fac1) - std::sqrt(1 + fac2 * fac2));
0120 double fact2 = (pT1 / (magCharge * bZ) * (asinh(fac1) - asinh(fac2)));
0121 LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: Factor eT = " << eT << " " << eT1 << " " << eT2
0122 << ", pT = " << pT << " " << pT1 << ", pZ = " << pZ << " " << pZ1
0123 << ", time = " << tStep << " " << ts1 << ", Factors " << fac0 << " " << fac1 << " "
0124 << fac2 << ", Final ones " << fact1 << " " << fact2;
0125 G4ThreeVector ez(0, 0, 1), et(std::sqrt(0.5), std::sqrt(0.5), 0);
0126 G4ThreeVector displacement = (fact1 * ez + fact2 * et) * CLHEP::m;
0127
0128 LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: Initial " << initialPosition << " Displacement "
0129 << displacement << " Final " << (initialPosition + displacement);
0130 aTrack->SetPosition(initialPosition + displacement);
0131 }
0132 }
