From 81e52298b34b8ca91eb2b30f37240d065256b8ac Mon Sep 17 00:00:00 2001
From: Adrien Matta <matta@lpccaen.in2p3.fr>
Date: Wed, 23 Feb 2022 17:24:16 +0100
Subject: [PATCH] * Progress on Phase space
---
NPLib/Physics/NPPhaseSpace.cxx | 36 ++++++++++++++++++++----
NPLib/Physics/NPPhaseSpace.h | 19 +++++++++----
NPSimulation/Process/BeamReaction.cc | 41 ++++++++++++++++++++++------
3 files changed, 77 insertions(+), 19 deletions(-)
diff --git a/NPLib/Physics/NPPhaseSpace.cxx b/NPLib/Physics/NPPhaseSpace.cxx
index 9599230fc..02d1cd99f 100644
--- a/NPLib/Physics/NPPhaseSpace.cxx
+++ b/NPLib/Physics/NPPhaseSpace.cxx
@@ -53,30 +53,34 @@ void NPL::PhaseSpace::ReadConfigurationFile(NPL::InputParser parser){
vector<string> token = {"Beam","Target","Daughters","ExcitationEnergies","Fermi"};
fDaughters.clear();
-
+ fMasses.clear();
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(token)){
int v = NPOptionManager::getInstance()->GetVerboseLevel();
NPOptionManager::getInstance()->SetVerboseLevel(0);
fBeam.ReadConfigurationFile(parser);
NPOptionManager::getInstance()->SetVerboseLevel(v);
-
fBeamEnergy= fBeam.GetEnergy();
fTarget= GetParticle(blocks[i]->GetString("Target"),parser);
+ fTargetLV=TLorentzVector(0,0,0,fTarget.Mass());
vector<string> vDaughters= blocks[i]->GetVectorString("Daughters");
fExcitations= blocks[i]->GetVectorDouble("ExcitationEnergies","MeV");
-
- for(auto d : vDaughters){
+ for(auto d:vDaughters){
fDaughters.push_back(GetParticle(d,parser));
+ fMasses.push_back(fDaughters.rend()->Mass());
}
+
+ SetBeamLV(0,0,1,fBeam.GetGamma()*fBeam.Mass());
+
}
else{
cout << "ERROR: check your input file formatting \033[0m" << endl;
exit(1);
}
+
+ }
}
-}
////////////////////////////////////////////////////////////////////////////////
Particle NPL::PhaseSpace::GetParticle(string name, NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("DefineParticle",name);
@@ -109,3 +113,25 @@ Particle NPL::PhaseSpace::GetParticle(string name, NPL::InputParser parser){
return (NPL::Particle());
}
+////////////////////////////////////////////////////////////////////////////////
+bool NPL::PhaseSpace::SetBeamLV(const TLorentzVector& LV){
+ fBeamLV=LV;
+ fTotalLV = fBeamLV+fTargetLV;
+ static std::string option;
+ option="";
+ if(fFermi)
+ option="Fermi";
+ return fPhaseSpace.SetDecay(fTotalLV,fMasses.size(),&fMasses[0],option.c_str());
+
+ }
+////////////////////////////////////////////////////////////////////////////////
+bool NPL::PhaseSpace::SetBeamLV(double dirX,double dirY,double dirZ,double Kinetic){
+ fBeam.SetKineticEnergy(Kinetic);
+ double E = fBeam.GetGamma()*fBeam.Mass();
+ double pc=sqrt(E*E-fBeam.Mass()*fBeam.Mass());
+ TVector3 p(dirX,dirY,dirZ);
+ p.Unit();
+ p*=pc;
+ TLorentzVector beamLV(p,E);
+ return SetBeamLV(fBeamLV);
+ }
diff --git a/NPLib/Physics/NPPhaseSpace.h b/NPLib/Physics/NPPhaseSpace.h
index 81447eade..737361c23 100644
--- a/NPLib/Physics/NPPhaseSpace.h
+++ b/NPLib/Physics/NPPhaseSpace.h
@@ -63,13 +63,22 @@ namespace NPL{
public:
// Getters and Setters
- void SetBeamEnergy(const double& eBeam) {fBeamEnergy = eBeam; initializePrecomputeVariable();}
+ bool SetBeamLV(const TLorentzVector& LV);
+ bool SetBeamLV(double dirX,double dirY,double dirZ,double Kinetic);
Particle* GetBeam(){return &fBeam;};
-
+ double Generate(){return fPhaseSpace.Generate();}
+ TLorentzVector* GetDecayLV(unsigned int& n){return fPhaseSpace.GetDecay(n);}
+ size_t GetDecaySize(){return fDaughters.size();}
+ Particle* GetParticle(unsigned int& n){return &fDaughters[n];}
+ double GetExcitation(unsigned int& n){return fExcitations[n];}
+
private: // intern precompute variable
- void initializePrecomputeVariable();
- TLorentzVector fInitialEnergyImpulsion;
- std::vector<double> fmasses;
+ void Init();
+ TLorentzVector fBeamLV;
+ TLorentzVector fTargetLV;
+ TLorentzVector fTotalLV;
+
+ std::vector<double> fMasses;
ClassDef(PhaseSpace,0)
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 509f6955f..a09f5fccc 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -197,7 +197,7 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
if (m_shoot && m_length < m_StepSize) {
if(m_ReactionType==QFS){
//if ( m_QFS.IsAllowed() ) {
- return true;
+ return true;
//}
//else{
// m_shoot=false;
@@ -505,16 +505,16 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
int j=0;
m_QFS.SetIsAllowed(false);
while(!m_QFS.IsAllowed()){
- // Shoot internal momentum for the removed cluster
- // if a momentum Sigma is given then shoot in 3 indep. Gaussians
- // if input files are given for distributions use them instead
- m_QFS.ShootInternalMomentum();
+ // Shoot internal momentum for the removed cluster
+ // if a momentum Sigma is given then shoot in 3 indep. Gaussians
+ // if input files are given for distributions use them instead
+ m_QFS.ShootInternalMomentum();
- // Go from CM to Lab
- m_QFS.KineRelativistic(Theta1, Phi1, TKE1, Theta2, Phi2, TKE2);
+ // Go from CM to Lab
+ m_QFS.KineRelativistic(Theta1, Phi1, TKE1, Theta2, Phi2, TKE2);
- j++;
- if(j>100) cout<< "ERROR: too many iteration and QFS kinematical conditions not allowed"<<endl;
+ j++;
+ if(j>100) cout<< "ERROR: too many iteration and QFS kinematical conditions not allowed"<<endl;
}
//---------------------------------------------------------
@@ -666,6 +666,29 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
}// end QFS
+ //////////////////////////
+ // Phase space Case //
+ //////////////////////////
+ else if(m_ReactionType==PhaseSpace){
+ // Prepare beam LV
+ if(m_PhaseSpace.SetBeamLV(pdirection.x(),pdirection.y(),pdirection.z(),energy)){
+ for(unsigned int i = 0,size=m_PhaseSpace.GetDecaySize() ; i < size ; i++){
+ int d_Z = m_QFS.GetParticle1()->GetZ();
+ int d_A = m_QFS.GetParticle1()->GetA();
+
+ static G4IonTable* IonTable
+ = G4ParticleTable::GetParticleTable()->GetIonTable();
+
+ G4ParticleDefinition* dName;
+
+ if (d_Z == 0 && d_A == 1) // neutron is special case
+ dName = G4Neutron::Definition();
+ else
+ dName = IonTable->GetIon(d_Z, d_A);
+ }
+ }
+ }
+
////////////////////
// Fusion Case //
////////////////////
--
GitLab