diff --git a/SimG4Core/CustomPhysics/BuildFile.xml b/SimG4Core/CustomPhysics/BuildFile.xml
index df15a9aeea4ce..f805ec2132a3f 100644
--- a/SimG4Core/CustomPhysics/BuildFile.xml
+++ b/SimG4Core/CustomPhysics/BuildFile.xml
@@ -1,11 +1,16 @@
+
+
+
+
+
diff --git a/SimG4Core/CustomPhysics/interface/APrimePhysics.h b/SimG4Core/CustomPhysics/interface/APrimePhysics.h
new file mode 100644
index 0000000000000..698d091dc3006
--- /dev/null
+++ b/SimG4Core/CustomPhysics/interface/APrimePhysics.h
@@ -0,0 +1,40 @@
+#ifndef SIMG4CORE_CUSTOMPHYSICS_APRIMEPHYSICS_H
+#define SIMG4CORE_CUSTOMPHYSICS_APRIMEPHYSICS_H
+
+// Geant4
+#include "G4VPhysicsConstructor.hh"
+
+class APrimePhysics : public G4VPhysicsConstructor {
+public:
+ /**
+ * Class constructor.
+ * @param name The name of the physics.
+ */
+ APrimePhysics(double APMass, const G4String& scalefile, const G4double cxBias, const G4String& name = "APrime");
+
+ /**
+ * Class destructor.
+ */
+ ~APrimePhysics() override;
+
+ /**
+ * Construct particles.
+ */
+ void ConstructParticle() override;
+
+ /**
+ * Construct the process.
+ */
+ void ConstructProcess() override;
+
+private:
+ /**
+ * Definition of the APrime particle.
+ */
+ G4ParticleDefinition* aprimeDef_;
+ double apmass;
+ G4String mgfile;
+ G4double biasFactor;
+};
+
+#endif
diff --git a/SimG4Core/CustomPhysics/interface/DBremWatcher.h b/SimG4Core/CustomPhysics/interface/DBremWatcher.h
new file mode 100644
index 0000000000000..3428e872dc85f
--- /dev/null
+++ b/SimG4Core/CustomPhysics/interface/DBremWatcher.h
@@ -0,0 +1,45 @@
+#ifndef SimG4Core_DBremWatcher_H
+#define SimG4Core_DBremWatcher_H
+
+#include "SimG4Core/Notification/interface/Observer.h"
+#include "SimG4Core/Notification/interface/BeginOfTrack.h"
+#include "SimG4Core/Notification/interface/EndOfTrack.h"
+#include "SimG4Core/Notification/interface/BeginOfEvent.h"
+#include "SimG4Core/Notification/interface/BeginOfRun.h"
+#include "SimG4Core/Notification/interface/EndOfEvent.h"
+#include "SimG4Core/Watcher/interface/SimProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "G4ThreeVector.hh"
+
+#include
+#include
+
+class DBremWatcher : public SimProducer,
+ public Observer,
+ public Observer,
+ public Observer,
+ public Observer,
+ public Observer {
+public:
+ DBremWatcher(edm::ParameterSet const &p);
+ ~DBremWatcher() override;
+ void update(const BeginOfTrack *trk) override;
+ void update(const BeginOfEvent *event) override;
+ void update(const EndOfEvent *event) override;
+ void update(const BeginOfRun *run) override;
+ void update(const EndOfTrack *trk) override;
+ void produce(edm::Event &, const edm::EventSetup &) override;
+
+private:
+ std::vector pdgs_;
+ int MotherId;
+ float m_weight;
+ double biasFactor;
+ bool foundbrem;
+ G4ThreeVector aPrimeTraj;
+ G4ThreeVector finaltraj;
+ G4ThreeVector VertexPos;
+ float f_energy;
+};
+
+#endif
diff --git a/SimG4Core/CustomPhysics/interface/G4APrime.h b/SimG4Core/CustomPhysics/interface/G4APrime.h
new file mode 100644
index 0000000000000..13f9feb65ab07
--- /dev/null
+++ b/SimG4Core/CustomPhysics/interface/G4APrime.h
@@ -0,0 +1,41 @@
+/**
+ * @file G4APrime.h
+ * @brief Class creating the A' particle in Geant.
+ * @author Michael Revering, University of Minnesota
+ */
+
+#ifndef G4APrime_h
+#define G4APrime_h
+
+// Geant
+#include "G4ParticleDefinition.hh"
+
+class G4APrime : public G4ParticleDefinition {
+private:
+ static G4APrime* theAPrime;
+
+ G4APrime(const G4String& Name,
+ G4double mass,
+ G4double width,
+ G4double charge,
+ G4int iSpin,
+ G4int iParity,
+ G4int iConjugation,
+ G4int iIsospin,
+ G4int iIsospin3,
+ G4int gParity,
+ const G4String& pType,
+ G4int lepton,
+ G4int baryon,
+ G4int encoding,
+ G4bool stable,
+ G4double lifetime,
+ G4DecayTable* decaytable);
+
+ ~G4APrime() override;
+
+public:
+ static G4APrime* APrime(double apmass = 1000);
+};
+
+#endif
diff --git a/SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h b/SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h
new file mode 100644
index 0000000000000..39f41b4faa817
--- /dev/null
+++ b/SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h
@@ -0,0 +1,40 @@
+/**
+ * @file G4muDarkBremsstrahlung.h
+ * @brief Class providing the Dark Bremsstrahlung process class.
+ * @author Michael Revering, University of Minnesota
+ */
+
+#ifndef G4muDarkBremsstrahlung_h
+#define G4muDarkBremsstrahlung_h
+
+// Geant
+#include "G4VEmProcess.hh"
+
+class G4Material;
+
+class G4muDarkBremsstrahlung : public G4VEmProcess {
+public:
+ G4muDarkBremsstrahlung(const G4String& scalefile, const G4double biasFactor, const G4String& name = "muDBrem");
+
+ ~G4muDarkBremsstrahlung() override;
+
+ G4bool IsApplicable(const G4ParticleDefinition& p) override;
+
+ void PrintInfo() override;
+
+ void SetMethod(std::string method_in);
+
+ G4bool IsEnabled();
+ void SetEnable(bool active);
+ G4muDarkBremsstrahlung& operator=(const G4muDarkBremsstrahlung& right) = delete;
+ G4muDarkBremsstrahlung(const G4muDarkBremsstrahlung&) = delete;
+
+protected:
+ void InitialiseProcess(const G4ParticleDefinition*) override;
+ G4bool isInitialised;
+ const G4String& mgfile;
+ const G4double cxBias;
+ G4bool isEnabled;
+};
+
+#endif
diff --git a/SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlungModel.h b/SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlungModel.h
new file mode 100644
index 0000000000000..27ecbca0fa152
--- /dev/null
+++ b/SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlungModel.h
@@ -0,0 +1,98 @@
+/**
+ * @file G4muDarkBremsstrahlungModel.h
+ * @brief Class provided to simulate the dark brem cross section and interaction.
+ * @author Michael Revering, University of Minnesota
+ */
+
+#ifndef G4muDarkBremsstrahlungModel_h
+#define G4muDarkBremsstrahlungModel_h
+
+// Geant
+#include "G4VEmModel.hh"
+#include "G4Material.hh"
+#include "G4Element.hh"
+#include "G4DataVector.hh"
+#include "G4ParticleChangeForLoss.hh"
+
+// ROOT
+#include "TLorentzVector.h"
+
+struct ParamsForChi {
+ double AA;
+ double ZZ;
+ double MMA;
+ double EE0;
+};
+struct frame {
+ TLorentzVector* fEl;
+ TLorentzVector* cm;
+ G4double E;
+};
+
+class G4Element;
+class G4ParticleChangeForLoss;
+
+class G4muDarkBremsstrahlungModel : public G4VEmModel {
+public:
+ G4muDarkBremsstrahlungModel(const G4String& scalefile,
+ const G4double biasFactor,
+ const G4ParticleDefinition* p = nullptr,
+ const G4String& nam = "eDBrem");
+
+ ~G4muDarkBremsstrahlungModel() override;
+
+ void Initialise(const G4ParticleDefinition*, const G4DataVector&) override;
+
+ G4double ComputeCrossSectionPerAtom(
+ const G4ParticleDefinition*, G4double tkin, G4double Z, G4double, G4double cut, G4double maxE = DBL_MAX) override;
+
+ G4DataVector* ComputePartialSumSigma(const G4Material* material, G4double tkin, G4double cut);
+
+ void SampleSecondaries(std::vector*,
+ const G4MaterialCutsCouple*,
+ const G4DynamicParticle*,
+ G4double tmin,
+ G4double maxEnergy) override;
+
+ void LoadMG();
+
+ void MakePlaceholders();
+
+ void SetMethod(std::string);
+
+ frame GetMadgraphData(double E0);
+ G4muDarkBremsstrahlungModel& operator=(const G4muDarkBremsstrahlungModel& right) = delete;
+ G4muDarkBremsstrahlungModel(const G4muDarkBremsstrahlungModel&) = delete;
+
+protected:
+ const G4Element* SelectRandomAtom(const G4MaterialCutsCouple* couple);
+
+private:
+ void SetParticle(const G4ParticleDefinition* p);
+
+ static G4double chi(double t, void* pp);
+
+ static G4double DsigmaDx(double x, void* pp);
+
+protected:
+ const G4String& mgfile;
+ const G4double cxBias;
+ const G4ParticleDefinition* particle;
+ G4ParticleDefinition* theAPrime;
+ G4ParticleChangeForLoss* fParticleChange;
+ G4double MA;
+ G4bool isMuon;
+
+private:
+ G4double highKinEnergy;
+ G4double lowKinEnergy;
+ G4double probsup;
+ G4bool isInitialised;
+ std::string method;
+ G4bool mg_loaded;
+ std::map > mgdata;
+ std::vector > energies;
+ std::vector partialSumSigma;
+};
+
+#endif
diff --git a/SimG4Core/CustomPhysics/python/DarkBrem_SIM_cfi.py b/SimG4Core/CustomPhysics/python/DarkBrem_SIM_cfi.py
new file mode 100644
index 0000000000000..ca37d54184918
--- /dev/null
+++ b/SimG4Core/CustomPhysics/python/DarkBrem_SIM_cfi.py
@@ -0,0 +1,23 @@
+#To use, add the following to the python configuration:
+#process.load('SimG4Core.CustomPhysics.DarkBrem_SIM_cfi')
+#process.g4SimHits.Physics.type = 'SimG4Core/Physics/CustomPhysics'
+#process.g4SimHits.Physics = cms.PSet(
+#process.g4SimHits.Physics,
+#process.customPhysicsSetup
+#)
+#process.g4SimHits.Watchers = cms.VPSet(cms.PSet(
+# DBremWatcher = cms.PSet(
+# PDGCodes = cms.untracked.vint32([9994]),
+# DBremBiasFactor = process.customPhysicsSetup.DBremBiasFactor
+# ),
+# type = cms.string('DBremWatcher')
+#) )
+
+import FWCore.ParameterSet.Config as cms
+
+customPhysicsSetup = cms.PSet(
+ DBrem = cms.untracked.bool(True),
+ DBremMass = cms.untracked.double(1000.0), #Mass in MeV
+ DBremScaleFile = cms.untracked.string("root://cmseos.fnal.gov//store/user/revering/MuDBrem_Cu_mA1p0.root"),
+ DBremBiasFactor = cms.untracked.double(100)
+)
diff --git a/SimG4Core/CustomPhysics/src/APrimePhysics.cc b/SimG4Core/CustomPhysics/src/APrimePhysics.cc
new file mode 100644
index 0000000000000..575b4aa596b72
--- /dev/null
+++ b/SimG4Core/CustomPhysics/src/APrimePhysics.cc
@@ -0,0 +1,36 @@
+#include "SimG4Core/CustomPhysics/interface/APrimePhysics.h"
+#include "SimG4Core/CustomPhysics/interface/G4APrime.h"
+#include "SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h"
+// Geant 4
+#include "G4Electron.hh"
+#include "G4MuonMinus.hh"
+#include "G4MuonPlus.hh"
+#include "G4ProcessManager.hh"
+#include "G4SystemOfUnits.hh"
+
+APrimePhysics::APrimePhysics(double APMass, const G4String& scalefile, const G4double cxBias, const G4String& name)
+ : G4VPhysicsConstructor(name), aprimeDef_(nullptr) {
+ apmass = APMass;
+ mgfile = scalefile;
+ biasFactor = cxBias;
+}
+
+APrimePhysics::~APrimePhysics() {}
+
+void APrimePhysics::ConstructParticle() {
+ /**
+ * Insert A-prime into the Geant4 particle table.
+ * For now we flag it as stable.
+ */
+ aprimeDef_ = G4APrime::APrime(apmass);
+ //aprimeDef->SetProcessManager(new G4ProcessManager(aprimeDef));
+}
+
+void APrimePhysics::ConstructProcess() {
+ G4ParticleDefinition* muonminus = G4MuonMinus::MuonMinusDefinition();
+ G4ParticleDefinition* muonplus = G4MuonPlus::MuonPlusDefinition();
+ G4ProcessManager* pmplus = muonplus->GetProcessManager();
+ G4ProcessManager* pmminus = muonminus->GetProcessManager();
+ pmplus->AddDiscreteProcess(new G4muDarkBremsstrahlung(mgfile, biasFactor), 6);
+ pmminus->AddDiscreteProcess(new G4muDarkBremsstrahlung(mgfile, biasFactor), 6);
+}
diff --git a/SimG4Core/CustomPhysics/src/CustomPhysics.cc b/SimG4Core/CustomPhysics/src/CustomPhysics.cc
index ac502c89fe49a..51dbfbda36b1a 100644
--- a/SimG4Core/CustomPhysics/src/CustomPhysics.cc
+++ b/SimG4Core/CustomPhysics/src/CustomPhysics.cc
@@ -4,6 +4,7 @@
#include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysics.h"
#include "SimG4Core/PhysicsLists/interface/CMSHadronPhysicsFTFP_BERT.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "SimG4Core/CustomPhysics/interface/APrimePhysics.h"
#include "G4DecayPhysics.hh"
#include "G4EmExtraPhysics.hh"
@@ -18,6 +19,7 @@ CustomPhysics::CustomPhysics(const edm::ParameterSet& p) : PhysicsList(p) {
int ver = p.getUntrackedParameter("Verbosity", 0);
bool tracking = p.getParameter("TrackingCut");
bool ssPhys = p.getUntrackedParameter("ExoticaPhysicsSS", false);
+ bool dbrem = p.getUntrackedParameter("DBrem", false);
double timeLimit = p.getParameter("MaxTrackTime") * ns;
edm::LogInfo("PhysicsList") << "You are using the simulation engine: "
<< "FTFP_BERT_EMM for regular particles \n"
@@ -52,7 +54,11 @@ CustomPhysics::CustomPhysics(const edm::ParameterSet& p) : PhysicsList(p) {
}
// Custom Physics
- if (ssPhys) {
+ if (dbrem) {
+ RegisterPhysics(new APrimePhysics(p.getUntrackedParameter("DBremMass"),
+ p.getUntrackedParameter("DBremScaleFile"),
+ p.getUntrackedParameter("DBremBiasFactor")));
+ } else if (ssPhys) {
RegisterPhysics(new CustomPhysicsListSS("custom", p));
} else {
RegisterPhysics(new CustomPhysicsList("custom", p));
diff --git a/SimG4Core/CustomPhysics/src/DBremWatcher.cc b/SimG4Core/CustomPhysics/src/DBremWatcher.cc
new file mode 100644
index 0000000000000..2d3b14576b6c7
--- /dev/null
+++ b/SimG4Core/CustomPhysics/src/DBremWatcher.cc
@@ -0,0 +1,148 @@
+#include "SimG4Core/CustomPhysics/interface/DBremWatcher.h"
+
+#include "SimG4Core/Notification/interface/BeginOfTrack.h"
+#include "SimG4Core/Notification/interface/BeginOfEvent.h"
+#include "SimG4Core/Notification/interface/BeginOfRun.h"
+#include "SimG4Core/Notification/interface/EndOfEvent.h"
+#include "SimG4Core/Notification/interface/TrackInformation.h"
+#include "SimG4Core/Watcher/interface/SimWatcher.h"
+#include "SimG4Core/Notification/interface/Observer.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "G4Region.hh"
+#include "G4RegionStore.hh"
+#include "G4LogicalVolumeStore.hh"
+#include "G4ProductionCuts.hh"
+#include "G4ProcessTable.hh"
+#include "G4ProcessManager.hh"
+#include "G4MuonMinus.hh"
+#include "SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h"
+#include "G4Track.hh"
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4VProcess.hh"
+#include "G4VParticleChange.hh"
+#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+#include
+
+DBremWatcher::DBremWatcher(edm::ParameterSet const& p) {
+ edm::ParameterSet ps = p.getParameter("DBremWatcher");
+ biasFactor = ps.getUntrackedParameter("DBremBiasFactor", 1);
+ m_weight = 0;
+ foundbrem = false;
+ finaltraj = G4ThreeVector(0, 0, 0);
+ aPrimeTraj = G4ThreeVector(0, 0, 0);
+ VertexPos = G4ThreeVector(0, 0, 0);
+ f_energy = 0;
+ produces("DBremEventWeight");
+ produces("DBremLocationX");
+ produces("DBremLocationY");
+ produces("DBremLocationZ");
+ //produces("DBremMaterial");
+ produces("DBremAngle");
+ produces("DBremInitialEnergy");
+ produces("DBremFinalEnergy");
+ produces("BiasFactor");
+ pdgs_ = ps.getUntrackedParameter>("PDGCodes");
+ edm::LogInfo("DBremWatcher") << "DBremWatcher:: Save Sim Track if PDG code "
+ << "is one from the list of " << pdgs_.size() << " items";
+ for (unsigned int k = 0; k < pdgs_.size(); ++k)
+ edm::LogInfo("DBremWatcher") << "[" << k << "] " << pdgs_[k];
+}
+
+DBremWatcher::~DBremWatcher() {}
+
+void DBremWatcher::update(const BeginOfTrack* trk) {
+ G4Track* theTrack = (G4Track*)((*trk)());
+ TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
+ if (trkInfo) {
+ int pdg = theTrack->GetDefinition()->GetPDGEncoding();
+ G4ThreeVector Vpos = theTrack->GetVertexPosition();
+ const G4VProcess* TrPro = theTrack->GetCreatorProcess();
+ if (TrPro != nullptr) {
+ if ((theTrack->GetCreatorProcess()->GetProcessName()) == "muDBrem") {
+ if (std::find(pdgs_.begin(), pdgs_.end(), pdg) == pdgs_.end()) {
+ //Found the deflected muon track
+ trkInfo->storeTrack(true);
+ if (!theTrack->IsGoodForTracking()) {
+ theTrack->SetGoodForTrackingFlag(true);
+ }
+ f_energy = theTrack->GetTotalEnergy();
+ foundbrem = true;
+ finaltraj = theTrack->GetMomentum();
+ } else {
+ m_weight = theTrack->GetWeight();
+ }
+ }
+ }
+ if (std::find(pdgs_.begin(), pdgs_.end(), pdg) != pdgs_.end()) {
+ //Found an A'
+ trkInfo->storeTrack(true);
+ VertexPos = Vpos;
+ aPrimeTraj = theTrack->GetMomentum();
+ LogDebug("DBremWatcher") << "Save SimTrack the Track " << theTrack->GetTrackID() << " Type "
+ << theTrack->GetDefinition()->GetParticleName() << " Momentum "
+ << theTrack->GetMomentum() / MeV << " MeV/c";
+ }
+ }
+}
+
+void DBremWatcher::update(const BeginOfRun* run) {}
+
+void DBremWatcher::update(const BeginOfEvent* event) {
+ G4String pname = "muDBrem";
+ G4ProcessTable* ptable = G4ProcessTable::GetProcessTable();
+ G4bool state = true;
+ ptable->SetProcessActivation(pname, state);
+ foundbrem = false;
+}
+
+void DBremWatcher::update(const EndOfEvent* event) {}
+
+void DBremWatcher::update(const EndOfTrack* trk) {
+ G4Track* theTrack = (G4Track*)((*trk)());
+ TrackInformation* trkInfo = (TrackInformation*)(theTrack->GetUserInformation());
+ const G4VProcess* TrPro = theTrack->GetCreatorProcess();
+ if (trkInfo && TrPro != nullptr) {
+ int pdg = theTrack->GetDefinition()->GetPDGEncoding();
+
+ if (std::find(pdgs_.begin(), pdgs_.end(), pdg) == pdgs_.end() &&
+ (theTrack->GetCreatorProcess()->GetProcessName()) == "muDBrem") {
+ trkInfo->storeTrack(true);
+ }
+ }
+}
+
+void DBremWatcher::produce(edm::Event& fEvent, const edm::EventSetup&) {
+ if (foundbrem) {
+ std::unique_ptr weight = std::make_unique(m_weight);
+ fEvent.put(std::move(weight), "DBremEventWeight");
+ std::unique_ptr vtxposx = std::make_unique(VertexPos.x());
+ std::unique_ptr vtxposy = std::make_unique(VertexPos.y());
+ std::unique_ptr vtxposz = std::make_unique(VertexPos.z());
+ fEvent.put(std::move(vtxposx), "DBremLocationX");
+ fEvent.put(std::move(vtxposy), "DBremLocationY");
+ fEvent.put(std::move(vtxposz), "DBremLocationZ");
+ std::unique_ptr finalE = std::make_unique(f_energy / GeV);
+ fEvent.put(std::move(finalE), "DBremFinalEnergy");
+ float deflectionAngle = -1;
+ float initialEnergy = sqrt(pow(finaltraj.x() + aPrimeTraj.x(), 2) + pow(finaltraj.y() + aPrimeTraj.y(), 2) +
+ pow(finaltraj.z() + aPrimeTraj.z(), 2) + pow(0.1056 * GeV, 2));
+ G4ThreeVector mother(
+ finaltraj.x() + aPrimeTraj.x(), finaltraj.y() + aPrimeTraj.y(), finaltraj.z() + aPrimeTraj.z());
+ deflectionAngle = mother.angle(finaltraj);
+ std::unique_ptr dAngle = std::make_unique(deflectionAngle);
+ std::unique_ptr initialE = std::make_unique(initialEnergy / GeV);
+ fEvent.put(std::move(dAngle), "DBremAngle");
+ fEvent.put(std::move(initialE), "DBremInitialEnergy");
+ std::unique_ptr bias = std::make_unique(biasFactor);
+ fEvent.put(std::move(bias), "BiasFactor");
+ } else {
+ std::unique_ptr weight = std::make_unique(0.);
+ fEvent.put(std::move(weight), "DBremEventWeight");
+ }
+}
+
+#include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
+
+DEFINE_SIMWATCHER(DBremWatcher);
diff --git a/SimG4Core/CustomPhysics/src/G4APrime.cc b/SimG4Core/CustomPhysics/src/G4APrime.cc
new file mode 100644
index 0000000000000..34064773e5fe3
--- /dev/null
+++ b/SimG4Core/CustomPhysics/src/G4APrime.cc
@@ -0,0 +1,82 @@
+#include "SimG4Core/CustomPhysics/interface/G4APrime.h"
+#include "G4SystemOfUnits.hh"
+
+G4APrime* G4APrime::theAPrime = nullptr;
+
+G4APrime::G4APrime(const G4String& aName,
+ G4double mass,
+ G4double width,
+ G4double charge,
+ G4int iSpin,
+ G4int iParity,
+ G4int iConjugation,
+ G4int iIsospin,
+ G4int iIsospin3,
+ G4int gParity,
+ const G4String& pType,
+ G4int lepton,
+ G4int baryon,
+ G4int encoding,
+ G4bool stable,
+ G4double lifetime,
+ G4DecayTable* decaytable)
+ : G4ParticleDefinition(aName,
+ mass,
+ width,
+ charge,
+ iSpin,
+ iParity,
+ iConjugation,
+ iIsospin,
+ iIsospin3,
+ gParity,
+ pType,
+ lepton,
+ baryon,
+ encoding,
+ stable,
+ lifetime,
+ decaytable) {}
+
+G4APrime::~G4APrime() {}
+
+G4APrime* G4APrime::APrime(double apmass) {
+ if (!theAPrime) {
+ const G4String& name = "A^1";
+ G4double mass = apmass * MeV;
+ G4double width = 0.;
+ G4double charge = 0;
+ G4int iSpin = 0;
+ G4int iParity = 0;
+ G4int iConjugation = 0;
+ G4int iIsospin = 0;
+ G4int iIsospin3 = 0;
+ G4int gParity = 0;
+ const G4String& pType = "APrime";
+ G4int lepton = 0;
+ G4int baryon = 0;
+ G4int encoding = 9994;
+ G4bool stable = true;
+ G4double lifetime = -1;
+ G4DecayTable* decaytable = nullptr;
+
+ theAPrime = new G4APrime(name,
+ mass,
+ width,
+ charge,
+ iSpin,
+ iParity,
+ iConjugation,
+ iIsospin,
+ iIsospin3,
+ gParity,
+ pType,
+ lepton,
+ baryon,
+ encoding,
+ stable,
+ lifetime,
+ decaytable);
+ }
+ return theAPrime;
+}
diff --git a/SimG4Core/CustomPhysics/src/G4muDarkBremsstrahlung.cc b/SimG4Core/CustomPhysics/src/G4muDarkBremsstrahlung.cc
new file mode 100644
index 0000000000000..5d35d9e0e1b5f
--- /dev/null
+++ b/SimG4Core/CustomPhysics/src/G4muDarkBremsstrahlung.cc
@@ -0,0 +1,48 @@
+#include "SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h"
+#include "SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlungModel.h"
+#include "SimG4Core/CustomPhysics/interface/G4APrime.h"
+
+//Geant 4
+#include "G4MuonMinus.hh"
+#include "G4MuonPlus.hh"
+#include "G4LossTableManager.hh"
+
+using namespace std;
+
+G4muDarkBremsstrahlung::G4muDarkBremsstrahlung(const G4String& scalefile,
+ const G4double biasFactor,
+ const G4String& name)
+ : G4VEmProcess(name), isInitialised(false), mgfile(scalefile), cxBias(biasFactor) {
+ G4int subtype = 40;
+ SetProcessSubType(subtype);
+ SetSecondaryParticle(G4APrime::APrime());
+}
+
+G4muDarkBremsstrahlung::~G4muDarkBremsstrahlung() {}
+
+G4bool G4muDarkBremsstrahlung::IsApplicable(const G4ParticleDefinition& p) {
+ return (&p == G4MuonPlus::MuonPlus() || &p == G4MuonMinus::MuonMinus());
+}
+
+void G4muDarkBremsstrahlung::InitialiseProcess(const G4ParticleDefinition*) {
+ if (!isInitialised) {
+ AddEmModel(0, new G4muDarkBremsstrahlungModel(mgfile, cxBias));
+
+ isInitialised = true;
+ isEnabled = true;
+ }
+}
+
+void G4muDarkBremsstrahlung::PrintInfo() {}
+
+void G4muDarkBremsstrahlung::SetMethod(std::string method_in) {
+ ((G4muDarkBremsstrahlungModel*)EmModel(1))->SetMethod(method_in);
+ return;
+}
+
+G4bool G4muDarkBremsstrahlung::IsEnabled() { return isEnabled; }
+
+void G4muDarkBremsstrahlung::SetEnable(bool state) {
+ isEnabled = state;
+ return;
+}
diff --git a/SimG4Core/CustomPhysics/src/G4muDarkBremsstrahlungModel.cc b/SimG4Core/CustomPhysics/src/G4muDarkBremsstrahlungModel.cc
new file mode 100644
index 0000000000000..a1b6a53945020
--- /dev/null
+++ b/SimG4Core/CustomPhysics/src/G4muDarkBremsstrahlungModel.cc
@@ -0,0 +1,454 @@
+#include "SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlungModel.h"
+#include "SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h"
+#include "SimG4Core/CustomPhysics/interface/G4APrime.h"
+
+//Geant 4
+#include "G4ProcessTable.hh"
+#include "G4ProcTblElement.hh"
+#include "G4MuonMinus.hh"
+#include "G4MuonPlus.hh"
+#include "G4ProductionCutsTable.hh"
+#include "G4SystemOfUnits.hh"
+//Root
+#include "TFile.h"
+#include "TTree.h"
+//gsl
+#include
+#include
+
+using namespace std;
+
+G4muDarkBremsstrahlungModel::G4muDarkBremsstrahlungModel(const G4String& scalefile,
+ const G4double biasFactor,
+ const G4ParticleDefinition* p,
+ const G4String& nam)
+ : G4VEmModel(nam),
+ mgfile(scalefile),
+ cxBias(biasFactor),
+ particle(nullptr),
+ isMuon(true),
+ probsup(1.0),
+ isInitialised(false),
+ method("forward_only"),
+ mg_loaded(false) {
+ if (p) {
+ SetParticle(p);
+ } //Verify that the particle is a muon.
+ theAPrime = G4APrime::APrime();
+ MA = G4APrime::APrime()->GetPDGMass() / CLHEP::GeV; //Get the A' mass.
+
+ highKinEnergy = HighEnergyLimit();
+ lowKinEnergy = LowEnergyLimit();
+ fParticleChange = nullptr;
+}
+
+G4muDarkBremsstrahlungModel::~G4muDarkBremsstrahlungModel() {
+ size_t n = partialSumSigma.size();
+ if (n > 0) {
+ for (size_t i = 0; i < n; i++) {
+ delete partialSumSigma[i];
+ }
+ }
+}
+
+void G4muDarkBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p) {
+ particle = p;
+
+ if ((p == G4MuonPlus::MuonPlus()) || (p == G4MuonMinus::MuonMinus())) {
+ isMuon = true;
+ } else {
+ isMuon = false;
+ }
+}
+
+void G4muDarkBremsstrahlungModel::Initialise(const G4ParticleDefinition* p, const G4DataVector& cuts) {
+ if (p) {
+ SetParticle(p);
+ }
+
+ highKinEnergy = HighEnergyLimit();
+ lowKinEnergy = LowEnergyLimit();
+ const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
+
+ if (theCoupleTable) {
+ G4int numOfCouples = theCoupleTable->GetTableSize();
+ G4int nn = partialSumSigma.size();
+ G4int nc = cuts.size();
+ if (nn > 0) {
+ for (G4int ii = 0; ii < nn; ii++) {
+ G4DataVector* a = partialSumSigma[ii];
+ if (a)
+ delete a;
+ }
+ partialSumSigma.clear();
+ }
+ if (numOfCouples > 0) {
+ for (G4int i = 0; i < numOfCouples; i++) {
+ G4double cute = DBL_MAX;
+ if (i < nc)
+ cute = cuts[i];
+ const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
+ const G4Material* material = couple->GetMaterial();
+ G4DataVector* dv = ComputePartialSumSigma(material, 0.5 * highKinEnergy, std::min(cute, 0.25 * highKinEnergy));
+ partialSumSigma.push_back(dv);
+ }
+ }
+ }
+
+ if (isInitialised)
+ return;
+ fParticleChange = GetParticleChangeForLoss();
+ isInitialised = true;
+}
+
+void G4muDarkBremsstrahlungModel::LoadMG()
+//Parses a Madgraph root file to extract the kinetic energy fraction and pt of the outgoing electron in each event. Loads the two numbers from every event into a map of vectors of pairs (mgdata). Map is keyed by energy, vector pairs are energy fraction + pt. Also creates an list of energies and placeholders (energies), so that different energies can be looped separately.
+{
+ TFile* f = TFile::Open(mgfile.c_str());
+ TTree* tree = (TTree*)f->Get("Events");
+ TLorentzVector* mvec = new TLorentzVector();
+ TLorentzVector* avec = new TLorentzVector();
+ TLorentzVector* nvec = new TLorentzVector();
+ tree->SetBranchAddress("IncidentParticle", &mvec);
+ tree->SetBranchAddress("APrime", &avec);
+ tree->SetBranchAddress("Nucleus", &nvec);
+ int entries = tree->GetEntries();
+ int start = int(G4UniformRand() * entries);
+ for (int i = start; i < int(start + entries / 10.); i++) {
+ if (i < entries) {
+ tree->GetEntry(i);
+ } else {
+ tree->GetEntry(i - entries);
+ }
+ frame evnt;
+ evnt.fEl = (TLorentzVector*)mvec->Clone();
+ evnt.cm = (TLorentzVector*)avec->Clone();
+ *evnt.cm = *evnt.cm + *evnt.fEl;
+ TLorentzVector* ebeam = (TLorentzVector*)nvec->Clone();
+ *ebeam = *ebeam + *evnt.cm;
+ evnt.E = round(ebeam->Z() * 10.) / 10.;
+ if (mgdata.count(evnt.E) == 0) {
+ mgdata[evnt.E];
+ }
+ mgdata[evnt.E].push_back(evnt);
+ }
+ f->Close();
+}
+
+void G4muDarkBremsstrahlungModel::MakePlaceholders() {
+ //Need to do this to set up random sampling of mg distributions
+ for (const auto& iter : mgdata) {
+ energies.push_back(std::make_pair(iter.first, iter.second.size()));
+ }
+
+ for (uint64_t i = 0; i < energies.size(); i++) {
+ energies[i].second = int(G4UniformRand() * mgdata[energies[i].first].size());
+ }
+}
+
+frame G4muDarkBremsstrahlungModel::GetMadgraphData(double E0)
+//Gets the energy fraction and Pt from the imported LHE data. E0 should be in GeV, returns the total energy and Pt in GeV. Scales from the closest imported beam energy above the given value (scales down to avoid biasing issues).
+{
+ double samplingE = energies[0].first;
+ frame cmdata;
+ uint64_t i = 0;
+ bool pass = false;
+ //G4double Mel = 5.1E-04;
+ //Cycle through imported beam energies until the closest one above is found, or the max is reached.
+ while (!pass) {
+ i++;
+ samplingE = energies[i].first;
+ if ((E0 <= samplingE) || (i >= energies.size())) {
+ pass = true;
+ }
+ }
+
+ if (i == energies.size()) {
+ i--;
+ } //Decrement if the loop hit the maximum, to prevent a segfault.
+ //energies is a vector of pairs. energies[i].first holds the imported beam energy,
+ //energies[i].second holds the place as we loop through different energy sampling files.
+ //Need to loop around if we hit the end, when the size of mgdata is smaller than our index
+ //on energies[i].second.
+ if (energies[i].second >= double(mgdata[energies[i].first].size())) {
+ energies[i].second = 0;
+ }
+
+ //Get the lorentz vectors from the index given by the placeholder.
+ cmdata = mgdata[energies[i].first].at(energies[i].second);
+
+ //Increment the placeholder.
+ energies[i].second = energies[i].second + 1;
+
+ return cmdata;
+}
+
+G4double G4muDarkBremsstrahlungModel::DsigmaDx(double x, void* pp) {
+ ParamsForChi* params = (ParamsForChi*)pp;
+
+ //G4double Mel = 5.1E-04;
+ G4double Mmu = 0.106;
+
+ G4double beta = sqrt(1 - (params->MMA) * (params->MMA) / (params->EE0) / (params->EE0));
+ G4double num = 1. - x + x * x / 3.;
+ G4double denom = (params->MMA) * (params->MMA) * (1. - x) / x + Mmu * Mmu * x;
+ G4double DsDx = beta * num / denom;
+
+ return DsDx;
+}
+
+G4double G4muDarkBremsstrahlungModel::chi(double t, void* pp) {
+ ParamsForChi* params = (ParamsForChi*)pp;
+ /* Reminder II:
+ * params->AA;
+ * params->ZZ;
+ * params->MMA;
+ * params->EE0;
+ */
+ G4double Mel = 5.1E-04;
+ G4double MUp = 2.79;
+ G4double Mpr = 0.938;
+
+ G4double d = 0.164 / pow((params->AA), 2. / 3.);
+ G4double ap = 773.0 / Mel / pow((params->ZZ), 2. / 3.);
+ G4double a = 111.0 / Mel / pow((params->ZZ), 1. / 3.);
+ G4double G2el = (params->ZZ) * (params->ZZ) * a * a * a * a * t * t / (1.0 + a * a * t) / (1.0 + a * a * t) /
+ (1.0 + t / d) / (1.0 + t / d);
+ G4double G2in = (params->ZZ) * ap * ap * ap * ap * t * t / (1.0 + ap * ap * t) / (1.0 + ap * ap * t) /
+ (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) /
+ (1.0 + t / 0.71) / (1.0 + t / 0.71) / (1.0 + t / 0.71) *
+ (1.0 + t * (MUp * MUp - 1.0) / 4.0 / Mpr / Mpr) * (1.0 + t * (MUp * MUp - 1.0) / 4.0 / Mpr / Mpr);
+ G4double G2 = G2el + G2in;
+ G4double ttmin = (params->MMA) * (params->MMA) * (params->MMA) * (params->MMA) / 4.0 / (params->EE0) / (params->EE0);
+ // G4double ttmin = lowerLimit(x,theta,p);
+ G4double Under = G2 * (t - ttmin) / t / t;
+
+ return Under;
+}
+
+G4double G4muDarkBremsstrahlungModel::ComputeCrossSectionPerAtom(
+ const G4ParticleDefinition*, G4double E0, G4double Z, G4double A, G4double cut, G4double)
+// Calculates the cross section per atom in GEANT4 internal units. Uses WW approximation to find the total cross section, performing numerical integrals over x and theta.
+{
+ G4double cross = 0.0;
+ if (E0 < keV || E0 < cut) {
+ return cross;
+ }
+
+ E0 = E0 / CLHEP::GeV; //Change energy to GeV.
+ //G4double Mel = 5.1E-04;
+ G4double Mmu = 0.106;
+ if (E0 < 2. * MA)
+ return 0.;
+
+ //begin: chi-formfactor calculation
+ gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
+
+ G4double result, error;
+ G4double tmin = MA * MA * MA * MA / (4. * E0 * E0);
+ G4double tmax = MA * MA;
+
+ gsl_function F;
+ ParamsForChi alpha = {1.0, 1.0, 1.0, 1.0};
+ F.function = &G4muDarkBremsstrahlungModel::chi;
+ F.params = α
+ alpha.AA = A;
+ alpha.ZZ = Z;
+ alpha.MMA = MA;
+ alpha.EE0 = E0;
+
+ //Integrate over chi.
+ gsl_integration_qags(&F, tmin, tmax, 0, 1e-7, 1000, w, &result, &error);
+
+ G4double ChiRes = result;
+ gsl_integration_workspace_free(w);
+
+ //Integrate over x. Can use log approximation instead, which falls off at high A' mass.
+ gsl_integration_workspace* dxspace = gsl_integration_workspace_alloc(1000);
+ gsl_function G;
+ G.function = &DsigmaDx;
+ G.params = α
+ G4double xmin = 0;
+ G4double xmax = 1;
+ if ((Mmu / E0) > (MA / E0))
+ xmax = 1 - Mmu / E0;
+ else
+ xmax = 1 - MA / E0;
+ G4double res, err;
+
+ gsl_integration_qags(&G, xmin, xmax, 0, 1e-7, 1000, dxspace, &res, &err);
+
+ G4double DsDx = res;
+ gsl_integration_workspace_free(dxspace);
+
+ G4double GeVtoPb = 3.894E08;
+ G4double alphaEW = 1.0 / 137.0;
+ G4double epsilBench = 1;
+
+ cross = GeVtoPb * 4. * alphaEW * alphaEW * alphaEW * epsilBench * epsilBench * ChiRes * DsDx * CLHEP::picobarn;
+ if (cross < 0.) {
+ cross = 0.;
+ }
+ E0 = E0 * CLHEP::GeV;
+
+ cross = cross * cxBias;
+ return cross;
+}
+
+G4DataVector* G4muDarkBremsstrahlungModel::ComputePartialSumSigma(const G4Material* material,
+ G4double kineticEnergy,
+ G4double cut)
+
+// Build the table of cross section per element.
+//The table is built for MATERIALS.
+// This table is used by DoIt to select randomly an element in the material.
+{
+ G4int nElements = material->GetNumberOfElements();
+ const G4ElementVector* theElementVector = material->GetElementVector();
+ const G4double* theAtomNumDensityVector = material->GetAtomicNumDensityVector();
+ G4DataVector* dv = new G4DataVector();
+ G4double cross = 0.0;
+
+ for (G4int i = 0; i < nElements; i++) {
+ cross += theAtomNumDensityVector[i] *
+ ComputeCrossSectionPerAtom(
+ particle, kineticEnergy, (*theElementVector)[i]->GetZ(), (*theElementVector)[i]->GetA(), cut);
+ dv->push_back(cross);
+ }
+ return dv;
+}
+
+void G4muDarkBremsstrahlungModel::SampleSecondaries(std::vector* vdp,
+ const G4MaterialCutsCouple* couple,
+ const G4DynamicParticle* dp,
+ G4double tmin,
+ G4double maxEnergy)
+//Simulates the emission of a dark photon + electron. Gets an energy fraction and Pt from madgraph files. Scales the energy so that the fraction of kinectic energy is constant, keeps the Pt constant. If the Pt is larger than the new energy, that event is skipped, and a new one is taken from the file.
+{
+ //Deactivate the process after one dark brem. Needs to be reactivated in the end of event action. If this is in the stepping action instead, more than one brem can occur within each step.
+ G4bool state = false;
+ G4String pname = "muDBrem";
+ G4ProcessTable* ptable = G4ProcessTable::GetProcessTable();
+ ptable->SetProcessActivation(pname, state);
+ if (mg_loaded == false) {
+ LoadMG();
+ MakePlaceholders(); //Setup the placeholder offsets for getting data.
+ mg_loaded = true;
+ }
+ G4double E0 = dp->GetTotalEnergy();
+ G4double tmax = min(maxEnergy, E0);
+ if (tmin >= tmax) {
+ return;
+ } // limits of the energy sampling
+ G4double Mmu = 0.106;
+ E0 = E0 / CLHEP::GeV; //Convert the energy to GeV, the units used in the sampling files.
+
+ frame data = GetMadgraphData(E0);
+ double EAcc, Pt, P, PhiAcc;
+ if (method == "forward_only") {
+ EAcc = (data.fEl->E() - Mmu) / (data.E - Mmu - MA) * (E0 - Mmu - MA);
+ EAcc = Mmu + EAcc;
+ Pt = data.fEl->Pt();
+ P = sqrt(EAcc * EAcc - Mmu * Mmu);
+ PhiAcc = data.fEl->Phi();
+ int i = 0;
+ while (Pt * Pt + Mmu * Mmu > EAcc * EAcc) //Skip events until the Pt is less than the energy.
+ {
+ i++;
+ data = GetMadgraphData(E0);
+ EAcc = (data.fEl->E() - Mmu) / (data.E - Mmu - MA) * (E0 - Mmu - MA);
+ EAcc = Mmu + EAcc;
+ Pt = data.fEl->Pt();
+ P = sqrt(EAcc * EAcc - Mmu * Mmu);
+ PhiAcc = data.fEl->Phi();
+
+ if (i > 10000) {
+ break;
+ }
+ }
+ } else if (method == "cm_scaling") {
+ TLorentzVector* el = new TLorentzVector(data.fEl->X(), data.fEl->Y(), data.fEl->Z(), data.fEl->E());
+ double ediff = data.E - E0;
+ TLorentzVector* newcm = new TLorentzVector(data.cm->X(), data.cm->Y(), data.cm->Z() - ediff, data.cm->E() - ediff);
+ el->Boost(-1. * data.cm->BoostVector());
+ el->Boost(newcm->BoostVector());
+ double newE = (data.fEl->E() - Mmu) / (data.E - Mmu - MA) * (E0 - Mmu - MA);
+ el->SetE(newE);
+ EAcc = el->E();
+ Pt = el->Pt();
+ P = el->P();
+ PhiAcc = el->Phi();
+ } else {
+ EAcc = E0;
+ P = dp->GetTotalMomentum();
+ Pt = sqrt(dp->Get4Momentum().px() * dp->Get4Momentum().px() + dp->Get4Momentum().py() * dp->Get4Momentum().py());
+ PhiAcc = dp->Get4Momentum().phi();
+ }
+
+ EAcc = EAcc * CLHEP::GeV; //Change the energy back to MeV, the internal GEANT unit.
+
+ G4double muon_mass_mev = 105.6;
+ G4double momentum = sqrt(EAcc * EAcc - muon_mass_mev * muon_mass_mev); //Muon momentum in MeV.
+ G4ThreeVector newDirection;
+ double ThetaAcc = std::asin(Pt / P);
+ newDirection.set(std::sin(ThetaAcc) * std::cos(PhiAcc), std::sin(ThetaAcc) * std::sin(PhiAcc), std::cos(ThetaAcc));
+ newDirection.rotateUz(dp->GetMomentumDirection());
+ newDirection.setMag(momentum);
+ // create g4dynamicparticle object for the dark photon.
+ G4ThreeVector direction = (dp->GetMomentum() - newDirection);
+ G4DynamicParticle* dphoton = new G4DynamicParticle(theAPrime, direction);
+ vdp->push_back(dphoton);
+
+ // energy of primary
+ G4double finalKE = EAcc - muon_mass_mev;
+
+ // stop tracking and create new secondary instead of primary
+ bool makeSecondary = true;
+ if (makeSecondary) {
+ fParticleChange->ProposeTrackStatus(fStopAndKill);
+ fParticleChange->SetProposedKineticEnergy(0.0);
+ G4DynamicParticle* mu =
+ new G4DynamicParticle(const_cast(particle), newDirection.unit(), finalKE);
+ vdp->push_back(mu);
+ // continue tracking
+ } else {
+ fParticleChange->SetProposedMomentumDirection(newDirection.unit());
+ fParticleChange->SetProposedKineticEnergy(finalKE);
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+const G4Element* G4muDarkBremsstrahlungModel::SelectRandomAtom(const G4MaterialCutsCouple* couple) {
+ // select randomly 1 element within the material
+
+ const G4Material* material = couple->GetMaterial();
+ G4int nElements = material->GetNumberOfElements();
+ const G4ElementVector* theElementVector = material->GetElementVector();
+
+ const G4Element* elm = nullptr;
+
+ if (1 < nElements) {
+ --nElements;
+ G4DataVector* dv = partialSumSigma[couple->GetIndex()];
+ G4double rval = G4UniformRand() * ((*dv)[nElements]);
+
+ elm = (*theElementVector)[nElements];
+ for (G4int i = 0; i < nElements; ++i) {
+ if (rval <= (*dv)[i]) {
+ elm = (*theElementVector)[i];
+ break;
+ }
+ }
+ } else {
+ elm = (*theElementVector)[0];
+ }
+
+ SetCurrentElement(elm);
+ return elm;
+}
+
+void G4muDarkBremsstrahlungModel::SetMethod(std::string method_in) {
+ method = method_in;
+ return;
+}