Merge pull request #26092 from civanch/monopole_ionisation_fixed

Backport fix of monopole ionisation

SimG4Core/PhysicsLists/interface/CMSmplIonisation.h

@@ -0,0 +1,66 @@
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class header file
+//
+//
+// File name:     CMSmplIonisation
+//
+// Author:        Vladimir Ivanchenko copied from Geant4 10.5p01
+//
+// Creation date: 02.03.2019
+//
+// Class Description:
+//
+// This class manages the ionisation process for a magnetic monopole
+// it inherites from G4VContinuousDiscreteProcess via G4VEnergyLossProcess.
+// Magnetic charge of the monopole should be defined in the constructor of 
+// the process, unless it is assumed that it is classic Dirac monopole with 
+// the charge 67.5*eplus. The name of the particle should be "monopole".
+//
+
+// -------------------------------------------------------------------
+//
+
+#ifndef CMSmplIonisation_h
+#define CMSmplIonisation_h 1
+
+#include "G4VEnergyLossProcess.hh"
+#include "globals.hh"
+#include "G4VEmModel.hh"
+
+class G4Material;
+class G4VEmFluctuationModel;
+
+class CMSmplIonisation : public G4VEnergyLossProcess
+{
+
+public:
+
+  explicit CMSmplIonisation(G4double mCharge = 0.0, 
+                           const G4String& name = "mplIoni");
+
+  ~CMSmplIonisation() override;
+
+  G4bool IsApplicable(const G4ParticleDefinition& p) override;
+
+  G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
+                                    const G4Material*, G4double cut) final;
+
+  // Print out of the class parameters
+  void PrintInfo() override;
+
+protected:
+
+  void InitialiseEnergyLossProcess(const G4ParticleDefinition*,
+                                   const G4ParticleDefinition*) override;
+
+private:
+
+  G4double    magneticCharge;
+  G4bool      isInitialised;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+#endif

SimG4Core/PhysicsLists/interface/CMSmplIonisationWithDeltaModel.h

@@ -0,0 +1,114 @@
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class header file
+//
+//
+// File name:     CMSmplIonisationWithDeltaModel
+//
+// Author:        Vladimir Ivanchenko copied from Geant4 10.5p01
+//
+// Creation date: 02.03.2019
+//
+// Class Description:
+//
+// Implementation of model of energy loss of the magnetic monopole
+
+// -------------------------------------------------------------------
+//
+
+#ifndef CMSmplIonisationWithDeltaModel_h
+#define CMSmplIonisationWithDeltaModel_h 1
+
+#include "G4VEmModel.hh"
+#include "G4VEmFluctuationModel.hh"
+#include <vector>
+
+class G4ParticleChangeForLoss;
+
+class CMSmplIonisationWithDeltaModel : public G4VEmModel, public G4VEmFluctuationModel
+{
+
+public:
+
+  explicit CMSmplIonisationWithDeltaModel(G4double mCharge, 
+                                const G4String& nam = "mplIonisationWithDelta");
+
+  ~CMSmplIonisationWithDeltaModel() override;
+
+  void Initialise(const G4ParticleDefinition*, 
+                          const G4DataVector&) override;
+
+  G4double ComputeDEDXPerVolume(const G4Material*,
+                                        const G4ParticleDefinition*,
+                                        G4double kineticEnergy,
+                                        G4double cutEnergy) override;
+
+  virtual G4double ComputeCrossSectionPerElectron(
+                                 const G4ParticleDefinition*,
+                                 G4double kineticEnergy,
+                                 G4double cutEnergy,
+                                 G4double maxEnergy);
+
+  G4double ComputeCrossSectionPerAtom(
+                                 const G4ParticleDefinition*,
+                                 G4double kineticEnergy,
+                                 G4double Z, G4double A,
+                                 G4double cutEnergy,
+                                 G4double maxEnergy) override;
+
+  void SampleSecondaries(std::vector<G4DynamicParticle*>*,
+                                 const G4MaterialCutsCouple*,
+                                 const G4DynamicParticle*,
+                                 G4double tmin,
+                                 G4double maxEnergy) override;
+
+
+  G4double SampleFluctuations(const G4MaterialCutsCouple*,
+                                      const G4DynamicParticle*,
+                                      G4double tmax,
+                                      G4double length,
+                                      G4double meanLoss) override;
+
+  G4double Dispersion(const G4Material*,
+                              const G4DynamicParticle*,
+                              G4double tmax,
+                              G4double length) override;
+
+  G4double MinEnergyCut(const G4ParticleDefinition*,
+                                const G4MaterialCutsCouple* couple) override;
+
+  void SetParticle(const G4ParticleDefinition* p);
+
+protected:
+
+  G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
+                                      G4double kinEnergy) override;
+
+private:
+
+  G4double ComputeDEDXAhlen(const G4Material* material, G4double bg2, 
+                            G4double cut);
+
+  const G4ParticleDefinition* monopole;
+  G4ParticleDefinition*       theElectron;
+  G4ParticleChangeForLoss*    fParticleChange;
+
+  G4double mass;
+  G4double magCharge;
+  G4double twoln10;
+  G4double betalow;
+  G4double betalim;
+  G4double beta2lim;
+  G4double bg2lim;
+  G4double chargeSquare;
+  G4double dedxlim;
+  G4int    nmpl;
+  G4double pi_hbarc2_over_mc2;
+
+  static std::vector<G4double>* dedx0;
+};
+
+#endif
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

SimG4Core/PhysicsLists/src/CMSMonopolePhysics.cc

@@ -1,14 +1,14 @@
 #include "SimG4Core/PhysicsLists/interface/CMSMonopolePhysics.h"
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "SimG4Core/PhysicsLists/interface/CMSmplIonisation.h"
 
 #include "G4ParticleDefinition.hh"
 #include "G4ProcessManager.hh"
 
 #include "G4StepLimiter.hh"
 #include "G4Transportation.hh"
 #include "G4MonopoleTransportation.hh"
-#include "G4mplIonisation.hh"
-#include "G4mplIonisationWithDeltaModel.hh"
+#include "G4mplIonisationModel.hh"
 #include "G4hMultipleScattering.hh"
 #include "G4hIonisation.hh"
 #include "G4hhIonisation.hh"
@@ -148,13 +148,12 @@ void CMSMonopolePhysics::ConstructProcess() {
 	++idx;
       }
       if(magn != 0.0) {
-	G4mplIonisation* mplioni = new G4mplIonisation(magn);
+	CMSmplIonisation* mplioni = new CMSmplIonisation(magn);
 	mplioni->SetDEDXBinning(nbin);
 	mplioni->SetMinKinEnergy(emin);
 	mplioni->SetMaxKinEnergy(emax);
 	if (!deltaRay) {
-	  G4mplIonisationWithDeltaModel* mod = 
-	    new G4mplIonisationWithDeltaModel(magn,"PAI");
+	  G4mplIonisationModel* mod = new G4mplIonisationModel(magn,"PAI");
 	  mplioni->AddEmModel(0,mod,mod);
 	}
 	pmanager->AddProcess(mplioni, -1, idx, idx);

SimG4Core/PhysicsLists/src/CMSmplIonisation.cc

@@ -0,0 +1,100 @@
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:     CMSmplIonisation
+//
+// Author:        Vladimir Ivanchenko copied from Geant4 10.5p01
+//
+// Creation date: 02.03.2019
+//
+//
+// -------------------------------------------------------------------
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+#include "SimG4Core/PhysicsLists/interface/CMSmplIonisation.h"
+#include "SimG4Core/PhysicsLists/interface/CMSmplIonisationWithDeltaModel.h"
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4Electron.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+using namespace std;
+
+CMSmplIonisation::CMSmplIonisation(G4double mCharge, const G4String& name)
+  : G4VEnergyLossProcess(name),
+    magneticCharge(mCharge),
+    isInitialised(false)
+{
+  // By default classical magnetic charge is used
+  if(magneticCharge == 0.0) { magneticCharge = eplus*0.5/fine_structure_const; }
+
+  SetVerboseLevel(0);
+  SetProcessSubType(fIonisation);
+  SetStepFunction(0.2, 1*mm);
+  SetSecondaryParticle(G4Electron::Electron());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+CMSmplIonisation::~CMSmplIonisation()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool CMSmplIonisation::IsApplicable(const G4ParticleDefinition&)
+{
+  return true;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double CMSmplIonisation::MinPrimaryEnergy(const G4ParticleDefinition* mpl,
+                                           const G4Material*,
+                                           G4double cut)
+{
+  G4double x = 0.5*cut/electron_mass_c2;
+  G4double mass  = mpl->GetPDGMass();
+  G4double ratio = electron_mass_c2/mass;
+  G4double gam   = x*ratio + std::sqrt((1. + x)*(1. + x*ratio*ratio));
+  return mass*(gam - 1.0);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void CMSmplIonisation::InitialiseEnergyLossProcess(const G4ParticleDefinition* p,
+                                                   const G4ParticleDefinition*)
+{
+  if(isInitialised) { return; }
+
+  SetBaseParticle(nullptr);
+
+  // monopole model is responsible both for energy loss and fluctuations
+  CMSmplIonisationWithDeltaModel* ion =
+    new CMSmplIonisationWithDeltaModel(magneticCharge,"PAI");
+  ion->SetParticle(p);
+
+  // define size of dedx and range tables
+  G4double emin  = std::min(MinKinEnergy(),ion->LowEnergyLimit());
+  G4double emax  = std::max(MaxKinEnergy(),ion->HighEnergyLimit());
+  ion->SetLowEnergyLimit(emin);
+  ion->SetHighEnergyLimit(emax);
+  SetMinKinEnergy(emin);
+  SetMaxKinEnergy(emax);
+  AddEmModel(1,ion,ion);
+
+  isInitialised = true;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void CMSmplIonisation::PrintInfo()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+