diff --git a/src/GlueXBernardConversionProcess.cc b/src/GlueXBernardConversionProcess.cc index 4c6d624..e93e521 100644 --- a/src/GlueXBernardConversionProcess.cc +++ b/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; +} diff --git a/src/GlueXBernardConversionProcess.hh b/src/GlueXBernardConversionProcess.hh index aff5e5c..f868572 100644 --- a/src/GlueXBernardConversionProcess.hh +++ b/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; diff --git a/src/GlueXPhysicsList.cc b/src/GlueXPhysicsList.cc index 9813f20..c59b49e 100644 --- a/src/GlueXPhysicsList.cc +++ b/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 "