From fa8b4123a6b243327cc99a3b60a909a6849dd5d0 Mon Sep 17 00:00:00 2001
From: Adrien Matta <matta@lpccaen.in2p3.fr>
Date: Fri, 6 Sep 2024 16:50:55 +0200
Subject: [PATCH] * Inelastic break up now working
---
NPSimulation/Process/BeamReaction.cc | 182 +++++++++++++++++++++++++++
NPSimulation/Process/BeamReaction.hh | 4 +-
2 files changed, 185 insertions(+), 1 deletion(-)
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index b4bd84c6c..3f9bb927c 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -68,6 +68,8 @@ void NPS::BeamReaction::ReadConfiguration() {
if (input.GetAllBlocksWithToken("TwoBodyReaction").size() > 0)
m_ReactionType = TwoBody;
+ if (input.GetAllBlocksWithToken("InelasticBreakup").size() > 0)
+ m_ReactionType = InelasticBreakup;
else if (input.GetAllBlocksWithToken("QFSReaction").size() > 0)
m_ReactionType = QFS;
else if (input.GetAllBlocksWithToken("PhaseSpace").size() > 0)
@@ -89,6 +91,20 @@ void NPS::BeamReaction::ReadConfiguration() {
}
}
+ // Inelastic Breakup
+ if (m_ReactionType == InelasticBreakup) {
+ m_InelasticBreakup.ReadConfigurationFile(input);
+ m_BeamName = NPL::ChangeNameToG4Standard(m_InelasticBreakup.GetParticleBeam()->GetName());
+ if (m_InelasticBreakup.GetParticleLight()->GetName() != "") {
+ m_active = true;
+ m_ReactionConditions = new TReactionConditions();
+ AttachReactionConditions();
+ if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
+ RootOutput::getInstance()->GetTree()->Branch("ReactionConditions", "TReactionConditions",
+ &m_ReactionConditions);
+ }
+ }
+
// QFS
else if (m_ReactionType == QFS) {
m_QFS.ReadConfigurationFile(input);
@@ -209,6 +225,9 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
std::cout << "Two body reaction not allowed" << std::endl;
}
}
+ else if (m_ReactionType == InelasticBreakup) {
+ return true;
+ }
else if (m_ReactionType == Fusion) {
return true;
}
@@ -415,6 +434,169 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
m_ReactionConditions->SetMomentumDirectionZ(momentum_kine3_world.z());
m_ReactionConditions->SetMomentumDirectionZ(momentum_kine4_world.z());
+ } // end if TwoBodyReaction
+ ////////////////////////////////////////
+ // Inelastic Breakup Reaction Case /////
+ ////////////////////////////////////////
+ if (m_ReactionType == InelasticBreakup) {
+
+ static G4IonTable* IonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
+ //////Define the kind of particle to shoot////////
+ // Particle Light
+ G4ParticleDefinition* LightName;
+ int LightZ = m_InelasticBreakup.GetParticleLight()->GetZ();
+ int LightA = m_InelasticBreakup.GetParticleLight()->GetA();
+
+ if (LightZ == 0 && LightA == 1) {
+ LightName = G4Neutron::Definition();
+ }
+ else {
+ LightName = IonTable->GetIon(LightZ, LightA);
+ }
+
+ // Particle Heavy
+ G4int HeavyZ = m_InelasticBreakup.GetParticleHeavy()->GetZ();
+ G4int HeavyA = m_InelasticBreakup.GetParticleHeavy()->GetA();
+ G4ParticleDefinition* HeavyName = IonTable->GetIon(HeavyZ, HeavyA);
+
+ // Particle TargetLike
+ G4int TargetZ = m_InelasticBreakup.GetParticleTarget()->GetZ();
+ G4int TargetA = m_InelasticBreakup.GetParticleTarget()->GetA();
+ G4ParticleDefinition* TargetName = IonTable->GetIon(TargetZ, TargetA);
+
+ // Set the Energy of the reaction
+ m_InelasticBreakup.SetBeamEnergy(reac_energy);
+
+ double Beam_theta = pdirection.theta();
+ double Beam_phi = pdirection.phi();
+
+ ///////////////////////////
+ ///// Beam Parameters /////
+ ///////////////////////////
+ m_ReactionConditions->SetBeamParticleName(PrimaryTrack->GetParticleDefinition()->GetParticleName());
+
+ m_ReactionConditions->SetBeamReactionEnergy(reac_energy);
+ m_ReactionConditions->SetVertexPositionX(worldPosition.x());
+ m_ReactionConditions->SetVertexPositionY(worldPosition.y());
+ m_ReactionConditions->SetVertexPositionZ(worldPosition.z());
+
+ G4ThreeVector U(1, 0, 0);
+ G4ThreeVector V(0, 1, 0);
+ G4ThreeVector ZZ(0, 0, 1);
+ m_ReactionConditions->SetBeamEmittanceTheta(PrimaryTrack->GetMomentumDirection().theta() / deg);
+ m_ReactionConditions->SetBeamEmittancePhi(PrimaryTrack->GetMomentumDirection().phi() / deg);
+ m_ReactionConditions->SetBeamEmittanceThetaX(PrimaryTrack->GetMomentumDirection().angle(U) / deg);
+ m_ReactionConditions->SetBeamEmittancePhiY(PrimaryTrack->GetMomentumDirection().angle(V) / deg);
+
+ //////////////////////////////////////////////////////////
+ ///// Build rotation matrix to go from the incident //////
+ ///// beam frame to the "world" frame //////
+ //////////////////////////////////////////////////////////
+
+ // G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
+ // cos(Beam_theta) * sin(Beam_phi),
+ // -sin(Beam_theta));
+ // G4ThreeVector col2(-sin(Beam_phi),
+ // cos(Beam_phi),
+ // 0);
+ // G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
+ // sin(Beam_theta) * sin(Beam_phi),
+ // cos(Beam_theta));
+ // G4RotationMatrix BeamToWorld(col1, col2, col3);
+
+ /////////////////////////////////////////////////////////////////
+ ///// Angles for emitted particles following Cross Section //////
+ ///// Angles are in the beam frame //////
+ /////////////////////////////////////////////////////////////////
+ double phi = G4RandFlat::shoot() * 2. * pi;
+
+ //////////////////////////////////////////////////
+ ///// Momentum and angles from kinematics /////
+ ///// Angles are in the beam frame /////
+ //////////////////////////////////////////////////
+ // Variable where to store results
+ double ThetaLight, PhiLight, EnergyLight, ThetaTarget, EnergyTarget, ThetaHeavy, EnergyHeavy;
+
+ // Compute Kinematic using previously defined ThetaCM
+ m_InelasticBreakup.GenerateEvent(ThetaLight, PhiLight, EnergyLight, ThetaTarget, EnergyTarget, ThetaHeavy,
+ EnergyHeavy);
+ // Momentum in beam frame for light particle
+ G4ThreeVector momentum_kineLight_beam(sin(ThetaLight) * cos(phi), sin(ThetaLight) * sin(phi), cos(ThetaLight));
+ // Momentum in World frame //to go from the incident beam frame to the "world"
+ // frame
+ G4ThreeVector momentum_kineLight_world = momentum_kineLight_beam;
+ momentum_kineLight_world.rotate(Beam_theta,
+ V); // rotation of Beam_theta on Y axis
+ momentum_kineLight_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Momentum in beam frame for light particle
+ G4ThreeVector momentum_kineTarget_beam(sin(ThetaTarget) * cos(phi), sin(ThetaTarget) * sin(phi), cos(ThetaTarget));
+ // Momentum in World frame //to go from the incident beam frame to the "world"
+ // frame
+ G4ThreeVector momentum_kineTarget_world = momentum_kineTarget_beam;
+ momentum_kineTarget_world.rotate(Beam_theta,
+ V); // rotation of Beam_theta on Y axis
+ momentum_kineTarget_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Momentum in beam frame for heavy particle
+ G4ThreeVector momentum_kineHeavy_beam(sin(ThetaHeavy) * cos(phi + pi), sin(ThetaHeavy) * sin(phi + pi),
+ cos(ThetaHeavy));
+ // Momentum in World frame
+ G4ThreeVector momentum_kineHeavy_world = momentum_kineHeavy_beam;
+ momentum_kineHeavy_world.rotate(Beam_theta,
+ V); // rotation of Beam_theta on Y axis
+ momentum_kineHeavy_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Emitt secondary
+ G4DynamicParticle particleLight(LightName, momentum_kineLight_world, EnergyLight);
+ fastStep.CreateSecondaryTrack(particleLight, localPosition, time);
+
+ G4DynamicParticle particleTarget(LightName, momentum_kineTarget_world, EnergyTarget);
+ fastStep.CreateSecondaryTrack(particleTarget, localPosition, time);
+
+ G4DynamicParticle particleHeavy(HeavyName, momentum_kineHeavy_world, EnergyHeavy);
+ fastStep.CreateSecondaryTrack(particleHeavy, localPosition, time);
+
+ ///////////////////////////////////////
+ ///// Emitted particle Parameters /////
+ ///////////////////////////////////////
+ // Names 3 and 4//
+ m_ReactionConditions->SetParticleName(LightName->GetParticleName());
+ m_ReactionConditions->SetParticleName(TargetName->GetParticleName());
+ m_ReactionConditions->SetParticleName(HeavyName->GetParticleName());
+ // Angle 3 and 4 //
+ m_ReactionConditions->SetTheta(ThetaLight / deg);
+ m_ReactionConditions->SetTheta(ThetaTarget / deg);
+ m_ReactionConditions->SetTheta(ThetaHeavy / deg);
+
+ m_ReactionConditions->SetPhi(phi / deg);
+ if ((phi + pi) / deg > 360)
+ m_ReactionConditions->SetPhi((phi - pi) / deg);
+ else
+ m_ReactionConditions->SetPhi((phi + pi) / deg);
+
+ // Energy 3 and 4 //
+ m_ReactionConditions->SetKineticEnergy(EnergyLight);
+ m_ReactionConditions->SetKineticEnergy(EnergyTarget);
+ m_ReactionConditions->SetKineticEnergy(EnergyHeavy);
+ // ThetaCM and Ex//
+ m_ReactionConditions->SetThetaCM(m_InelasticBreakup.GetThetaCM() / deg);
+ m_ReactionConditions->SetExcitationEnergy3(0);
+ m_ReactionConditions->SetExcitationEnergy3(0);
+ m_ReactionConditions->SetExcitationEnergy4(0);
+ // Momuntum X 3 and 4 //
+ m_ReactionConditions->SetMomentumDirectionX(momentum_kineLight_world.x());
+ m_ReactionConditions->SetMomentumDirectionX(momentum_kineTarget_world.x());
+ m_ReactionConditions->SetMomentumDirectionX(momentum_kineHeavy_world.x());
+ // Momuntum Y 3 and 4 //
+ m_ReactionConditions->SetMomentumDirectionY(momentum_kineLight_world.y());
+ m_ReactionConditions->SetMomentumDirectionY(momentum_kineTarget_world.y());
+ m_ReactionConditions->SetMomentumDirectionY(momentum_kineHeavy_world.y());
+ // Momuntum Z 3 and 4 //
+ m_ReactionConditions->SetMomentumDirectionZ(momentum_kineLight_world.z());
+ m_ReactionConditions->SetMomentumDirectionZ(momentum_kineTarget_world.z());
+ m_ReactionConditions->SetMomentumDirectionZ(momentum_kineHeavy_world.z());
+
} // end if TwoBodyReaction
// QFS
diff --git a/NPSimulation/Process/BeamReaction.hh b/NPSimulation/Process/BeamReaction.hh
index 9d1bf2d57..b2f802f30 100644
--- a/NPSimulation/Process/BeamReaction.hh
+++ b/NPSimulation/Process/BeamReaction.hh
@@ -28,6 +28,7 @@
#include "G4AblaInterface.hh"
#include "G4Fragment.hh"
#include "G4VFastSimulationModel.hh"
+#include "NPInelasticBreakup.h"
#include "NPPhaseSpace.h"
#include "NPQFS.h"
#include "NPReaction.h"
@@ -35,7 +36,7 @@
#include "TReactionConditions.h"
class G4VPhysicalVolume;
namespace NPS {
- enum ReactionType { TwoBody, QFS, PhaseSpace, Fusion };
+ enum ReactionType { TwoBody, InelasticBreakup, QFS, PhaseSpace, Fusion };
class BeamReaction : public G4VFastSimulationModel {
public:
@@ -51,6 +52,7 @@ namespace NPS {
private:
NPL::Reaction m_Reaction;
+ NPL::InelasticBreakup m_InelasticBreakup;
NPL::QFS m_QFS;
NPL::PhaseSpace m_PhaseSpace;
string m_BeamName;
--
GitLab