diff --git a/NPLib/Physics/NPInelasticBreakUp.cxx b/NPLib/Physics/NPInelasticBreakUp.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..baad4d998e7da4d03e4da542e1774809cbf2ffbe
--- /dev/null
+++ b/NPLib/Physics/NPInelasticBreakUp.cxx
@@ -0,0 +1,659 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * *
+ * Original Author : Adrien MAtta contact: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : April 2024 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include <cmath>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include "NPCore.h"
+#include "NPFunction.h"
+#include "NPInelasticBreakup.h"
+#include "NPOptionManager.h"
+
+// Use CLHEP System of unit and Physical Constant
+#include "NPGlobalSystemOfUnits.h"
+#include "NPPhysicalConstants.h"
+#ifdef NP_SYSTEM_OF_UNITS_H
+using namespace NPUNITS;
+#endif
+
+// ROOT
+#include "TF1.h"
+
+ClassImp(InelasticBreakup)
+
+ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ InelasticBreakup::InelasticBreakup() {
+ //------------- Default Constructor -------------
+
+ //
+ fBeamEnergy = 0;
+ fThetaCM = 0;
+ fExcitation1 = 0;
+ fExcitation3 = 0;
+ fExcitation4 = 0;
+ fQValue = 0;
+ fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+ initializePrecomputeVariable();
+
+ fCrossSectionHist = NULL;
+ fExcitationEnergyHist = NULL;
+ fDoubleDifferentialCrossSectionHist = NULL;
+
+ fshoot3 = true;
+ fshoot4 = true;
+ fUseExInGeant4 = true;
+
+ fLabCrossSection = false; // flag if the provided cross-section is in the lab or not
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// This constructor aim to provide a fast way to instantiate a reaction without input file
+// The string should be of the form "A(b,c)D@E" with E the ernegy of the beam in MeV
+InelasticBreakup::InelasticBreakup(string reaction) {
+ // Instantiate the parameter to default
+ // Analyse the reaction and extract revelant information
+ string A, b, c, D, E;
+ unsigned int i = 0;
+ for (; i < reaction.length(); i++) {
+ if (reaction.compare(i, 1, "(") != 0)
+ A.push_back(reaction[i]);
+ else
+ break;
+ }
+
+ i++;
+ for (; i < reaction.length(); i++) {
+ if (reaction.compare(i, 1, ",") != 0)
+ b.push_back(reaction[i]);
+ else
+ break;
+ }
+
+ i++;
+ for (; i < reaction.length(); i++) {
+ if (reaction.compare(i, 1, ")") != 0)
+ c.push_back(reaction[i]);
+ else
+ break;
+ }
+
+ i++;
+ for (; i < reaction.length(); i++) {
+ if (reaction.compare(i, 1, "@") != 0)
+ D.push_back(reaction[i]);
+ else
+ break;
+ }
+
+ i++;
+ for (; i < reaction.length(); i++) {
+ E.push_back(reaction[i]);
+ }
+
+ fParticle1 = Beam(A);
+ fParticle2 = Particle(b);
+ fParticle3 = Particle(c);
+ fParticle4 = Particle(D);
+ fBeamEnergy = atof(E.c_str());
+ fThetaCM = 0;
+ fExcitation1 = 0;
+ fExcitation3 = 0;
+ fExcitation4 = 0;
+ fQValue = 0;
+ fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+ initializePrecomputeVariable();
+
+ // do that to avoid warning from multiple Hist with same name... int offset = 0;
+ int offset = 0;
+ while (gDirectory->FindObjectAny(Form("EnergyHist_%i", offset)) != 0)
+ ++offset;
+
+ fCrossSectionHist = new TH1F(Form("EnergyHist_%i", offset), "InelasticBreakup_CS", 1, 0, 180);
+ fDoubleDifferentialCrossSectionHist = NULL;
+
+ fshoot3 = true;
+ fshoot4 = true;
+
+ fLabCrossSection = false;
+
+ initializePrecomputeVariable();
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+InelasticBreakup::~InelasticBreakup() {}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+bool InelasticBreakup::CheckKinematic() {
+ double theta = fThetaCM;
+ /*if (m1 > m2){
+ theta = 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;
+
+ fEnergyImpulsionLab_3 = fEnergyImpulsionCM_3;
+ fEnergyImpulsionLab_3.Boost(0, 0, fBetaCM);
+ fEnergyImpulsionLab_4 = fEnergyImpulsionCM_4;
+ fEnergyImpulsionLab_4.Boost(0, 0, fBetaCM);
+
+ if (fabs(fTotalEnergyImpulsionLab.E() - (fEnergyImpulsionLab_3.E() + fEnergyImpulsionLab_4.E())) > 1e-6) {
+ cout << "Problem with energy conservation" << endl;
+ return false;
+ }
+ else {
+ // cout << "Kinematic OK" << endl;
+ return true;
+ }
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+double InelasticBreakup::ShootRandomThetaCM() {
+ double theta; // CM
+ if (fDoubleDifferentialCrossSectionHist) {
+ // Take a slice in energy
+ TAxis* Y = fDoubleDifferentialCrossSectionHist->GetYaxis();
+ int binY;
+
+ // Those test are there for the tail event of the energy distribution
+ // In case the energy is outside the range of the 2D histo we take the
+ // closest available CS
+ if (Y->FindBin(fBeamEnergy) > Y->GetLast())
+ binY = Y->GetLast() - 1;
+
+ else if (Y->FindBin(fBeamEnergy) < Y->GetFirst())
+ binY = Y->GetFirst() + 1;
+
+ else
+ binY = Y->FindBin(fBeamEnergy);
+
+ TH1D* Proj = fDoubleDifferentialCrossSectionHist->ProjectionX("proj", binY, binY);
+ SetThetaCM(theta = Proj->GetRandom() * deg);
+ }
+ else if (fLabCrossSection && fCrossSectionHist) {
+ double thetalab = -1;
+ double energylab = -1;
+ while (energylab < 0) {
+ thetalab = fCrossSectionHist->GetRandom() * deg; // shoot in lab
+ energylab = EnergyLabFromThetaLab(thetalab); // get corresponding energy
+ }
+ theta = EnergyLabToThetaCM(energylab, thetalab); // transform to theta CM
+ SetThetaCM(theta);
+ }
+ else if (fCrossSectionHist) {
+ // When root perform a Spline interpolation to shoot random number out of
+ // the distribution, it can over shoot and output a number larger that 180
+ // this lead to an additional signal at 0-4 deg Lab, especially when using a
+ // flat distribution.
+ // This fix it.
+ theta = 181;
+ if (theta / deg > 180)
+ theta = fCrossSectionHist->GetRandom();
+ // cout << " Shooting Random ThetaCM " << theta << endl;
+ SetThetaCM(theta * deg);
+ }
+ else {
+ NPL::SendErrorAndExit("NPL::InelasticBreakup", "No cross section provided, add relevant token to input file.");
+ }
+
+ return theta;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void InelasticBreakup::ShootRandomExcitationEnergy() {
+ if (fExcitationEnergyHist) {
+ SetExcitation4(fExcitationEnergyHist->GetRandom());
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void InelasticBreakup::KineRelativistic(double& ThetaLab3, double& KineticEnergyLab3, double& ThetaLab4,
+ double& KineticEnergyLab4) {
+ // 2-body relativistic kinematics: direct + inverse
+ // EnergieLab3,4 : lab energy in MeV of the 2 ejectiles
+ // ThetaLab3,4 : angles in rad
+ // case of inverse kinematics
+
+ double theta = fThetaCM;
+ if (m1 > m2) {
+ theta = 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;
+
+ fEnergyImpulsionLab_3 = fEnergyImpulsionCM_3;
+ fEnergyImpulsionLab_3.Boost(0, 0, fBetaCM);
+ fEnergyImpulsionLab_4 = fEnergyImpulsionCM_4;
+ fEnergyImpulsionLab_4.Boost(0, 0, fBetaCM);
+
+ // Angle in the lab frame
+ ThetaLab3 = fEnergyImpulsionLab_3.Angle(fEnergyImpulsionLab_1.Vect());
+ if (ThetaLab3 < 0)
+ ThetaLab3 += M_PI;
+
+ ThetaLab4 = fEnergyImpulsionLab_4.Angle(fEnergyImpulsionLab_1.Vect());
+ if (fabs(ThetaLab3) < 1e-6)
+ ThetaLab3 = 0;
+ ThetaLab4 = fabs(ThetaLab4);
+ if (fabs(ThetaLab4) < 1e-6)
+ ThetaLab4 = 0;
+
+ // Kinetic Energy in the lab frame
+ KineticEnergyLab3 = fEnergyImpulsionLab_3.E() - m3;
+ KineticEnergyLab4 = fEnergyImpulsionLab_4.E() - m4;
+
+ // test for total energy conversion
+ if (fabs(fTotalEnergyImpulsionLab.E() - (fEnergyImpulsionLab_3.E() + fEnergyImpulsionLab_4.E())) > 1e-6)
+ cout << "Problem for energy conservation" << endl;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+double InelasticBreakup::ReconstructRelativistic(double EnergyLab, double ThetaLab, double PhiLab) {
+ // EnergyLab in MeV
+ // ThetaLab in rad
+ double E3 = m3 + EnergyLab;
+ double p_Lab_3 = sqrt(E3 * E3 - m3 * m3);
+ fEnergyImpulsionLab_3 =
+ TLorentzVector(p_Lab_3 * sin(ThetaLab) * cos(PhiLab), p_Lab_3 * sin(PhiLab), p_Lab_3 * cos(ThetaLab), E3);
+ fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3;
+
+ double Eex = fEnergyImpulsionLab_4.Mag() - fParticle4.Mass();
+
+ return Eex;
+}
+
+TLorentzVector InelasticBreakup::LorentzAfterInelasticBreakup(double EnergyLab, double ThetaLab, double PhiLab) {
+ double E3 = m3 + EnergyLab;
+ double p_Lab_3 = sqrt(E3 * E3 - m3 * m3);
+ fEnergyImpulsionLab_3 =
+ TLorentzVector(p_Lab_3 * sin(ThetaLab) * cos(PhiLab), p_Lab_3 * sin(PhiLab), p_Lab_3 * cos(ThetaLab), E3);
+ fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3;
+ return fEnergyImpulsionLab_4;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Return ThetaCM
+double InelasticBreakup::EnergyLabToThetaCM(double EnergyLab, double ThetaLab) {
+ double E3 = m3 + EnergyLab;
+ double p_Lab_3 = sqrt(E3 * E3 - m3 * m3);
+
+ fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3 * sin(ThetaLab), 0, p_Lab_3 * cos(ThetaLab), E3);
+ fEnergyImpulsionCM_3 = fEnergyImpulsionLab_3;
+ fEnergyImpulsionCM_3.Boost(0, 0, -fBetaCM);
+
+ double ThetaCM = M_PI - fEnergyImpulsionCM_1.Angle(fEnergyImpulsionCM_3.Vect());
+
+ return (ThetaCM);
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Return EnergyLab
+double InelasticBreakup::EnergyLabFromThetaLab(double ThetaLab) {
+ // ThetaLab in rad
+
+ // IMPORTANT NOTICE: This function is not suitable for reaction A(c,d)B
+ // where M(c) < M(d), i.e. (p,d) or (d,3He) since in this case the same
+ // angle observed in the lab corresponds to two different energies.
+ //
+ // If this situation is encountered here:
+ // 1- the final energy will be determined randomly with a 50%-50%
+ // split for either energies.
+ // 2- Any angle outside of the allowed range (in this case) will return -1.
+ // 3- The user lab distribution will be used between angles {0,max} and the
+ // spatial distribution of the high energy and low energy branch
+ // will be the same.
+ // 4- Practically, in an experiment using a downstream spectrometer
+ // only one of the energy branches is considered.
+ //
+ // !! If both of the branches are needed the user SHOULD use the center-of-mass distribution.
+
+ //
+ // This calculation uses the formalism from J.B Marion and F.C Young
+ // (Book: Nucler InelasticBreakup Analysis, Graphs and Tables)
+
+ // Treat Exeptions
+ if (fBeamEnergy == 0)
+ return m4 * fQValue / (m3 + m4);
+
+ double A, B, C, D, ThetaLabMax = 181 * deg;
+ double Q = fQValue;
+ double T1 = fBeamEnergy;
+ double TT = fBeamEnergy + Q;
+ double sign = +1;
+
+ A = m1 * m4 * (T1 / TT) / (m1 + m2) / (m3 + m4);
+ B = m1 * m3 * (T1 / TT) / (m1 + m2) / (m3 + m4);
+ C = m2 * m3 * (1 + (m1 * Q) / (m2 * TT)) / (m1 + m2) / (m3 + m4);
+ D = m2 * m4 * (1 + (m1 * Q) / (m2 * TT)) / (m1 + m2) / (m3 + m4);
+
+ if (fabs(A + B + C + D - 1) > 1e-6 or fabs(A * C - B * D) > 1e-6) {
+ cout << " InelasticBreakup balance is wrong in NPInelasticBreakup object." << endl;
+ exit(-1);
+ }
+
+ if (B > D) {
+ ThetaLabMax = asin(sqrt(D / B));
+ if (gRandom->Rndm() < 0.5)
+ sign = -1;
+ }
+ if (ThetaLab > ThetaLabMax)
+ return -1;
+
+ double cosine = cos(ThetaLab);
+ double sine2 = pow(sin(ThetaLab), 2);
+ double factor = sqrt(D / B - sine2);
+ double EnergyLab = TT * B * pow(cosine + sign * factor, 2);
+
+ // cout << " Angle/energy: " << ThetaLab/deg << " " << EnergyLab << endl ;
+ // cin.get();
+
+ return EnergyLab;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void InelasticBreakup::Print() const {
+ // Print informations concerning the reaction
+
+ cout << "InelasticBreakup : " << fParticle2.GetName() << "(" << fParticle1.GetName() << "," << fParticle3.GetName()
+ << ")" << fParticle4.GetName() << " @ " << fBeamEnergy << " MeV" << endl;
+
+ cout << "Exc Particle 1 = " << fExcitation1 << " MeV" << endl;
+ cout << "Exc Particle 3 = " << fExcitation3 << " MeV" << endl;
+ cout << "Exc Particle 4 = " << fExcitation4 << " MeV" << endl;
+ cout << "Qgg = " << fQValue << " MeV" << endl;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+void InelasticBreakup::ReadConfigurationFile(string Path) {
+ ifstream InelasticBreakupFile;
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + Path;
+ InelasticBreakupFile.open(Path.c_str());
+ if (!InelasticBreakupFile.is_open()) {
+ InelasticBreakupFile.open(StandardPath.c_str());
+ if (InelasticBreakupFile.is_open()) {
+ Path = StandardPath;
+ }
+ else {
+ cout << "InelasticBreakup File " << Path << " not found" << endl;
+ exit(1);
+ }
+ }
+ NPL::InputParser parser(Path);
+ ReadConfigurationFile(parser);
+}
+////////////////////////////////////////////////////////////////////////////////
+Particle InelasticBreakup::GetParticle(string name, NPL::InputParser parser) {
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("DefineParticle", name);
+ unsigned int size = blocks.size();
+ if (size == 0)
+ return NPL::Particle(name);
+ else if (size == 1) {
+ cout << " -- User defined nucleus " << name << " -- " << endl;
+ vector<string> token = {"SubPart", "BindingEnergy"};
+ if (blocks[0]->HasTokenList(token)) {
+ NPL::Particle N(name, blocks[0]->GetVectorString("SubPart"), blocks[0]->GetDouble("BindingEnergy", "MeV"));
+ if (blocks[0]->HasToken("ExcitationEnergy"))
+ N.SetExcitationEnergy(blocks[0]->GetDouble("ExcitationEnergy", "MeV"));
+ if (blocks[0]->HasToken("SpinParity"))
+ N.SetSpinParity(blocks[0]->GetString("SpinParity").c_str());
+ if (blocks[0]->HasToken("Spin"))
+ N.SetSpin(blocks[0]->GetDouble("Spin", ""));
+ if (blocks[0]->HasToken("Parity"))
+ N.SetParity(blocks[0]->GetString("Parity").c_str());
+ if (blocks[0]->HasToken("LifeTime"))
+ N.SetLifeTime(blocks[0]->GetDouble("LifeTime", "s"));
+
+ cout << " -- -- -- -- -- -- -- -- -- -- --" << endl;
+ return N;
+ }
+ }
+ else {
+ NPL::SendErrorAndExit("NPL::InelasticBreakup", "Too many nuclei define with the same name");
+ }
+
+ return (NPL::Particle());
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
+void InelasticBreakup::ReadConfigurationFile(NPL::InputParser parser) {
+
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("TwoBodyInelasticBreakup");
+ if (blocks.size() > 0 && NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "\033[1;35m//// Two body reaction found " << endl;
+
+ vector<string> token1 = {"Beam", "Target", "Light", "Heavy"};
+ vector<string> token2 = {"Beam", "Target", "Particle3", "Particle4"};
+ double CSHalfOpenAngleMin = 0 * deg;
+ double CSHalfOpenAngleMax = 180 * deg;
+ for (unsigned int i = 0; i < blocks.size(); i++) {
+ if (blocks[i]->HasTokenList(token1)) {
+ int v = NPOptionManager::getInstance()->GetVerboseLevel();
+ NPOptionManager::getInstance()->SetVerboseLevel(0);
+ fParticle1.ReadConfigurationFile(parser);
+ NPOptionManager::getInstance()->SetVerboseLevel(v);
+
+ fBeamEnergy = fParticle1.GetEnergy();
+ fParticle2 = GetParticle(blocks[i]->GetString("Target"), parser);
+ fParticle3 = GetParticle(blocks[i]->GetString("Light"), parser);
+ fParticle4 = GetParticle(blocks[i]->GetString("Heavy"), parser);
+ }
+ else if (blocks[i]->HasTokenList(token2)) {
+ fParticle1.SetVerboseLevel(0);
+ fParticle1.ReadConfigurationFile(parser);
+ fBeamEnergy = fParticle1.GetEnergy();
+
+ fParticle2 = GetParticle(blocks[i]->GetString("Target"), parser);
+ fParticle3 = GetParticle(blocks[i]->GetString("Particle3"), parser);
+ fParticle4 = GetParticle(blocks[i]->GetString("Particle4"), parser);
+ }
+ else {
+ cout << "ERROR: check your input file formatting \033[0m" << endl;
+ exit(1);
+ }
+
+ if (blocks[i]->HasToken("ExcitationEnergyBeam")) {
+ fExcitation1 = blocks[i]->GetDouble("ExcitationEnergyBeam", "MeV");
+ }
+ else if (blocks[i]->HasToken("ExcitationEnergy1")) {
+ fExcitation1 = blocks[i]->GetDouble("ExcitationEnergy1", "MeV");
+ }
+
+ if (blocks[i]->HasToken("ExcitationEnergyLight"))
+ fExcitation3 = blocks[i]->GetDouble("ExcitationEnergyLight", "MeV");
+ else if (blocks[i]->HasToken("ExcitationEnergy3"))
+ fExcitation3 = blocks[i]->GetDouble("ExcitationEnergy3", "MeV");
+
+ if (blocks[i]->HasToken("ExcitationEnergyHeavy"))
+ fExcitation4 = blocks[i]->GetDouble("ExcitationEnergyHeavy", "MeV");
+ else if (blocks[i]->HasToken("ExcitationEnergy4"))
+ fExcitation4 = blocks[i]->GetDouble("ExcitationEnergy4", "MeV");
+
+ if (blocks[i]->HasToken("ExcitationEnergyDistribution")) {
+ vector<string> file = blocks[i]->GetVectorString("ExcitationEnergyDistribution");
+ fExcitationEnergyHist = Read1DProfile(file[0], file[1]);
+ fExcitation4 = 0;
+ }
+
+ if (blocks[i]->HasToken("CrossSectionPath")) {
+ vector<string> file = blocks[i]->GetVectorString("CrossSectionPath");
+ TH1F* CStemp = Read1DProfile(file[0], file[1]);
+
+ // multiply CStemp by sin(theta)
+ TF1* fsin = new TF1("sin", Form("1/(sin(x*%f/180.))", M_PI), 0, 180);
+ CStemp->Divide(fsin, 1);
+ SetCrossSectionHist(CStemp);
+ delete fsin;
+ }
+
+ if (blocks[i]->HasToken("LabCrossSectionPath")) {
+ fLabCrossSection = true;
+
+ vector<string> file = blocks[i]->GetVectorString("LabCrossSectionPath");
+ TH1F* CStemp = Read1DProfile(file[0], file[1]);
+
+ // multiply CStemp by sin(theta)
+ TF1* fsin = new TF1("sin", Form("1/(sin(x*%f/180.))", M_PI), 0, 180);
+ CStemp->Divide(fsin, 1);
+ SetCrossSectionHist(CStemp);
+ delete fsin;
+ }
+
+ if (blocks[i]->HasToken("DoubleDifferentialCrossSectionPath")) {
+ vector<string> file = blocks[i]->GetVectorString("DoubleDifferentialCrossSectionPath");
+ TH2F* CStemp = Read2DProfile(file[0], file[1]);
+
+ // multiply CStemp by sin(theta)
+ // X axis is theta CM
+ // Y axis is beam energy
+ // Division affect only X axis
+ TF1* fsin = new TF1("sin", Form("1/(sin(x*%f/180.))", M_PI), 0, 180);
+ CStemp->Divide(fsin, 1);
+
+ SetDoubleDifferentialCrossSectionHist(CStemp);
+ delete fsin;
+ }
+
+ if (blocks[i]->HasToken("HalfOpenAngleMin")) {
+ CSHalfOpenAngleMin = blocks[i]->GetDouble("HalfOpenAngleMin", "deg");
+ }
+ if (blocks[i]->HasToken("HalfOpenAngleMax")) {
+ CSHalfOpenAngleMax = blocks[i]->GetDouble("HalfOpenAngleMax", "deg");
+ }
+ if (blocks[i]->HasToken("Shoot3")) {
+ fshoot3 = blocks[i]->GetInt("Shoot3");
+ }
+ if (blocks[i]->HasToken("Shoot4")) {
+ fshoot4 = blocks[i]->GetInt("Shoot4");
+ }
+ if (blocks[i]->HasToken("ShootHeavy")) {
+ fshoot4 = blocks[i]->GetInt("ShootHeavy");
+ }
+ if (blocks[i]->HasToken("ShootLight")) {
+ fshoot3 = blocks[i]->GetInt("ShootLight");
+ }
+ if (blocks[i]->HasToken("UseExInGeant4")) {
+ // This option will not change the Ex of the produced ion in G4 Tracking
+ // This is to be set to true when using a Ex distribution without decay
+ // Otherwise the Ion Table size grew four ech event slowing down the
+ // simulation
+ fUseExInGeant4 = blocks[i]->GetInt("UseExInGeant4");
+ }
+ }
+ SetCSAngle(CSHalfOpenAngleMin / deg, CSHalfOpenAngleMax / deg);
+ initializePrecomputeVariable();
+ cout << "\033[0m";
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void InelasticBreakup::initializePrecomputeVariable() {
+
+ if (fBeamEnergy < 0)
+ fBeamEnergy = 0;
+
+ if (fExcitation1 >= 0)
+ fParticle1.SetExcitationEnergy(fExcitation1); // Write over the beam excitation energy
+
+ // fParticle1.GetExcitationEnergy() is a copy of fExcitation1
+ m1 = fParticle1.Mass() + fParticle1.GetExcitationEnergy(); // in case of isomeric state,
+ m2 = fParticle2.Mass(); // Target
+ m3 = fParticle3.Mass() + fExcitation3;
+ m4 = fParticle4.Mass() + fExcitation4;
+ fQValue = m1 + m2 - m3 - m4;
+
+ s = m1 * m1 + m2 * m2 + 2 * m2 * (fBeamEnergy + m1);
+ fTotalEnergyImpulsionCM = TLorentzVector(0, 0, 0, sqrt(s));
+ fEcm = sqrt(s) - m1 - m2;
+
+ ECM_1 = (s + m1 * m1 - m2 * m2) / (2 * sqrt(s));
+ ECM_2 = (s + m2 * m2 - m1 * m1) / (2 * sqrt(s));
+ ECM_3 = (s + m3 * m3 - m4 * m4) / (2 * sqrt(s));
+ ECM_4 = (s + m4 * m4 - m3 * m3) / (2 * sqrt(s));
+
+ pCM_1 = sqrt(ECM_1 * ECM_1 - m1 * m1);
+ pCM_2 = sqrt(ECM_2 * ECM_2 - m2 * m2);
+ pCM_3 = sqrt(ECM_3 * ECM_3 - m3 * m3);
+ pCM_4 = sqrt(ECM_4 * ECM_4 - m4 * m4);
+
+ fImpulsionLab_1 = TVector3(0, 0, sqrt(fBeamEnergy * fBeamEnergy + 2 * fBeamEnergy * m1));
+ fImpulsionLab_2 = TVector3(0, 0, 0);
+
+ fEnergyImpulsionLab_1 = TLorentzVector(fImpulsionLab_1, m1 + fBeamEnergy);
+ fEnergyImpulsionLab_2 = TLorentzVector(fImpulsionLab_2, m2);
+
+ fTotalEnergyImpulsionLab = fEnergyImpulsionLab_1 + fEnergyImpulsionLab_2;
+
+ fBetaCM = fTotalEnergyImpulsionLab.Beta();
+
+ fEnergyImpulsionCM_1 = fEnergyImpulsionLab_1;
+ fEnergyImpulsionCM_1.Boost(0, 0, -fBetaCM);
+
+ fEnergyImpulsionCM_2 = fEnergyImpulsionLab_2;
+ fEnergyImpulsionCM_2.Boost(0, 0, -fBetaCM);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void InelasticBreakup::SetParticle3(double EnergyLab, double ThetaLab) {
+ double p3 = sqrt(pow(EnergyLab, 2) + 2 * m3 * EnergyLab);
+
+ fEnergyImpulsionLab_3 = TLorentzVector(p3 * sin(ThetaLab), 0, p3 * cos(ThetaLab), EnergyLab + m3);
+ fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3;
+
+ fParticle3.SetEnergyImpulsion(fEnergyImpulsionLab_3);
+ fParticle4.SetEnergyImpulsion(fEnergyImpulsionLab_4);
+
+ fThetaCM = EnergyLabToThetaCM(EnergyLab, ThetaLab);
+ fExcitation4 = ReconstructRelativistic(EnergyLab, ThetaLab);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+Double_t InelasticBreakup::GetTotalCrossSection() const {
+ Double_t stot = fCrossSectionHist->Integral("width"); // take bin width into account (in deg!)
+ stot *= M_PI / 180; // correct so that bin width is in rad
+ stot *= 2 * M_PI; // integration over phi
+
+ return stot;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void InelasticBreakup::SetCSAngle(double CSHalfOpenAngleMin, double CSHalfOpenAngleMax) {
+ if (fCrossSectionHist) {
+ for (int i = 0; i < fCrossSectionHist->GetNbinsX(); i++) {
+ if (fCrossSectionHist->GetBinCenter(i) > CSHalfOpenAngleMax ||
+ fCrossSectionHist->GetBinCenter(i) < CSHalfOpenAngleMin) {
+ fCrossSectionHist->SetBinContent(i, 0);
+ }
+ }
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Check whenever the reaction is allowed at the given energy
+bool InelasticBreakup::IsAllowed(double Energy) {
+ double AvailableEnergy = Energy + fParticle1.Mass() + fParticle2.Mass();
+ double RequiredEnergy = fParticle3.Mass() + fParticle4.Mass();
+
+ if (AvailableEnergy > RequiredEnergy)
+ return true;
+ else
+ return false;
+}
diff --git a/NPLib/Physics/NPInelasticBreakUp.h b/NPLib/Physics/NPInelasticBreakUp.h
new file mode 100644
index 0000000000000000000000000000000000000000..1a3dd89014ec1af30f75948cb3d4c93d80759fd5
--- /dev/null
+++ b/NPLib/Physics/NPInelasticBreakUp.h
@@ -0,0 +1,265 @@
+#ifndef NPInelasticBreakUp_h
+#define NPInelasticBreakUp_h
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * *
+ * Original Author : Adrien MAtta contact: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : April 2024 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+// C++ header
+#include <string>
+
+// NPL
+#include "NPBeam.h"
+#include "NPInputParser.h"
+#include "NPParticle.h"
+using namespace NPL;
+
+// ROOT header
+#include "NPReaction.h"
+#include "TCanvas.h"
+#include "TGraph.h"
+#include "TH1F.h"
+#include "TLorentzRotation.h"
+#include "TLorentzVector.h"
+#include "TRandom.h"
+#include "TVector3.h"
+
+using namespace std;
+
+namespace NPL {
+
+ class InelasticBreakup {
+
+ public: // Constructors and Destructors
+ InelasticBreakup();
+ ~InelasticBreakup();
+
+ public: // Various Method
+ Particle GetParticle(string name, NPL::InputParser parser);
+ void ReadConfigurationFile(string Path);
+ void ReadConfigurationFile(NPL::InputParser);
+
+ private:
+ TRandom* RandGen;
+
+ private:
+ int fVerboseLevel;
+
+ private: // use for Monte Carlo simulation
+ bool fshoot3;
+ bool fshoot4;
+ bool fshoot5;
+ bool fUseExInGeant4;
+ bool fLabCrossSection;
+
+ private:
+ Beam fParticle1; // Beam
+ Particle fParticle2; // Target
+ Particle fParticle3; // Light ejectile
+ Particle fParticle4; // Heavy ejectile
+ double fQValue; // Q-value in MeV
+ double fEcm; // Ecm in MeV
+ double fBeamEnergy; // Beam energy in MeV
+ double fThetaCM; // Center-of-mass angle in radian
+ double fExcitation1; // Excitation energy in MeV of the beam, useful for isomers
+ double fExcitation3; // Excitation energy in MeV
+ double fExcitation4; // Excitation energy in MeV
+ TH1F* fCrossSectionHist; // Differential cross section in CM frame
+ TH2F* fDoubleDifferentialCrossSectionHist; // Diff. CS CM frame vs Beam E
+ TH1F* fExcitationEnergyHist; // Distribution of Excitation energy
+ public:
+ // Getters and Setters
+ void SetBeamEnergy(const double& eBeam) {
+ fBeamEnergy = eBeam;
+ initializePrecomputeVariable();
+ }
+ void SetEcm(const double& Ecm) {
+ fEcm = Ecm;
+ s = pow(Ecm + m1 + m2, 2);
+ fBeamEnergy = (s - m1 * m1 - m2 * m2) / (2 * m2) - m1;
+ initializePrecomputeVariable();
+ }
+ void SetThetaCM(const double& angle) {
+ fThetaCM = angle;
+ initializePrecomputeVariable();
+ }
+ void SetExcitation1(const double& exci) {
+ fExcitation1 = exci;
+ initializePrecomputeVariable();
+ }
+ void SetExcitation3(const double& exci) {
+ fExcitation3 = exci;
+ initializePrecomputeVariable();
+ }
+ void SetExcitation4(const double& exci) {
+ fExcitation4 = exci;
+ initializePrecomputeVariable();
+ }
+ // For retro compatibility
+ void SetExcitationBeam(const double& exci) {
+ fExcitation1 = exci;
+ initializePrecomputeVariable();
+ }
+ void SetExcitationLight(const double& exci) {
+ fExcitation3 = exci;
+ initializePrecomputeVariable();
+ }
+ void SetExcitationHeavy(const double& exci) {
+ fExcitation4 = exci;
+ initializePrecomputeVariable();
+ }
+ void SetVerboseLevel(const int& verbose) { fVerboseLevel = verbose; }
+ void SetCrossSectionHist(TH1F* CrossSectionHist) {
+ delete fCrossSectionHist;
+ fCrossSectionHist = CrossSectionHist;
+ }
+
+ void SetDoubleDifferentialCrossSectionHist(TH2F* CrossSectionHist) {
+ fDoubleDifferentialCrossSectionHist = CrossSectionHist;
+ }
+ double GetBeamEnergy() const { return fBeamEnergy; }
+ double GetThetaCM() const { return fThetaCM; }
+ double GetExcitation1() const { return fExcitation1; }
+ double GetExcitation3() const { return fExcitation3; }
+ double GetExcitation4() const { return fExcitation4; }
+ double GetQValue() const { return fQValue; }
+ double GetEcm() const { return fEcm; }
+ Particle* GetParticle1() { return &fParticle1; }
+ Particle* GetParticle2() { return &fParticle2; }
+ Particle* GetParticle3() { return &fParticle3; }
+ Particle* GetParticle4() { return &fParticle4; }
+ Particle* GetNucleus1() { return GetParticle1(); }
+ Particle* GetNucleus2() { return GetParticle2(); }
+ Particle* GetNucleus3() { return GetParticle3(); }
+ Particle* GetNucleus4() { return GetParticle4(); }
+
+ TH1F* GetCrossSectionHist() const { return fCrossSectionHist; }
+ int GetVerboseLevel() const { return fVerboseLevel; }
+ bool GetShoot3() const { return fshoot3; }
+ bool GetShoot4() const { return fshoot4; }
+ bool GetUseExInGeant4() const { return fUseExInGeant4; }
+
+ public:
+ // Modify the CS histo so cross section shoot is within ]HalfOpenAngleMin,HalfOpenAngleMax[
+ void SetCSAngle(double CSHalfOpenAngleMin, double CSHalfOpenAngleMax);
+
+ private: // intern precompute variable
+ void initializePrecomputeVariable();
+ double m1;
+ double m2;
+ double m3;
+ double m4;
+ double m5;
+
+ // Lorents Vector
+ TLorentzVector fEnergyImpulsionLab_1;
+ TLorentzVector fEnergyImpulsionLab_2;
+ TLorentzVector fEnergyImpulsionLab_3;
+ TLorentzVector fEnergyImpulsionLab_4;
+ TLorentzVector fTotalEnergyImpulsionLab;
+
+ TLorentzVector fEnergyImpulsionCM_1;
+ TLorentzVector fEnergyImpulsionCM_2;
+ TLorentzVector fEnergyImpulsionCM_3;
+ TLorentzVector fEnergyImpulsionCM_4;
+ TLorentzVector fTotalEnergyImpulsionCM;
+
+ // Impulsion Vector3
+ TVector3 fImpulsionLab_1;
+ TVector3 fImpulsionLab_2;
+ TVector3 fImpulsionLab_3;
+ TVector3 fImpulsionLab_4;
+
+ TVector3 fImpulsionCM_1;
+ TVector3 fImpulsionCM_2;
+ TVector3 fImpulsionCM_3;
+ TVector3 fImpulsionCM_4;
+
+ // CM Energy composante & CM impulsion norme
+ Double_t ECM_1;
+ Double_t ECM_2;
+ Double_t ECM_3;
+ Double_t ECM_4;
+ Double_t pCM_1;
+ Double_t pCM_2;
+ Double_t pCM_3;
+ Double_t pCM_4;
+
+ // Mandelstam variable
+ Double_t s;
+
+ // Center of Mass Kinematic
+ Double_t fBetaCM;
+
+ public:
+ TLorentzVector GetEnergyImpulsionLab_1() const { return fEnergyImpulsionLab_1; }
+ TLorentzVector GetEnergyImpulsionLab_2() const { return fEnergyImpulsionLab_2; }
+ TLorentzVector GetEnergyImpulsionLab_3() const { return fEnergyImpulsionLab_3; }
+ TLorentzVector GetEnergyImpulsionLab_4() const { return fEnergyImpulsionLab_4; }
+
+ public: // Kinematics
+ // Check that the reaction is alowed
+ bool CheckKinematic();
+ bool IsLabCrossSection() { return fLabCrossSection; };
+
+ // Use fCrossSectionHist to shoot a Random ThetaCM and set fThetaCM to this value
+ double ShootRandomThetaCM();
+
+ // Use fExcitationEnergyHist to shoot a Random Excitation energy and set fExcitation4 to this value
+ void ShootRandomExcitationEnergy();
+
+ // Compute ThetaLab and EnergyLab for product of reaction
+ void KineRelativistic(double& ThetaLab3, double& KineticEnergyLab3, double& ThetaLab4, double& KineticEnergyLab4);
+
+ // Return Excitation Energy
+ double ReconstructRelativistic(double EnergyLab, double ThetaLab, double PhiLab = 0);
+
+ // Return Lorentz Vector of Heavy Recoil after reaction
+ TLorentzVector LorentzAfterInelasticBreakup(double EnergyLab, double ThetaLab, double PhiLab = 0);
+
+ // Return ThetaCM
+ // EnergyLab: energy measured in the laboratory frame
+ // ExcitationEnergy: excitation energy previously calculated.
+ double EnergyLabToThetaCM(double EnergyLab, double ThetaLab);
+ // Return theoretical EnergyLab, useful for a random distribution in the lab frame
+ // ThetaLab: angle measured in the laboratory frame
+ // This uses the fExcitation4 and fQValue both set previously
+ double EnergyLabFromThetaLab(double ThetaLab);
+
+ // Check whenever the reaction is allowed at the given energy
+ bool IsAllowed(double Energy);
+
+ void SetParticle3(double EnergyLab, double ThetaLab);
+
+ double GetP_CM_1() { return pCM_1; }
+ double GetP_CM_2() { return pCM_2; }
+ double GetP_CM_3() { return pCM_3; }
+ double GetP_CM_4() { return pCM_4; }
+ double GetE_CM_1() { return ECM_1; }
+ double GetE_CM_2() { return ECM_2; }
+ double GetE_CM_3() { return ECM_3; }
+ double GetE_CM_4() { return ECM_4; }
+
+ // Print private paremeter
+ void Print() const;
+
+ ClassDef(InelasticBreakup, 0)
+ };
+} // namespace NPL
+#endif