From 02f30ca2087b57dd845f35f01bc24846f606688e Mon Sep 17 00:00:00 2001
From: Adrien Matta <matta@lpccaen.in2p3.fr>
Date: Tue, 1 Oct 2024 09:03:15 +0200
Subject: [PATCH] * adding option to prevent direct reaction to occure in
inelastic breakup
---
NPLib/Physics/NPInelasticBreakup.cxx | 1 +
NPLib/Physics/NPReaction.cxx | 25 ++++++++++++-------------
NPLib/Physics/NPReaction.h | 2 ++
3 files changed, 15 insertions(+), 13 deletions(-)
diff --git a/NPLib/Physics/NPInelasticBreakup.cxx b/NPLib/Physics/NPInelasticBreakup.cxx
index 1714ed6d6..30ca3bae5 100644
--- a/NPLib/Physics/NPInelasticBreakup.cxx
+++ b/NPLib/Physics/NPInelasticBreakup.cxx
@@ -190,6 +190,7 @@ void InelasticBreakup::ReadConfigurationFile(NPL::InputParser parser) {
fReaction.SetParticle3(GetParticle(blocks[i]->GetString("Target"), parser));
fReaction.SetParticle4(GetParticle(blocks[i]->GetString("Heavy"), parser));
fReaction.SetBeamEnergy(fBeamEnergy);
+ fReaction.PreventDirect();
fReaction.initializePrecomputeVariable();
fParticle3 = GetParticle(blocks[i]->GetString("Light"), parser);
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index e79397f08..baa83d5c1 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -87,7 +87,7 @@ ClassImp(Reaction)
fshoot3 = true;
fshoot4 = true;
fUseExInGeant4 = true;
-
+ fallow_direct = true;
fLabCrossSection = false; // flag if the provided cross-section is in the lab or not
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -261,9 +261,9 @@ void Reaction::KineRelativistic(double& ThetaLab3, double& KineticEnergyLab3, do
// case of inverse kinematics
double theta = fThetaCM;
- if (m1 > m2){
+ if (m1 > m2 && fallow_direct) {
theta = M_PI - fThetaCM;
- //fThetaCM = M_PI - fThetaCM;
+ // fThetaCM = M_PI - fThetaCM;
}
fEnergyImpulsionCM_3 = TLorentzVector(pCM_3 * sin(theta), 0, pCM_3 * cos(theta), ECM_3);
fEnergyImpulsionCM_4 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_3;
@@ -309,7 +309,7 @@ double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab, doub
return Eex;
}
-TLorentzVector Reaction::LorentzAfterReaction(double EnergyLab, double ThetaLab, double PhiLab){
+TLorentzVector Reaction::LorentzAfterReaction(double EnergyLab, double ThetaLab, double PhiLab) {
double E3 = m3 + EnergyLab;
double p_Lab_3 = sqrt(E3 * E3 - m3 * m3);
fEnergyImpulsionLab_3 =
@@ -753,20 +753,20 @@ TGraph* Reaction::GetThetaLabVersusThetaCM(double AngleStep_CM) {
return (fAngleLine);
}
////////////////////////////////////////////////////////////////////////////////////////////
-TGraph* Reaction::GetJacobian(double AngleStep_CM){
+TGraph* Reaction::GetJacobian(double AngleStep_CM) {
vector<double> vx;
vector<double> vy;
double theta3, E3, theta4, E4;
double E3_CM = GetE_CM_3();
double P3_CM = GetP_CM_3();
- double B = P3_CM/E3_CM;
+ double B = P3_CM / E3_CM;
double Mass1 = fParticle1.Mass();
double Mass2 = fParticle2.Mass();
double BeamEnergy = GetBeamEnergy();
- double beta = sqrt(BeamEnergy*BeamEnergy+2*Mass1*BeamEnergy)/(BeamEnergy+Mass1+Mass2);
- double r = beta/B;
- double gamma = 1./sqrt(1-beta*beta);
+ double beta = sqrt(BeamEnergy * BeamEnergy + 2 * Mass1 * BeamEnergy) / (BeamEnergy + Mass1 + Mass2);
+ double r = beta / B;
+ double gamma = 1. / sqrt(1 - beta * beta);
double Jacobian_num;
double Jacobian_denum;
@@ -775,10 +775,10 @@ TGraph* Reaction::GetJacobian(double AngleStep_CM){
SetThetaCM(angle * deg);
KineRelativistic(theta3, E3, theta4, E4);
- Jacobian_num = abs(gamma*(1.+r*cos(fThetaCM)));
+ Jacobian_num = abs(gamma * (1. + r * cos(fThetaCM)));
- Jacobian_denum = pow(gamma,2)*pow(r+cos(fThetaCM),2) + pow(sin(fThetaCM),2);
- Jacobian_denum = pow(Jacobian_denum,3./2.);
+ Jacobian_denum = pow(gamma, 2) * pow(r + cos(fThetaCM), 2) + pow(sin(fThetaCM), 2);
+ Jacobian_denum = pow(Jacobian_denum, 3. / 2.);
Jacobian = Jacobian_num / Jacobian_denum;
@@ -789,7 +789,6 @@ TGraph* Reaction::GetJacobian(double AngleStep_CM){
TGraph* gJaco = new TGraph(vx.size(), &vx[0], &vy[0]);
return gJaco;
-
}
////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 62926e9ce..2c801a048 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -75,6 +75,7 @@ namespace NPL {
bool fshoot4;
bool fUseExInGeant4;
bool fLabCrossSection;
+ bool fallow_direct;
private: // use to display the kinematical line
TGraph* fKineLine3;
@@ -169,6 +170,7 @@ namespace NPL {
bool GetShoot3() const { return fshoot3; }
bool GetShoot4() const { return fshoot4; }
bool GetUseExInGeant4() const { return fUseExInGeant4; }
+ void PreventDirect() { fallow_direct = false; }
public:
// Modify the CS histo so cross section shoot is within ]HalfOpenAngleMin,HalfOpenAngleMax[
--
GitLab