diff --git a/NPLib/Physics/NPInelasticBreakup.cxx b/NPLib/Physics/NPInelasticBreakup.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..97eb07b64a2022c4086ff5bc5aaf0bcbb2b275a6
--- /dev/null
+++ b/NPLib/Physics/NPInelasticBreakup.cxx
@@ -0,0 +1,248 @@
+/*****************************************************************************
+ * 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;
+ fQValue = 0;
+ fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+ initializePrecomputeVariable();
+
+ fCrossSectionHist = NULL;
+ fDoubleDifferentialCrossSectionHist = NULL;
+
+ fLabCrossSection = false; // flag if the provided cross-section is in the lab or not
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+InelasticBreakup::~InelasticBreakup() {}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void InelasticBreakup::GenerateEvent(double& ThetaLight, double& KineticEnergyLight, 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 (fParticle1.Mass() > fParticle2.Mass()) {
+ theta = M_PI - fThetaCM;
+ // fThetaCM = M_PI - fThetaCM;
+ }
+
+ // shoot light first
+ double thetaLight = gRandom->Gaus(fMeanAngleLight, fSigmaAngleLight);
+ double energyLight = gRandom->Gaus(fMeanEnergyLight, fSigmaEnergyLight);
+ double phiLight = gRandom->Uniform() * 2. * M_PI;
+ double totalELight = sqrt(fParticle3.Mass() * fParticle3.Mass() + energyLight * energyLight);
+ TLorentzVector LVLight(energyLight * sin(thetaLight) * cos(phiLight), energyLight * sin(thetaLight) * sin(phiLight),
+ energyLight * cos(thetaLight), totalELight);
+
+ // Substract light LV to Beam
+
+ // Two body kine elastic scattering with heavy
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+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", "MeanEnergyLight", "SigmaEnergyLight", "CrossSectionPath"};
+ 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);
+ // Read the beam block
+ fBeam.ReadConfigurationFile(parser);
+ NPOptionManager::getInstance()->SetVerboseLevel(v);
+
+ fBeamEnergy = fBeam.GetEnergy();
+ // set the particle
+ fReaction.SetParticle2(GetParticle(blocks[i]->GetString("Target"), parser));
+ fReaction.SetParticle3(GetParticle(blocks[i]->GetString("Target"), parser));
+ fReaction.SetParticle4(GetParticle(blocks[i]->GetString("Heavy"), parser));
+ fMeanEnergyLight = blocks[i]->GetDouble("MeanEnergyLight", "MeV");
+ fSigmaEnergyLight = blocks[i]->GetDouble("SigmaEnergyLight", "MeV");
+ fMeanAngleLight = blocks[i]->GetDouble("MeanAngleLight", "deg");
+ fSigmaAngleLight = blocks[i]->GetDouble("SigmaAngleLight", "deg");
+ }
+ else {
+ cout << "ERROR: check your input file formatting \033[0m" << endl;
+ exit(1);
+ }
+
+ if (blocks[i]->HasToken("ExcitationEnergyLight"))
+ fExcitationLight = blocks[i]->GetDouble("ExcitationEnergyLight", "MeV");
+
+ if (blocks[i]->HasToken("ExcitationEnergyHeavy"))
+ fExcitationHeavy = blocks[i]->GetDouble("ExcitationEnergyHeavy", "MeV");
+
+ 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");
+ }
+ }
+ SetCSAngle(CSHalfOpenAngleMin / deg, CSHalfOpenAngleMax / deg);
+ initializePrecomputeVariable();
+ cout << "\033[0m";
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+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);
+ }
+ }
+ }
+}
+
diff --git a/NPLib/Physics/NPInelasticBreakup.h b/NPLib/Physics/NPInelasticBreakup.h
new file mode 100644
index 0000000000000000000000000000000000000000..178276cb9278ab5682ae6c7685e56bb78c63204b
--- /dev/null
+++ b/NPLib/Physics/NPInelasticBreakup.h
@@ -0,0 +1,153 @@
+#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"
+#include "NPReaction.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 fLabCrossSection;
+
+ private:
+ NPL::Reaction fReaction; // two body part of the kinematic
+ Beam fBeam; // Beam
+ Beam fParticle1; // Beam
+ Particle fParticle2; // Target
+ Particle fParticle3; // Light ejectile
+ Particle fParticle4; // Heavy ejectile
+ double fMeanEnergyLight;
+ double fSigmaEnergyLight;
+ double fMeanAngleLight;
+ double fSigmaAngleLight;
+ double fExcitationHeavy;
+ double fExcitationLight;
+
+ double fQValue; // Q-value in MeV
+ double fBeamEnergy; // Beam energy in MeV
+ double fThetaCM; // Center-of-mass angle in radian
+ 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 SetThetaCM(const double& angle) {
+ fThetaCM = angle;
+ initializePrecomputeVariable();
+ }
+
+ void SetExcitationLight(const double& exci) {
+ fExcitationLight = exci;
+ initializePrecomputeVariable();
+ }
+ void SetExcitationHeavy(const double& exci) {
+ fExcitationHeavy = 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 GetQValue() const { return fQValue; }
+ 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; }
+
+ 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();
+
+ public: // Kinematics
+ // Check that the reaction is alowed
+ bool IsLabCrossSection() { return fLabCrossSection; };
+
+ // Use fCrossSectionHist to shoot a Random ThetaCM and set fThetaCM to this value
+ double ShootRandomThetaCM();
+
+ // Compute ThetaLab and EnergyLab for product of reaction
+ void GenerateEvent(double& ThetaLight, double& KineticEnergyLight, double& ThetaLab3, double& KineticEnergyLab3,
+ double& ThetaLab4, double& KineticEnergyLab4);
+
+ void SetParticle3(double EnergyLab, double ThetaLab);
+
+ ClassDef(InelasticBreakup, 0)
+ };
+} // namespace NPL
+#endif