Merge pull request #38901 from revering/DBremSim_10_6_X

Dark Bremsstrahlung backport to 10_6_X

SimG4Core/CustomPhysics/BuildFile.xml

@@ -1,11 +1,16 @@
 <export>
   <lib   name="1"/>
 </export>
-<use   name="FWCore/Framework"/>
-<use   name="FWCore/PluginManager"/>
-<use   name="FWCore/MessageLogger"/>
-<use   name="SimG4Core/Physics"/>
-<use   name="SimG4Core/PhysicsLists"/>
-<use   name="geant4core"/>
-<use   name="clhep"/>
-<use   name="boost"/>
+<use name="DataFormats/GeometryVector" source_only="1"/>
+<use name="FWCore/Framework"/>
+<use name="FWCore/MessageLogger"/>
+<use name="FWCore/PluginManager"/>
+<use name="SimG4Core/Notification"/>
+<use name="SimG4Core/Physics"/>
+<use name="SimG4Core/PhysicsLists"/>
+<use name="SimG4Core/Watcher"/>
+<use name="geant4core"/>
+<use name="clhep"/>
+<use name="boost"/>
+<use name="rootmath"/>
+<use name="root"/>

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

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 <vector>
+#include <tuple>
+
+class DBremWatcher : public SimProducer,
+                     public Observer<const BeginOfTrack *>,
+                     public Observer<const BeginOfEvent *>,
+                     public Observer<const BeginOfRun *>,
+                     public Observer<const EndOfEvent *>,
+                     public Observer<const EndOfTrack *> {
+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<int> pdgs_;
+  int MotherId;
+  float m_weight;
+  double biasFactor;
+  bool foundbrem;
+  G4ThreeVector aPrimeTraj;
+  G4ThreeVector finaltraj;
+  G4ThreeVector VertexPos;
+  float f_energy;
+};
+
+#endif

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

SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlung.h

@@ -0,0 +1,38 @@
+/**
+ * @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 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

SimG4Core/CustomPhysics/interface/G4muDarkBremsstrahlungModel.h

@@ -0,0 +1,100 @@
+/**
+ * @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 MMu;
+  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<G4DynamicParticle*>*,
+                         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;
+  G4double muonMass;
+  G4bool isMuon;
+
+private:
+  G4double highKinEnergy;
+  G4double lowKinEnergy;
+  G4double probsup;
+  G4bool isInitialised;
+  std::string method;
+  G4bool mg_loaded;
+  std::map<double, std::vector<frame> > mgdata;
+  std::vector<std::pair<double, int> > energies;
+  std::vector<G4DataVector*> partialSumSigma;
+};
+
+#endif

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) 
+)

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);
+}

SimG4Core/CustomPhysics/src/CustomPhysics.cc

@@ -3,6 +3,7 @@
 #include "SimG4Core/CustomPhysics/interface/CustomPhysicsListSS.h"
 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsLPM.h"
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "SimG4Core/CustomPhysics/interface/APrimePhysics.h"
 
 #include "G4DecayPhysics.hh"
 #include "G4EmExtraPhysics.hh"
@@ -21,6 +22,7 @@ CustomPhysics::CustomPhysics(const edm::ParameterSet& p) : PhysicsList(p) {
   int ver = p.getUntrackedParameter<int>("Verbosity", 0);
   bool tracking = p.getParameter<bool>("TrackingCut");
   bool ssPhys = p.getUntrackedParameter<bool>("ExoticaPhysicsSS", false);
+  bool dbrem = p.getUntrackedParameter<bool>("DBrem", false);
   double timeLimit = p.getParameter<double>("MaxTrackTime") * ns;
   edm::LogInfo("PhysicsList") << "You are using the simulation engine: "
                               << "FTFP_BERT_EMM for regular particles \n"
@@ -55,7 +57,11 @@ CustomPhysics::CustomPhysics(const edm::ParameterSet& p) : PhysicsList(p) {
   }
 
   // Custom Physics
-  if (ssPhys) {
+  if (dbrem) {
+    RegisterPhysics(new APrimePhysics(p.getUntrackedParameter<double>("DBremMass"),
+                                      p.getUntrackedParameter<std::string>("DBremScaleFile"),
+                                      p.getUntrackedParameter<double>("DBremBiasFactor")));
+  } else if (ssPhys) {
     RegisterPhysics(new CustomPhysicsListSS("custom", p));
   } else {
     RegisterPhysics(new CustomPhysicsList("custom", p));