Skip to content

Commit

Permalink
backport
Browse files Browse the repository at this point in the history
  • Loading branch information
@civanch
civanch committed Mar 4, 2019
1 parent c4c9bd4 commit 770e744
Show file tree
Hide file tree
Showing 5 changed files with 658 additions and 12 deletions.
69 changes: 69 additions & 0 deletions SimG4Core/PhysicsLists/interface/CMSmplIonisation.h
View file
@@ -0,0 +1,69 @@
//
// -------------------------------------------------------------------
//
// 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;

// print description in html
void ProcessDescription(std::ostream&) const override;

protected:

void InitialiseEnergyLossProcess(const G4ParticleDefinition*,
const G4ParticleDefinition*) override;

private:

G4double magneticCharge;
G4bool isInitialised;
};

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

#endif
114 changes: 114 additions & 0 deletions SimG4Core/PhysicsLists/interface/CMSmplIonisationWithDeltaModel.h
View file
@@ -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....
21 changes: 9 additions & 12 deletions SimG4Core/PhysicsLists/src/CMSMonopolePhysics.cc
View file
@@ -1,25 +1,22 @@
#include "SimG4Core/PhysicsLists/interface/CMSMonopolePhysics.h"
#include "SimG4Core/PhysicsLists/interface/MonopoleTransportation.h"
#include "SimG4Core/MagneticField/interface/ChordFinderSetter.h"
#include "SimG4Core/PhysicsLists/interface/CMSmplIonisation.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "G4ParticleDefinition.hh"
#include "G4ProcessManager.hh"

#include "G4StepLimiter.hh"
#include "G4mplIonisation.hh"
#include "G4mplIonisationWithDeltaModel.hh"
#include "G4hMultipleScattering.hh"
#include "G4hIonisation.hh"
#include "G4hhIonisation.hh"

#include "CLHEP/Units/GlobalSystemOfUnits.h"

CMSMonopolePhysics::CMSMonopolePhysics(const HepPDT::ParticleDataTable * pdt,
sim::ChordFinderSetter * cfs,
const edm::ParameterSet & p) :
G4VPhysicsConstructor("Monopole Physics"), chordFinderSetter(cfs) {

G4VPhysicsConstructor("Monopole Physics")
{
verbose = p.getUntrackedParameter<int>("Verbosity",0);
magCharge = p.getUntrackedParameter<int>("MonopoleCharge",1);
deltaRay = p.getUntrackedParameter<bool>("MonopoleDeltaRay",true);
Expand Down Expand Up @@ -97,8 +94,8 @@ void CMSMonopolePhysics::ConstructProcess() {
if(!pmanager) {
std::ostringstream o;
o << "Monopole without a Process Manager";
G4Exception("CMSMonopolePhysics::ConstructProcess()","",
FatalException,o.str().c_str());
throw edm::Exception( edm::errors::Configuration, o.str().c_str());
return;
}

G4double magn = mpl->MagneticCharge();
Expand All @@ -112,9 +109,9 @@ void CMSMonopolePhysics::ConstructProcess() {
}

if (magn != 0.0) {
pmanager->RemoveProcess(0);
pmanager->AddProcess(new MonopoleTransportation(mpl,chordFinderSetter,verbose),
-1, 0, 0);
G4int idxt(0);
pmanager->RemoveProcess(idxt);
pmanager->AddProcess(new MonopoleTransportation(mpl,verbose),-1,0,0);
}

if (mpl->GetPDGCharge() != 0.0) {
Expand All @@ -126,7 +123,7 @@ void CMSMonopolePhysics::ConstructProcess() {
ph->RegisterProcess(hioni, mpl);
}
if(magn != 0.0) {
G4mplIonisation* mplioni = new G4mplIonisation(magn);
CMSmplIonisation* mplioni = new CMSmplIonisation(magn);
ph->RegisterProcess(mplioni, mpl);
}
pmanager->AddDiscreteProcess(new G4StepLimiter());
Expand Down
114 changes: 114 additions & 0 deletions SimG4Core/PhysicsLists/src/CMSmplIonisation.cc
View file
@@ -0,0 +1,114 @@
//
// -------------------------------------------------------------------
//
// 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"
#include "G4EmParameters.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
G4EmParameters* param = G4EmParameters::Instance();
G4double emin = std::min(param->MinKinEnergy(),ion->LowEnergyLimit());
G4double emax = std::max(param->MaxKinEnergy(),ion->HighEnergyLimit());
G4int bin = G4lrint(param->NumberOfBinsPerDecade()*std::log10(emax/emin));
ion->SetLowEnergyLimit(emin);
ion->SetHighEnergyLimit(emax);
SetMinKinEnergy(emin);
SetMaxKinEnergy(emax);
SetDEDXBinning(bin);

SetEmModel(ion);
AddEmModel(1,ion,ion);

isInitialised = true;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void CMSmplIonisation::PrintInfo()
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void CMSmplIonisation::ProcessDescription(std::ostream& out) const
{
out << "No description available." << G4endl;
G4VEnergyLossProcess::ProcessDescription(out);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

0 comments on commit 770e744

Please sign in to comment.