From 8949be389542dd493c2b269b359d87a22be7d4ce Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Fri, 16 Oct 2020 18:26:18 +0200
Subject: [PATCH] * Preparing a simple fusion Process in npsimulation -
all in place - need to do the proper kinematic calculation
---
NPLib/Core/NPInputParser.cxx | 2 +-
NPSimulation/Process/BeamReaction.cc | 78 ++++++++++++++++++++++------
NPSimulation/Process/BeamReaction.hh | 15 +++++-
3 files changed, 75 insertions(+), 20 deletions(-)
diff --git a/NPLib/Core/NPInputParser.cxx b/NPLib/Core/NPInputParser.cxx
index 13196d004..7962e3e80 100644
--- a/NPLib/Core/NPInputParser.cxx
+++ b/NPLib/Core/NPInputParser.cxx
@@ -384,7 +384,7 @@ void NPL::InputParser::TreatAliases(){
name += alias[i]->GetMainValue();
std::string action = alias[i]->GetString("Action");
std::vector<std::string> value = alias[i]->GetVectorString("Value");
- if(action=="Inplace" && value.size()!=1)
+ if(action=="Replace" && value.size()!=1)
NPL::SendErrorAndExit("NPL::InputParser", "Inplace alias can only take one value");
// Scan all blocks for aliases
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 1d73194f8..d3003cf3c 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -61,15 +61,16 @@ NPS::BeamReaction::BeamReaction(G4String modelName)
void NPS::BeamReaction::ReadConfiguration() {
NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
- if(input.GetAllBlocksWithToken("TwoBodyReaction").size()>0) m_ReactionType ="TwoBodyReaction";
- else if(input.GetAllBlocksWithToken("QFSReaction").size()>0) m_ReactionType ="QFSReaction";
+ if(input.GetAllBlocksWithToken("TwoBodyReaction").size()>0) m_ReactionType =TwoBody;
+ else if(input.GetAllBlocksWithToken("QFSReaction").size()>0) m_ReactionType =QFS;
+ else if(input.GetAllBlocksWithToken("FusionReaction").size()>0) m_ReactionType =Fusion;
-
- if (m_ReactionType=="TwoBodyReaction" ) {
+ // Two body
+ if (m_ReactionType==TwoBody ) {
m_Reaction.ReadConfigurationFile(input);
m_BeamName = NPL::ChangeNameToG4Standard(m_Reaction.GetParticle1()->GetName());
if(m_Reaction.GetParticle3()->GetName() != ""){
- m_active = true;
+ m_active = true;
m_ReactionConditions = new TReactionConditions();
AttachReactionConditions();
if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
@@ -78,16 +79,29 @@ void NPS::BeamReaction::ReadConfiguration() {
}
}
- else if (m_ReactionType=="QFSReaction") {
+ // QFS
+ else if (m_ReactionType==QFS) {
m_QFS.ReadConfigurationFile(input);
m_BeamName = NPL::ChangeNameToG4Standard(m_QFS.GetParticleA()->GetName());
- m_active = true;
+ m_active = true;
m_ReactionConditions = new TReactionConditions();
AttachReactionConditions();
if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
RootOutput::getInstance()->GetTree()->Branch(
"ReactionConditions", "TReactionConditions", &m_ReactionConditions);
}
+
+ // Fusion
+ else if (m_ReactionType==Fusion) {
+ vector<InputBlock*> blocks=input.GetAllBlocksWithToken("FusionReaction");
+ m_BeamName = NPL::ChangeNameToG4Standard(blocks[0]->GetString("Beam"));
+ m_TargetNuclei = blocks[0]->GetString("Target");
+ m_FusionProduct = blocks[0]->GetString("Product");
+ m_FusionExcitation = blocks[0]->GetDouble("ExcitationEnergy","MeV");
+ m_active = true;
+ // not used
+ m_ReactionConditions = new TReactionConditions();
+ }
else {
m_active = false;
}
@@ -153,7 +167,7 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
m_StepSize = PrimaryTrack->GetStepLength();
// If the condition is met, the event is generated
if (m_shoot && m_length < m_StepSize) {
- if(m_ReactionType=="QFSReaction"){
+ if(m_ReactionType==QFS){
if ( m_QFS.IsAllowed() ) {
return true;
}
@@ -163,7 +177,7 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
}
}
- else if(m_ReactionType=="TwoBodyReaction"){
+ else if(m_ReactionType==TwoBody){
if ( m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy()) ) {
return true;
}
@@ -172,6 +186,9 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
std::cout << "Two body reaction not allowed" << std::endl;
}
}
+ else if(m_ReactionType==Fusion){
+ return true;
+ }
}
return false;
@@ -220,12 +237,11 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
fastStep.KillPrimaryTrack();
fastStep.SetPrimaryTrackPathLength(0.0);
+ ///////////////////////////////
+ // Two-Body Reaction Case /////
+ ///////////////////////////////
+ if (m_ReactionType==TwoBody) {
- if (m_ReactionType=="TwoBodyReaction" ) {
-
- ///////////////////////////////
- // Two-Body Reaction Case /////
- ///////////////////////////////
static G4IonTable* IonTable
= G4ParticleTable::GetParticleTable()->GetIonTable();
@@ -398,7 +414,9 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
}// end if TwoBodyReaction
- else if (m_ReactionType=="QFSReaction" ) {
+
+ // QFS
+ else if (m_ReactionType==QFS ) {
//////Define the kind of particle to shoot////////
// A --> T ==> B + (c -> T) => B + 1 + 2
@@ -576,9 +594,35 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
m_ReactionConditions->SetMomentumDirectionZ(momentum_kine2_world.z());
m_ReactionConditions->SetMomentumDirectionZ(momentum_kineB_world.z());
+ }// end QFS
+
+ ////////////////////
+ // Fusion Case //
+ ////////////////////
+ else if(m_ReactionType==Fusion){
+ // Set the end of the step conditions
+ fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0, pdirection);
+ fastStep.SetPrimaryTrackFinalPosition(worldPosition);
+ fastStep.SetTotalEnergyDeposited(0);
+ fastStep.SetPrimaryTrackFinalTime(time);// FIXME
+ fastStep.KillPrimaryTrack();
+ fastStep.SetPrimaryTrackPathLength(0.0);
+
+ static G4IonTable* IonTable
+ = G4ParticleTable::GetParticleTable()->GetIonTable();
+ //////Define the kind of particle to shoot////////
+ G4ParticleDefinition* Product;
+ NPL::Particle N(m_FusionProduct);
+ int PZ = N.GetZ();
+ int PA = N.GetA();
+ Product = IonTable->GetIon(PZ, PA, m_FusionExcitation*MeV);
+ // setup the daugter
+ G4ThreeVector momentum_dir(0,0,1);
+ double Energy=energy;
+ G4DynamicParticle particle(Product, momentum_dir, Energy);
+ fastStep.CreateSecondaryTrack(particle, localPosition, time);
+ }// end fusion
-
- }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Return the slow down beam in given thickness of material
diff --git a/NPSimulation/Process/BeamReaction.hh b/NPSimulation/Process/BeamReaction.hh
index 212d64f1c..e127c74c2 100644
--- a/NPSimulation/Process/BeamReaction.hh
+++ b/NPSimulation/Process/BeamReaction.hh
@@ -31,6 +31,12 @@
#include "TReactionConditions.h"
class G4VPhysicalVolume;
namespace NPS{
+ enum ReactionType{
+ TwoBody,
+ QFS,
+ Fusion
+ };
+
class BeamReaction : public G4VFastSimulationModel{
public:
BeamReaction (G4String, G4Region*);
@@ -47,7 +53,8 @@ namespace NPS{
NPL::Reaction m_Reaction;
NPL::QFS m_QFS;
string m_BeamName;
- string m_ReactionType;
+ int m_ReactionType;
+
bool m_active;// is the process active
bool m_shoot;
double m_StepSize;
@@ -57,7 +64,11 @@ namespace NPS{
double m_length;
int m_Parent_ID;
double SlowDownBeam(const G4ParticleDefinition* Beam, double IncidentEnergy, double Thickness,G4Material* Material);
-
+
+ private:// specific for the simple case of fusion
+ string m_TargetNuclei;
+ string m_FusionProduct;
+ double m_FusionExcitation;
private:
TReactionConditions* m_ReactionConditions;
--
GitLab