diff --git a/NPLib/Physics/CMakeLists.txt b/NPLib/Physics/CMakeLists.txt
index bb62685aeab5ac677d2903494c2636b34651e399..0c5817586a457fb4fe3f487d6244394836645d9a 100644
--- a/NPLib/Physics/CMakeLists.txt
+++ b/NPLib/Physics/CMakeLists.txt
@@ -2,6 +2,8 @@ add_custom_command(OUTPUT NPParticleDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts
add_custom_command(OUTPUT NPReactionDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPReaction.h NPReactionDict.cxx NPReaction.rootmap libNPPhysics.so NPReactionLinkDef.h DEPENDS NPReaction.h NPReactionLinkDef.h)
+add_custom_command(OUTPUT NPInelasticBreakupDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPInelasticBreakup.h NPInelasticBreakupDict.cxx NPInelasticBreakup.rootmap libNPPhysics.so NPInelasticBreakupLinkDef.h DEPENDS NPInelasticBreakup.h NPInelasticBreakupLinkDef.h)
+
add_custom_command(OUTPUT NPPhaseSpaceDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPPhaseSpace.h NPPhaseSpaceDict.cxx NPPhaseSpace.rootmap libNPPhysics.so NPPhaseSpaceLinkDef.h DEPENDS NPPhaseSpace.h NPPhaseSpaceLinkDef.h)
add_custom_command(OUTPUT NPQFSDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPQFS.h NPQFSDict.cxx NPQFS.rootmap libNPPhysics.so NPQFSLinkDef.h DEPENDS NPQFS.h NPQFSLinkDef.h)
@@ -26,7 +28,7 @@ add_custom_command(OUTPUT TFissionConditionsDict.cxx COMMAND ${CMAKE_BINARY_DIR}
add_custom_command(OUTPUT NPTimeStampDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPTimeStamp.h NPTimeStampDict.cxx NPTimeStamp.rootmap libNPPhysics.so NPTimeStampLinkDef.h DEPENDS NPTimeStamp.h NPTimeStampLinkDef.h)
-add_library(NPPhysics SHARED GEF.cxx NPFissionDecay.cxx NPDecay.cxx NPBeam.cxx NPEnergyLoss.cxx NPFunction.cxx NPParticle.cxx NPReaction.cxx NPTimeStamp.cxx NPPhaseSpace.cxx NPQFS.cxx NPParticleDict.cxx NPReactionDict.cxx NPTimeStampDict.cxx NPPhaseSpaceDict.cxx NPQFSDict.cxx NPEnergyLossDict.cxx )
+add_library(NPPhysics SHARED GEF.cxx NPInelasticBreakup.cxx NPInelasticBreakupDict.cxx NPFissionDecay.cxx NPDecay.cxx NPBeam.cxx NPEnergyLoss.cxx NPFunction.cxx NPParticle.cxx NPReaction.cxx NPTimeStamp.cxx NPPhaseSpace.cxx NPQFS.cxx NPParticleDict.cxx NPReactionDict.cxx NPTimeStampDict.cxx NPPhaseSpaceDict.cxx NPQFSDict.cxx NPEnergyLossDict.cxx )
target_link_libraries(NPPhysics ${ROOT_LIBRARIES} Physics NPCore)
add_library(NPInitialConditions SHARED TInitialConditions.cxx TInitialConditionsDict.cxx )
diff --git a/NPLib/Physics/NPInelasticBreakup.cxx b/NPLib/Physics/NPInelasticBreakup.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ee9690d8e0b2e3c8272bf40bc9f328fe438196fe
--- /dev/null
+++ b/NPLib/Physics/NPInelasticBreakup.cxx
@@ -0,0 +1,268 @@
+/*****************************************************************************
+ * 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();
+ fVerboseLevel = 1; // NPOptionManager::getInstance()->GetVerboseLevel();
+
+ 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& PhiLight, double& KineticEnergyLight,
+ double& ThetaLabTarget, double& KineticEnergyLabTarget, double& ThetaLabHeavy,
+ double& KineticEnergyLabHeavy) {
+ // 2-body relativistic kinematics: direct + inverse
+ // EnergieLabTarget,4 : lab energy in MeV of the 2 ejectiles
+ // ThetaLabTarget,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
+ ThetaLight = -1;
+ KineticEnergyLight = -1;
+ while (ThetaLight < 0 || KineticEnergyLight < 0) {
+ ThetaLight = gRandom->Gaus(fMeanAngleLight, fSigmaAngleLight);
+ KineticEnergyLight = gRandom->Gaus(fMeanEnergyLight, fSigmaEnergyLight);
+ }
+
+ PhiLight = gRandom->Uniform() * 2. * M_PI;
+ double totalELight = fParticle3.Mass() + KineticEnergyLight;
+ double momentumLight = sqrt(KineticEnergyLight * KineticEnergyLight + 2 * fParticle3.Mass() * KineticEnergyLight);
+ TLorentzVector lightLV(momentumLight * sin(ThetaLight) * cos(PhiLight),
+ momentumLight * sin(ThetaLight) * sin(PhiLight), momentumLight * cos(ThetaLight), totalELight);
+
+ // Substract light LV to Beam
+ double momentumBeam = sqrt(fBeamEnergy * fBeamEnergy + 2 * fBeam.Mass() * fBeamEnergy);
+ TLorentzVector beamLV(0, 0, momentumBeam, fBeamEnergy + fBeam.Mass());
+
+ auto newLV = beamLV - lightLV;
+
+ // Two body kine elastic scattering with heavy
+ fReaction.SetBeamEnergy(newLV.E() - newLV.Mag());
+ fReaction.ShootRandomThetaCM();
+ fReaction.KineRelativistic(ThetaLabTarget, KineticEnergyLabTarget, ThetaLabHeavy, KineticEnergyLabHeavy);
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+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("InelasticBreakup");
+ if (blocks.size() > 0 && NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "\033[1;35m//// Inelastic Breakup reaction found " << endl;
+
+ vector<string> token1 = {"Beam", "Target", "Light",
+ "Heavy", "MeanEnergyLight", "SigmaEnergyLight",
+ "MeanAngleLight", "SigmaAngleLight", "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.SetParticle1(fBeam);
+ fReaction.GetParticle1()->SetUp(blocks[i]->GetString("Heavy"));
+ fReaction.SetParticle2(GetParticle(blocks[i]->GetString("Target"), parser));
+ fReaction.SetParticle3(GetParticle(blocks[i]->GetString("Target"), parser));
+ fReaction.SetParticle4(GetParticle(blocks[i]->GetString("Heavy"), parser));
+ fReaction.SetBeamEnergy(fBeamEnergy);
+ fReaction.initializePrecomputeVariable();
+
+ fParticle3 = GetParticle(blocks[i]->GetString("Light"), 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);
+ fReaction.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);
+ 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..889ba3a727cb3d0312d75c3eacec0c84cfd43c33
--- /dev/null
+++ b/NPLib/Physics/NPInelasticBreakup.h
@@ -0,0 +1,141 @@
+#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; }
+ void SetThetaCM(const double& angle) { fThetaCM = angle; }
+
+ void SetExcitationLight(const double& exci) { fExcitationLight = exci; }
+ void SetExcitationHeavy(const double& exci) { fExcitationHeavy = exci; }
+ 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* GetParticleBeam() { return &fBeam; }
+ 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);
+
+ 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& PhiLight, double& KineticEnergyLight, double& ThetaLab3,
+ double& KineticEnergyLab3, double& ThetaLab4, double& KineticEnergyLab4);
+
+ void SetParticle3(double EnergyLab, double ThetaLab);
+
+ void SetParticle3(NPL::Particle p) { fParticle3 = p; };
+
+ ClassDef(InelasticBreakup, 0)
+ };
+} // namespace NPL
+#endif
diff --git a/NPLib/Physics/NPInelasticBreakupLinkDef.h b/NPLib/Physics/NPInelasticBreakupLinkDef.h
new file mode 100644
index 0000000000000000000000000000000000000000..00926376db86887caf81f16d993e67996c1d8c04
--- /dev/null
+++ b/NPLib/Physics/NPInelasticBreakupLinkDef.h
@@ -0,0 +1,3 @@
+#ifdef __CINT__
+#pragma link C++ class InelasticBreakup+;
+#endif
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 71388278fdde2cbc1408432c84c3071436e447df..62926e9ce52bb01f8e4c4419af3434ad23b1aaef 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -173,9 +173,9 @@ namespace NPL {
public:
// Modify the CS histo so cross section shoot is within ]HalfOpenAngleMin,HalfOpenAngleMax[
void SetCSAngle(double CSHalfOpenAngleMin, double CSHalfOpenAngleMax);
+ void initializePrecomputeVariable();
private: // intern precompute variable
- void initializePrecomputeVariable();
double m1;
double m2;
double m3;
@@ -243,7 +243,7 @@ namespace NPL {
// Return Excitation Energy
double ReconstructRelativistic(double EnergyLab, double ThetaLab, double PhiLab = 0);
-
+
// Return Lorentz Vector of Heavy Recoil after reaction
TLorentzVector LorentzAfterReaction(double EnergyLab, double ThetaLab, double PhiLab = 0);
@@ -261,6 +261,11 @@ namespace NPL {
void SetParticle3(double EnergyLab, double ThetaLab);
+ void SetParticle1(Beam p) { fParticle1 = p; };
+ void SetParticle2(Particle p) { fParticle2 = p; };
+ void SetParticle3(Particle p) { fParticle3 = p; };
+ void SetParticle4(Particle p) { fParticle4 = p; };
+
TGraph* GetKinematicLine3(double AngleStep_CM = 1);
TGraph* GetKinematicLine4(double AngleStep_CM = 1);
TGraph* GetBrhoLine3(double AngleStep_CM = 1);