Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Clean-up SIMP particles simulation #28814

Merged
merged 4 commits into from Jan 30, 2020
Merged

Clean-up SIMP particles simulation #28814

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
99d5d53
clean-up simp particles
civanch Jan 28, 2020
de44446
code format
civanch Jan 28, 2020
6e52c6b
code format
civanch Jan 28, 2020
3c87f42
avoid using methods changed between G4 releases
civanch Jan 29, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
35 changes: 9 additions & 26 deletions SimG4Core/CustomPhysics/interface/CMSQGSPSIMPBuilder.h
View file
Open in desktop
Expand Up @@ -3,44 +3,27 @@

#include "globals.hh"

#include "G4HadronElasticProcess.hh"
#include "G4HadronFissionProcess.hh"
#include "G4HadronCaptureProcess.hh"
#include "SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h"

#include "G4TheoFSGenerator.hh"
#include "G4PreCompoundModel.hh"
#include "G4GeneratorPrecompoundInterface.hh"
#include "G4QGSModel.hh"
#include "G4QGSParticipants.hh"
#include "G4QGSMFragmentation.hh"
#include "G4ExcitedStringDecay.hh"
#include "G4QuasiElasticChannel.hh"
#include "G4QGSModel.hh"

class CMSQGSPSIMPBuilder {
public:
CMSQGSPSIMPBuilder(G4bool quasiElastic = false);
virtual ~CMSQGSPSIMPBuilder();
class CMSSIMPInelasticProcess;
class G4QGSParticipants;
class G4QGSMFragmentation;
class G4ExcitedStringDecay;

class CMSQGSPSIMPBuilder {
public:
virtual void Build(G4HadronElasticProcess* aP);
virtual void Build(G4HadronFissionProcess* aP);
virtual void Build(G4HadronCaptureProcess* aP);
virtual void Build(CMSSIMPInelasticProcess* aP);
CMSQGSPSIMPBuilder();
~CMSQGSPSIMPBuilder();

void SetMinEnergy(G4double aM) { theMin = aM; }
void Build(CMSSIMPInelasticProcess* aP);

private:
G4TheoFSGenerator* theModel;
G4PreCompoundModel* thePreEquilib;
G4GeneratorPrecompoundInterface* theCascade;
G4QGSModel<G4QGSParticipants>* theStringModel;
G4ExcitedStringDecay* theStringDecay;
G4QuasiElasticChannel* theQuasiElastic;

G4QGSMFragmentation* theQGSM;

G4double theMin;
};

#endif
140 changes: 0 additions & 140 deletions SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h
View file
Open in desktop
Expand Up @@ -13,155 +13,15 @@ class CMSSIMPInelasticProcess : public G4HadronicProcess {

G4bool IsApplicable(const G4ParticleDefinition& aParticleType) override;

// register generator of secondaries
void RegisterMe(G4HadronicInteraction* a);

// get cross section per element
inline G4double GetElementCrossSection(const G4DynamicParticle* part,
const G4Element* elm,
const G4Material* mat = nullptr) {
G4double x = theCrossSectionDataStore->GetCrossSection(part, elm, mat);
if (x < 0.0) {
x = 0.0;
}
return x;
}

// obsolete method to get cross section per element
inline G4double GetMicroscopicCrossSection(const G4DynamicParticle* part,
const G4Element* elm,
const G4Material* mat = nullptr) {
return GetElementCrossSection(part, elm, mat);
}

// generic PostStepDoIt recommended for all derived classes
G4VParticleChange* PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) override;

// initialisation of physics tables and CMSSIMPInelasticProcessStore
void PreparePhysicsTable(const G4ParticleDefinition&) override;

// build physics tables and print out the configuration of the process
void BuildPhysicsTable(const G4ParticleDefinition&) override;

// dump physics tables
inline void DumpPhysicsTable(const G4ParticleDefinition& p) { theCrossSectionDataStore->DumpPhysicsTable(p); }

// add cross section data set
inline void AddDataSet(G4VCrossSectionDataSet* aDataSet) { theCrossSectionDataStore->AddDataSet(aDataSet); }

// access to the manager
inline G4EnergyRangeManager* GetManagerPointer() { return &theEnergyRangeManager; }

// get inverse cross section per volume
G4double GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*) override;

// access to the target nucleus
inline const G4Nucleus* GetTargetNucleus() const { return &targetNucleus; }

// G4ParticleDefinition* GetTargetDefinition();
inline const G4Isotope* GetTargetIsotope() { return targetNucleus.GetIsotope(); }

void ProcessDescription(std::ostream& outFile) const override;

protected:
// generic method to choose secondary generator
// recommended for all derived classes
inline G4HadronicInteraction* ChooseHadronicInteraction(G4double kineticEnergy,
G4Material* aMaterial,
G4Element* anElement) {
return theEnergyRangeManager.GetHadronicInteraction(kineticEnergy, aMaterial, anElement);
}

// access to the target nucleus
inline G4Nucleus* GetTargetNucleusPointer() { return &targetNucleus; }

public:
void BiasCrossSectionByFactor(G4double aScale);

// Energy-momentum non-conservation limits and reporting
inline void SetEpReportLevel(G4int level) { epReportLevel = level; }

inline void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double absoluteLevel) {
epCheckLevels.first = relativeLevel;
epCheckLevels.second = absoluteLevel;
levelsSetByProcess = true;
}

inline std::pair<G4double, G4double> GetEnergyMomentumCheckLevels() const { return epCheckLevels; }

// access to the cross section data store
inline G4CrossSectionDataStore* GetCrossSectionDataStore() { return theCrossSectionDataStore; }

inline void MultiplyCrossSectionBy(G4double factor) { aScaleFactor = factor; }

protected:
void DumpState(const G4Track&, const G4String&, G4ExceptionDescription&);

// obsolete method will be removed
inline const G4EnergyRangeManager& GetEnergyRangeManager() const { return theEnergyRangeManager; }

// obsolete method will be removed
inline void SetEnergyRangeManager(const G4EnergyRangeManager& value) { theEnergyRangeManager = value; }

// access to the chosen generator
inline G4HadronicInteraction* GetHadronicInteraction() const { return theInteraction; }

// access to the cross section data set
inline G4double GetLastCrossSection() { return theLastCrossSection; }

// fill result
void FillResult(G4HadFinalState* aR, const G4Track& aT);

// Check the result for catastrophic energy non-conservation
G4HadFinalState* CheckResult(const G4HadProjectile& thePro,
const G4Nucleus& targetNucleus,
G4HadFinalState* result) const;

// Check 4-momentum balance
void CheckEnergyMomentumConservation(const G4Track&, const G4Nucleus&);

private:
G4double XBiasSurvivalProbability();
G4double XBiasSecondaryWeight();

// hide assignment operator as private
CMSSIMPInelasticProcess& operator=(const CMSSIMPInelasticProcess& right) = delete;
CMSSIMPInelasticProcess(const CMSSIMPInelasticProcess&) = delete;

// Set E/p conservation check levels from environment variables
void GetEnergyMomentumCheckEnvvars();

protected:
G4HadProjectile thePro;

G4ParticleChange* theTotalResult;

G4int epReportLevel;

private:
G4ParticleDefinition* theParticle;

G4EnergyRangeManager theEnergyRangeManager;

G4HadronicInteraction* theInteraction;

G4CrossSectionDataStore* theCrossSectionDataStore;

G4Nucleus targetNucleus;

bool CMSSIMPInelasticProcess_debug_flag;

// Energy-momentum checking
std::pair<G4double, G4double> epCheckLevels;
G4bool levelsSetByProcess;

std::vector<G4VLeadingParticleBiasing*> theBias;

G4double theInitialNumberOfInteractionLength;

G4double aScaleFactor;
G4bool xBiasOn;
G4double theLastCrossSection;
};

#endif
34 changes: 3 additions & 31 deletions SimG4Core/CustomPhysics/interface/CMSSIMPInelasticXS.h
View file
Open in desktop
Expand Up @@ -3,17 +3,8 @@

#include "G4VCrossSectionDataSet.hh"
#include "globals.hh"
#include "G4ElementData.hh"
#include <vector>

const G4int MAXZINEL = 93;

class G4DynamicParticle;
class G4ParticleDefinition;
class G4Element;
class G4PhysicsVector;
class G4ComponentGGHadronNucleusXsc;
class G4HadronNucleonXsc;
class G4NeutronInelasticXS;

class CMSSIMPInelasticXS : public G4VCrossSectionDataSet {
public:
Expand All @@ -34,36 +25,17 @@ class CMSSIMPInelasticXS : public G4VCrossSectionDataSet {
const G4Element* elm,
const G4Material* mat) override;

G4Isotope* SelectIsotope(const G4Element*, G4double kinEnergy) override;

void BuildPhysicsTable(const G4ParticleDefinition&) override;

void CrossSectionDescription(std::ostream&) const override;

private:
void Initialise(G4int Z, G4DynamicParticle* dp = nullptr, const char* = nullptr);

G4PhysicsVector* RetrieveVector(std::ostringstream& in, G4bool warn);

G4double IsoCrossSection(G4double ekin, G4int Z, G4int A);

CMSSIMPInelasticXS& operator=(const CMSSIMPInelasticXS& right) = delete;
CMSSIMPInelasticXS(const CMSSIMPInelasticXS&) = delete;

G4ComponentGGHadronNucleusXsc* ggXsection;
G4HadronNucleonXsc* fNucleon;

const G4ParticleDefinition* proton;

G4ElementData data;
std::vector<G4PhysicsVector*> work;
std::vector<G4double> temp;
std::vector<G4double> coeff;

G4NeutronInelasticXS* nXsection;
const G4ParticleDefinition* neutron;
G4bool isInitialized;

static const G4int amin[MAXZINEL];
static const G4int amax[MAXZINEL];
};

#endif
58 changes: 19 additions & 39 deletions SimG4Core/CustomPhysics/src/CMSQGSPSIMPBuilder.cc
View file
Open in desktop
@@ -1,57 +1,37 @@

#include "SimG4Core/CustomPhysics/interface/CMSQGSPSIMPBuilder.h"
#include "SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h"

#include "G4SystemOfUnits.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"
#include "G4ProcessManager.hh"
#include "G4ExcitationHandler.hh"

#include "SimG4Core/CustomPhysics/interface/CMSSIMPInelasticXS.h"
#include "SimG4Core/CustomPhysics/interface/CMSSIMP.h"

CMSQGSPSIMPBuilder::CMSQGSPSIMPBuilder(G4bool quasiElastic) {
theMin = 12 * GeV;

theModel = new G4TheoFSGenerator("QGSP");

#include "G4TheoFSGenerator.hh"
#include "G4PreCompoundModel.hh"
#include "G4GeneratorPrecompoundInterface.hh"
#include "G4QGSParticipants.hh"
#include "G4QGSMFragmentation.hh"
#include "G4ExcitedStringDecay.hh"

CMSQGSPSIMPBuilder::CMSQGSPSIMPBuilder() {
theStringModel = new G4QGSModel<G4QGSParticipants>;
theStringDecay = new G4ExcitedStringDecay(theQGSM = new G4QGSMFragmentation);
theStringModel->SetFragmentationModel(theStringDecay);

theCascade = new G4GeneratorPrecompoundInterface;
thePreEquilib = new G4PreCompoundModel(new G4ExcitationHandler);
theCascade->SetDeExcitation(thePreEquilib);

theModel->SetTransport(theCascade);
theModel->SetHighEnergyGenerator(theStringModel);
if (quasiElastic) {
theQuasiElastic = new G4QuasiElasticChannel;
theModel->SetQuasiElasticChannel(theQuasiElastic);
} else {
theQuasiElastic = nullptr;
}
}

CMSQGSPSIMPBuilder::~CMSQGSPSIMPBuilder() {
delete theStringDecay;
delete theStringModel;
delete thePreEquilib;
delete theCascade;
if (theQuasiElastic)
delete theQuasiElastic;
delete theModel;
delete theQGSM;
}

void CMSQGSPSIMPBuilder::Build(G4HadronElasticProcess*) {}

void CMSQGSPSIMPBuilder::Build(G4HadronFissionProcess*) {}

void CMSQGSPSIMPBuilder::Build(G4HadronCaptureProcess*) {}

void CMSQGSPSIMPBuilder::Build(CMSSIMPInelasticProcess* aP) {
theModel->SetMinEnergy(theMin);
theModel->SetMaxEnergy(100 * TeV);
G4GeneratorPrecompoundInterface* theCascade = new G4GeneratorPrecompoundInterface;
G4PreCompoundModel* thePreEquilib = new G4PreCompoundModel();
theCascade->SetDeExcitation(thePreEquilib);

G4TheoFSGenerator* theModel = new G4TheoFSGenerator("QGSP");
theModel->SetTransport(theCascade);
theModel->SetHighEnergyGenerator(theStringModel);
theModel->SetMinEnergy(0.0);
theModel->SetMaxEnergy(100 * CLHEP::TeV);
aP->RegisterMe(theModel);
aP->AddDataSet(new CMSSIMPInelasticXS());
}