* add support in GlueXPhysicsList physics list initialization of

  the GlueXBernardConversionProcess to switch from e+e- to mu+mu-
  pair conversion, based on the first GENBEAM argument being
  'BHgen' for e+e- or 'BHmuons' for mu+mu-. [rtj]
* add a support to GlueXBeamConversionProcess to control the
  switching between e+e- and mu+mu- pair conversion at startup. [rtj]

Author: rjones30

src/GlueXBernardConversionProcess.cc

@@ -13,6 +13,11 @@
 #include "GlueXUserOptions.hh"
 #include "G4ParallelWorldProcess.hh"
 
+#include "G4Positron.hh"
+#include "G4Electron.hh"
+#include "G4MuonPlus.hh"
+#include "G4MuonMinus.hh"
+
 // If you set this flag to 1 then all beam photons that reach
 // the TPOL converter target will convert to e+e- pairs inside,
 // otherwise the standard pair conversion probabilities apply.
@@ -43,6 +48,7 @@ G4Mutex GlueXBernardConversionProcess::fMutex = G4MUTEX_INITIALIZER;
 int GlueXBernardConversionProcess::fStopBeamBeforeConverter = 0;
 int GlueXBernardConversionProcess::fStopBeamAfterConverter = 0;
 int GlueXBernardConversionProcess::fStopBeamAfterTarget = 0;
+int GlueXBernardConversionProcess::fLeptonPairFamily = 0;
 int GlueXBernardConversionProcess::fConfigured = 0;
 
 GlueXBernardConversionProcess::GlueXBernardConversionProcess(
@@ -98,6 +104,14 @@ GlueXBernardConversionProcess::GlueXBernardConversionProcess(
       {
          fStopBeamAfterTarget = 1;
       }
+      else if (genbeampars.find(1) != genbeampars.end() &&
+              (genbeampars[1] == "BHmuons_Bernard" ||
+               genbeampars[1] == "BHmuons_BERNARD" ||
+               genbeampars[1] == "bhmuons_bernard" ))
+      {
+         fStopBeamAfterTarget = 1;
+         fLeptonPairFamily = 1;
+      }
    }
 
    fConfigured = 1;
@@ -312,3 +326,19 @@ void GlueXBernardConversionProcess::GenerateConversionVertex(const G4Track &trac
       }
    }
 }
+
+G4ParticleDefinition *GlueXBernardConversionProcess::GetLepton(int charge) {
+   if (fLeptonPairFamily == 0) {
+      if (charge > 0)
+         return G4Positron::Definition();
+      else
+         return G4Electron::Definition();
+   }
+   else if (fLeptonPairFamily == 1) {
+      if (charge > 0)
+         return G4MuonPlus::Definition();
+      else
+         return G4MuonMinus::Definition();
+   }
+   return 0;
+}

src/GlueXBernardConversionProcess.hh

@@ -37,6 +37,8 @@ class GlueXBernardConversionProcess: public G4VEmProcess
    virtual G4double MinPrimaryEnergy(const G4ParticleDefinition*,
 				    const G4Material*) override;
 
+   virtual G4ParticleDefinition *GetLepton(int charge);
+
    // Print few lines of informations about the process: validity range,
    virtual void PrintInfo() override;
 
@@ -57,6 +59,7 @@ class GlueXBernardConversionProcess: public G4VEmProcess
    static int fStopBeamBeforeConverter;
    static int fStopBeamAfterConverter;
    static int fStopBeamAfterTarget;
+   static int fLeptonPairFamily;
 
  private:
    GlueXBernardConversionProcess operator=(GlueXBernardConversionProcess &src) = delete;

src/GlueXPhysicsList.cc

@@ -243,6 +243,8 @@ void GlueXPhysicsList::ConstructProcess()
             mgr->AddDiscreteProcess(fBernardConversion);
             G4VEmModel *model5D = new G4BetheHeitler5DModel();
             fBernardConversion->SetEmModel(model5D);
+            model5D->SetLeptonPair(fBernardConversion->GetLepton(+1),
+                                   fBernardConversion->GetLepton(-1));
 #else
             G4cerr << "Error in GlueXPhysicsList::ConstructProcess - "
                    << "BernardConversion process requested, but build "