diff --git a/NPSimulation/Process/ThreeBody.cc b/NPSimulation/Process/ThreeBody.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f8415e749b2f9520bd212ab1fbc6d06e8d617a87
--- /dev/null
+++ b/NPSimulation/Process/ThreeBody.cc
@@ -0,0 +1,532 @@
+/*****************************************************************************
+ * 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: Valerian Alcindor *
+ * contact address: valcindor@ikp.tu-darmstadt.de *
+ * *
+ * Creation Date : 2019 *
+ * Last update : August 2020 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Generation of sequential three body decay from resonant states *
+ * *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+#include "ThreeBody.hh"
+#include "G4Electron.hh"
+#include "G4Gamma.hh"
+#include "G4IonTable.hh"
+#include "G4ParticleTable.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4VPhysicalVolume.hh"
+#include "NPFunction.h"
+#include "NPInputParser.h"
+#include "NPOptionManager.h"
+#include "Particle.hh"
+#include "RootOutput.h"
+#include "TGenPhaseSpace.h"
+#include <iostream>
+#include <string>
+////////////////////////////////////////////////////////////////////////////////
+NPS::ThreeBody::ThreeBody(G4String modelName, G4Region* envelope)
+ : G4VFastSimulationModel(modelName, envelope) {
+ ReadConfiguration();
+ m_PreviousEnergy = 0;
+ m_PreviousLength = 0;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+NPS::ThreeBody::ThreeBody(G4String modelName)
+ : G4VFastSimulationModel(modelName) {}
+
+////////////////////////////////////////////////////////////////////////////////
+NPS::ThreeBody::~ThreeBody() {}
+
+////////////////////////////////////////////////////////////////////////////////
+void NPS::ThreeBody::AttachReactionConditions() {
+ // Reasssigned the branch address
+ if (RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
+ RootOutput::getInstance()->GetTree()->SetBranchAddress(
+ "ReactionConditions", &m_ReactionConditions);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void NPS::ThreeBody::ReadConfiguration() {
+ NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
+
+ vector<NPL::InputBlock*> blocks = input.GetAllBlocksWithToken("ThreeBody");
+ if (blocks.size() > 0) {
+ m_active = true;
+ m_ReactionConditions = new TReactionConditions();
+ AttachReactionConditions();
+ if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
+ RootOutput::getInstance()->GetTree()->Branch(
+ "ReactionConditions", "TReactionConditions", &m_ReactionConditions);
+
+ for (unsigned int i = 0; i < blocks.size(); i++) {
+ if (blocks[i]->HasToken("Type"))
+ m_type = blocks[i]->GetString("Type");
+
+ if (blocks[i]->HasToken("Beam"))
+ m_Beam = blocks[i]->GetString("Beam");
+
+ if (blocks[i]->HasToken("Target"))
+ m_TargetNucleus = blocks[i]->GetString("Target");
+
+ if (blocks[i]->HasToken("CompoundNucleusEnergy"))
+ m_CompoundNucleusEnergy
+ = blocks[i]->GetDouble("CompoundNucleusEnergy", "MeV");
+
+ if (blocks[i]->HasToken("CompoundNucleusStateWidth"))
+ m_CompoundNucleusStateWidth
+ = blocks[i]->GetDouble("CompoundNucleusStateWidth", "MeV");
+
+ if (blocks[i]->HasToken("IntermediaryNucleus"))
+ m_IntermediaryNucleus = blocks[i]->GetString("IntermediaryNucleus");
+
+ if (blocks[i]->HasToken("IntermediaryState"))
+ m_IntermediaryState = blocks[i]->GetDouble("IntermediaryState", "MeV");
+
+ if (blocks[i]->HasToken("IntermediaryStateWidth"))
+ m_IntermediaryStateWidth
+ = blocks[i]->GetDouble("IntermediaryStateWidth", "MeV");
+
+ if (blocks[i]->HasToken("Light1"))
+ m_Light1 = blocks[i]->GetString("Light1");
+
+ if (blocks[i]->HasToken("Heavy"))
+ m_Heavy = blocks[i]->GetString("Heavy");
+
+ if (blocks[i]->HasToken("Light2"))
+ m_Light2 = blocks[i]->GetString("Light2");
+
+ if (blocks[i]->HasToken("ExcitationEnergyHeavy"))
+ m_ExcitationEnergyHeavy
+ = blocks[i]->GetDouble("ExcitationEnergyHeavy", "MeV");
+
+ if (blocks[i]->HasToken("ShootLight1"))
+ m_ShootLight1 = blocks[i]->GetInt("ShootLight1");
+
+ if (blocks[i]->HasToken("ShootLight2"))
+ m_ShootLight2 = blocks[i]->GetInt("ShootLight2");
+
+ if (blocks[i]->HasToken("ShootHeavy"))
+ m_ShootHeavy = blocks[i]->GetInt("ShootHeavy");
+
+ m_UserEnergyCS = false;
+ if (blocks[i]->HasToken("UserEnergyCS")) {
+ m_UserEnergyCS = blocks[i]->GetInt("UserEnergyCS");
+
+ cout << m_UserEnergyCS << endl;
+
+ if (blocks[i]->HasToken("EnergyCSPath"))
+ m_EnergyCSPath = blocks[i]->GetString("EnergyCSPath");
+
+ if (blocks[i]->HasToken("EnergyCSName"))
+ m_EnergyCSName = blocks[i]->GetString("EnergyCSName");
+
+ if (m_UserEnergyCS) {
+ m_FileEnergyCS = new TFile(m_EnergyCSPath.c_str(), "open");
+ m_EnergyCS
+ = (TF1*)m_FileEnergyCS->FindObjectAny(m_EnergyCSName.c_str());
+ m_EnergyCSMax = m_EnergyCS->GetMaximum();
+ m_EnergyCS->SetNpx(1e6);
+ m_FileEnergyCS->Close();
+ }
+ } else {
+ cout << "ERROR: check your input file formatting \033[0m" << endl;
+ exit(1);
+ }
+ }
+ m_BeamName = NPL::ChangeNameToG4Standard(m_Beam);
+ m_N1 = NPL::Particle(m_Beam);
+ m_N2 = NPL::Particle(m_TargetNucleus);
+
+ m_N3 = NPL::Particle(m_Light1);
+ m_N4 = NPL::Particle(m_IntermediaryNucleus);
+
+ m_N5 = NPL::Particle(m_Light2);
+ m_N6 = NPL::Particle(m_Heavy);
+ m_end = false;
+ m_previousE = -1000;
+ fIntermediaryStateEnergy = GetIntermediaryStateEnergy(
+ m_IntermediaryState, m_IntermediaryStateWidth);
+ } else
+ m_active = false;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+G4bool NPS::ThreeBody::IsApplicable(const G4ParticleDefinition& particleType) {
+ NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
+
+ vector<NPL::InputBlock*> blocks = input.GetAllBlocksWithToken("ThreeBody");
+ if (!m_active)
+ return false;
+ if (!m_active)
+ return false;
+ else {
+ static std::string particleName;
+ particleName = particleType.GetParticleName();
+ if (particleName.find(m_BeamName) != std::string::npos
+ && blocks.size() > 0) {
+ return true;
+ }
+ }
+ return false;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+G4bool NPS::ThreeBody::ModelTrigger(const G4FastTrack& fastTrack) {
+ static bool shoot = false;
+ static double rand = 0;
+ const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
+ G4ThreeVector V = PrimaryTrack->GetMomentum().unit();
+ G4ThreeVector P = fastTrack.GetPrimaryTrackLocalPosition();
+ G4VSolid* solid = fastTrack.GetPrimaryTrack()
+ ->GetVolume()
+ ->GetLogicalVolume()
+ ->GetSolid();
+ double in = solid->DistanceToOut(P, V);
+ double out = solid->DistanceToOut(P, -V);
+ double ratio = in / (out + in);
+
+ if (out == 0) { // first step of current event
+ rand = G4RandFlat::shoot();
+ m_PreviousLength = m_StepSize;
+ m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
+ // Clear Previous Event
+ m_ReactionConditions->Clear();
+ shoot = true;
+ } else if (in == 0) { // last step
+ return true;
+ }
+
+ // If the condition is met, the event is generated
+ if (ratio < rand) {
+ // Reset the static for next event
+ if (shoot) {
+ // && m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy())) {
+ shoot = false;
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ // Record the situation of the current step
+ // so it can be used in the next one
+ if (!PrimaryTrack->GetStep()->IsLastStepInVolume()) {
+ m_PreviousLength = PrimaryTrack->GetStep()->GetStepLength();
+ m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
+ }
+
+ return false;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void NPS::ThreeBody::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) {
+ // Get the track info
+ const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
+ G4ThreeVector pdirection = PrimaryTrack->GetMomentum().unit();
+ G4ThreeVector localdir = fastTrack.GetPrimaryTrackLocalDirection();
+
+ G4ThreeVector worldPosition = PrimaryTrack->GetPosition();
+ G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition();
+
+ double energy = PrimaryTrack->GetKineticEnergy();
+ double time = PrimaryTrack->GetGlobalTime();
+
+ /////////////////////////
+ /// Beam Parameters /////
+ /////////////////////////
+ m_ReactionConditions->SetBeamParticleName(
+ PrimaryTrack->GetParticleDefinition()->GetParticleName());
+
+ G4ThreeVector U(1, 0, 0);
+ G4ThreeVector V(0, 1, 0);
+ G4ThreeVector ZZ(0, 0, 1);
+ m_ReactionConditions->SetBeamEmittanceTheta(
+ PrimaryTrack->GetMomentumDirection().theta() / deg);
+ m_ReactionConditions->SetBeamEmittancePhi(
+ PrimaryTrack->GetMomentumDirection().phi() / deg);
+ m_ReactionConditions->SetBeamEmittanceThetaX(
+ PrimaryTrack->GetMomentumDirection().angle(U) / deg);
+ m_ReactionConditions->SetBeamEmittancePhiY(
+ PrimaryTrack->GetMomentumDirection().angle(V) / deg);
+
+ // Randomize within the step
+ // Assume energy loss is linear within the step
+ // Assume no scattering
+ double rand = G4RandFlat::shoot();
+ double length = rand * (m_PreviousLength);
+ energy -= (1 - rand) * (m_PreviousEnergy - energy);
+ G4ThreeVector ldir = pdirection;
+ ldir *= length;
+ localPosition = localPosition - ldir;
+
+ m_ReactionConditions->SetBeamReactionEnergy(energy);
+ m_ReactionConditions->SetVertexPositionX(localPosition.x());
+ m_ReactionConditions->SetVertexPositionY(localPosition.y());
+ m_ReactionConditions->SetVertexPositionZ(localPosition.z());
+
+ // Set the end of the step conditions
+ fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0, pdirection);
+ fastStep.SetPrimaryTrackFinalPosition(worldPosition);
+ fastStep.SetTotalEnergyDeposited(0);
+ fastStep.SetPrimaryTrackFinalTime(time);
+ fastStep.KillPrimaryTrack();
+ fastStep.SetPrimaryTrackPathLength(0.0);
+
+ double Ecm = sqrt(m_N1.Mass() * m_N1.Mass() + m_N2.Mass() * m_N2.Mass()
+ + 2 * m_N2.Mass() * (energy + m_N1.Mass()))
+ - m_N1.Mass() - m_N2.Mass();
+ if (m_UserEnergyCS == 0) {
+ Ecm = m_CompoundNucleusEnergy;
+ }
+
+ bool Allowed = false;
+ if (Ecm < m_IntermediaryState)
+ Ecm = 0;
+ if (m_UserEnergyCS) {
+ double cs_value = m_EnergyCS->Eval(Ecm) / m_EnergyCSMax;
+ double rand_cs = G4RandFlat::shoot();
+ if (rand_cs < cs_value) {
+ Allowed = true;
+ }
+ } else if (Ecm > m_IntermediaryState) {
+ Allowed = true;
+ }
+
+ if (Allowed) {
+ /////////////////////////////////////////////////
+ //////Define the kind of particle to shoot////////
+ //////////////////////////////////////////////////
+
+ // Build the decaying particle four momenta vector:
+ double BeamEnergy = energy;
+
+ G4ThreeVector BeamMomentumG4V = pdirection;
+
+ // Before reaction
+ double m1 = m_N1.Mass();
+ double m2 = m_N2.Mass();
+
+ // After first emission
+ double m3 = m_N3.Mass();
+
+ double IntermediaryStateEnergy = 0.;
+ double SpSn = m_N4.GetBindingEnergy() - m_N6.GetBindingEnergy();
+ if (m_IntermediaryStateWidth == 0) {
+ IntermediaryStateEnergy = m_IntermediaryState;
+ } else {
+ while (true) {
+ IntermediaryStateEnergy = fIntermediaryStateEnergy->GetRandom();
+ if (IntermediaryStateEnergy > SpSn && IntermediaryStateEnergy < Ecm)
+ break;
+ }
+ }
+
+ m_N4.SetExcitationEnergy(IntermediaryStateEnergy);
+ double m4 = m_N4.Mass();
+
+ // After second emission
+ double m5 = m_N5.Mass();
+ double m6 = m_N6.Mass();
+
+ std::vector<NPL::Particle> DaughterNuclei;
+ DaughterNuclei.push_back(m_N3);
+ DaughterNuclei.push_back(m_N5);
+ DaughterNuclei.push_back(m_N6);
+
+ double T1 = BeamEnergy;
+ if (m_UserEnergyCS == 0) {
+ T1 = (pow((Ecm + m_N1.Mass() + m_N2.Mass()), 2.)
+ - (m_N1.Mass() * m_N1.Mass() + m_N2.Mass() * m_N2.Mass()))
+ / (2. * m_N2.Mass())
+ - m_N1.Mass();
+ }
+
+ TVector3 p1(BeamMomentumG4V.x(), BeamMomentumG4V.y(), BeamMomentumG4V.z());
+
+ double gamma1 = T1 / m1 + 1.;
+ double Beta1 = sqrt(1. - 1. / (gamma1 * gamma1));
+ p1.SetMag(gamma1 * m1 * Beta1);
+ TLorentzVector LV_1(p1.x(), p1.y(), p1.z(), T1 + m1);
+
+ TVector3 p2(0, 0, 0);
+ TLorentzVector LV_2(p2.x(), p2.y(), p2.z(), m2);
+
+ TLorentzVector LV_CN = LV_1 + LV_2;
+
+ /////////////////////////
+ // Sequential emission //
+ /////////////////////////
+
+ if (m_type == "sequential") {
+ // First particle emission
+ double T3, T4 = 0;
+ unsigned int nDecay = 2;
+ double* masses34 = new double[nDecay];
+ masses34[0] = m3 / 1000.;
+ masses34[1] = m4 / 1000.;
+ double* T34 = new double[nDecay];
+ double* theta34 = new double[nDecay];
+ double* phi34 = new double[nDecay];
+
+ GetDecay(LV_CN.Vect(), LV_CN.E(), 2, masses34, T34, theta34, phi34);
+ TVector3 p3(0, 0, 1);
+ p3.SetTheta(theta34[0]);
+ p3.SetPhi(phi34[0]);
+ T3 = T34[0];
+ double gamma3 = T3 / m3 + 1;
+ double Beta3 = sqrt(1 - 1 / (gamma3 * gamma3));
+ p3.SetMag(gamma3 * m3 * Beta3);
+ G4ThreeVector p3G4 = G4ThreeVector(p3.x(), p3.y(), p3.z());
+
+ TVector3 p4(0, 0, 1);
+ p4.SetTheta(theta34[1]);
+ p4.SetPhi(phi34[1]);
+ T4 = T34[1];
+ double gamma4 = T4 / m4 + 1;
+ double Beta4 = sqrt(1 - 1 / (gamma4 * gamma4));
+ p4.SetMag(gamma4 * m4 * Beta4);
+ TLorentzVector LV_4(p4.x(), p4.y(), p4.z(), T4 + m4);
+
+ // Second particle emission
+ double T5, T6 = 0;
+ nDecay = 2;
+ double* masses56 = new double[nDecay];
+ masses56[0] = m5 / 1000.;
+ masses56[1] = m6 / 1000.;
+ double* T56 = new double[nDecay];
+ double* theta56 = new double[nDecay];
+ double* phi56 = new double[nDecay];
+ GetDecay(LV_4.Vect(), LV_4.E(), 2, masses56, T56, theta56, phi56);
+ TVector3 p5(0, 0, 1);
+ p5.SetTheta(theta56[0]);
+ p5.SetPhi(phi56[0]);
+ T5 = T56[0];
+ double gamma5 = T5 / m5 + 1;
+ double Beta5 = sqrt(1 - 1 / (gamma5 * gamma5));
+ p5.SetMag(gamma5 * m5 * Beta5);
+ G4ThreeVector p5G4 = G4ThreeVector(p5.x(), p5.y(), p5.z());
+
+ TVector3 p6(0, 0, 1);
+ p6.SetTheta(theta56[1]);
+ p6.SetPhi(phi56[1]);
+ T6 = T56[1];
+ double gamma6 = T6 / m6 + 1;
+ double Beta6 = sqrt(1 - 1 / (gamma6 * gamma6));
+ p6.SetMag(gamma6 * m6 * Beta6);
+ G4ThreeVector p6G4 = G4ThreeVector(p6.x(), p6.y(), p6.z());
+
+ G4ParticleDefinition* Light1
+ = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(
+ m_N3.GetZ(), m_N3.GetA(), 0);
+
+ G4ParticleDefinition* Light2
+ = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(
+ m_N5.GetZ(), m_N5.GetA(), 0);
+
+ G4ParticleDefinition* Heavy
+ = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(
+ m_N6.GetZ(), m_N6.GetA(), m_ExcitationEnergyHeavy);
+
+ TLorentzVector LVproton1(p3.x(), p3.y(), p3.z(), T3 + m3);
+ TLorentzVector LVproton2(p5.x(), p5.y(), p5.z(), T5 + m5);
+ TLorentzVector LVHeavy(p6.x(), p6.y(), p6.z(), T6 + m6);
+
+ TLorentzVector All = LVproton1 + LVproton2 + LVHeavy;
+
+ // Emitt secondary
+ if (m_ShootLight1) {
+ G4DynamicParticle particleLight1(Light1, p3G4.unit(), T3);
+ fastStep.CreateSecondaryTrack(particleLight1, localPosition, time);
+ }
+
+ if (m_ShootLight2) {
+ G4DynamicParticle particleLight2(Light2, p5G4.unit(), T5);
+ fastStep.CreateSecondaryTrack(particleLight2, localPosition, time);
+ }
+ if (m_ShootHeavy) {
+ G4DynamicParticle particleHeavy(Heavy, p6G4.unit(), T6);
+ fastStep.CreateSecondaryTrack(particleHeavy, localPosition, time);
+ }
+ ///////////////////////////////////////
+ ///// Emitted particle Parameters /////
+ ///////////////////////////////////////
+ m_ReactionConditions->SetParticleName(Light1->GetParticleName());
+ m_ReactionConditions->SetKineticEnergy(T3);
+ m_ReactionConditions->SetTheta(theta34[0] / deg);
+ if ((phi34[0] + pi) / deg > 360)
+ m_ReactionConditions->SetPhi((phi34[0] - pi) / deg);
+ else
+ m_ReactionConditions->SetPhi((phi34[0] + pi) / deg);
+
+ m_ReactionConditions->SetParticleName(Light2->GetParticleName());
+ m_ReactionConditions->SetKineticEnergy(T5);
+ m_ReactionConditions->SetTheta(theta56[0] / deg);
+ if ((phi56[0] + pi) / deg > 360)
+ m_ReactionConditions->SetPhi((phi56[0] - pi) / deg);
+ else
+ m_ReactionConditions->SetPhi((phi56[0] + pi) / deg);
+
+ m_ReactionConditions->SetParticleName(Heavy->GetParticleName());
+ m_ReactionConditions->SetKineticEnergy(T6);
+ m_ReactionConditions->SetTheta(theta56[1] / deg);
+ if ((phi56[1] + pi) / deg > 360)
+ m_ReactionConditions->SetPhi((phi56[1] - pi) / deg);
+ else
+ m_ReactionConditions->SetPhi((phi56[1] + pi) / deg);
+ }
+ }
+
+ // Reinit for next event
+ m_PreviousEnergy = 0;
+ m_PreviousLength = 0;
+}
+
+void NPS::ThreeBody::GetDecay(TVector3 p1, double Etot1, unsigned int nDecay,
+ double* masses, double* T, double* theta,
+ double* phi) {
+ TGenPhaseSpace PhaseSpace;
+ TLorentzVector NucleiToDecay(p1.x() / 1000., p1.y() / 1000., p1.z() / 1000.,
+ Etot1 / 1000.);
+
+ if (!PhaseSpace.SetDecay(NucleiToDecay, nDecay, masses)) {
+ cout << "FORBIDDEN PHASE SPACE CHECK ENERGIES" << endl;
+ for (unsigned int i = 0; i < nDecay; i++) {
+ T[i] = (PhaseSpace.GetDecay(i)->E() - masses[i]) * 1000.;
+ theta[i] = PhaseSpace.GetDecay(i)->Theta();
+ phi[i] = PhaseSpace.GetDecay(i)->Phi();
+ }
+ exit(1);
+ } else {
+ PhaseSpace.Generate();
+ for (unsigned int i = 0; i < nDecay; i++) {
+ T[i] = (PhaseSpace.GetDecay(i)->E() - masses[i]) * 1000.;
+ theta[i] = PhaseSpace.GetDecay(i)->Theta();
+ phi[i] = PhaseSpace.GetDecay(i)->Phi();
+ }
+ }
+}
+
+TF1* NPS::ThreeBody::GetIntermediaryStateEnergy(double Er, double Width) {
+
+ TF1* f
+ = new TF1("f", "1. / (pow((x - [0]), 2.) + pow(([1] / 2.), 2.))", 0, 100);
+ f->SetParameter(0, Er);
+ f->SetParameter(1, Width);
+ f->SetNpx(1e6);
+ return f;
+}
diff --git a/NPSimulation/Process/ThreeBody.hh b/NPSimulation/Process/ThreeBody.hh
new file mode 100644
index 0000000000000000000000000000000000000000..91a3ec1be96d7c7926e36edbb7ebe5e88dd85c62
--- /dev/null
+++ b/NPSimulation/Process/ThreeBody.hh
@@ -0,0 +1,107 @@
+/*****************************************************************************
+ * 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 address: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : Octobre 2017 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Use to kill the beam track and replace it with the reaction product *
+ * *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+#ifndef ThreeBody_h
+#define ThreeBody_h
+
+#include "G4VFastSimulationModel.hh"
+#include "NPReaction.h"
+#include "PhysicsList.hh"
+#include "TReactionConditions.h"
+
+#include "TF1.h"
+#include "TFile.h"
+
+class G4VPhysicalVolume;
+namespace NPS {
+class ThreeBody : public G4VFastSimulationModel {
+public:
+ ThreeBody(G4String, G4Region*);
+ ThreeBody(G4String);
+ ~ThreeBody();
+
+public:
+ void ReadConfiguration();
+ G4bool IsApplicable(const G4ParticleDefinition&);
+ G4bool ModelTrigger(const G4FastTrack&);
+ void DoIt(const G4FastTrack&, G4FastStep&);
+
+private:
+ NPL::Reaction m_Reaction;
+ string m_BeamName;
+ double m_PreviousEnergy;
+ double m_PreviousLength;
+ bool m_active; // is the process active
+ double m_StepSize;
+
+private:
+ TReactionConditions* m_ReactionConditions;
+
+public:
+ void AttachReactionConditions();
+ void SetStepSize(double step) { m_StepSize = step; };
+
+public:
+ bool m_end;
+ double m_previousE;
+ std::string m_type;
+ std::string m_Beam;
+ std::string m_TargetNucleus;
+ double m_CompoundNucleusEnergy;
+ double m_CompoundNucleusStateWidth;
+ std::string m_IntermediaryNucleus;
+ double m_IntermediaryState;
+ double m_IntermediaryStateWidth;
+ std::string m_Light1;
+ std::string m_Heavy;
+ std::string m_Light2;
+ double m_ExcitationEnergyHeavy;
+ int m_ShootLight1;
+ int m_ShootLight2;
+ int m_ShootHeavy;
+ NPL::Particle m_N1;
+ NPL::Particle m_N2;
+
+ NPL::Particle m_N3;
+ NPL::Particle m_N4;
+
+ NPL::Particle m_N5;
+ NPL::Particle m_N6;
+
+ bool m_UserEnergyCS;
+ std::string m_EnergyCSPath;
+ std::string m_EnergyCSName;
+ TFile* m_FileEnergyCS;
+ TF1* m_EnergyCS;
+ double m_EnergyCSMax;
+ TF1* fIntermediaryStateEnergy;
+
+public: // Managing the decay
+ // Set everything for the decay
+ void GetDecay(TVector3 p, double Etot, unsigned int nDecay, double* masses,
+ double* T, double* theta, double* phi);
+ TF1* GetIntermediaryStateEnergy(double Er, double Width);
+};
+} // namespace NPS
+
+#endif
diff --git a/Projects/MultiParticleRES/.DS_Store b/Projects/MultiParticleRES/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..778fb5116f78b23a659c81f76fafdfb8fa17de68
Binary files /dev/null and b/Projects/MultiParticleRES/.DS_Store differ
diff --git a/Projects/MultiParticleRES/Inputs/.DS_Store b/Projects/MultiParticleRES/Inputs/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6
Binary files /dev/null and b/Projects/MultiParticleRES/Inputs/.DS_Store differ
diff --git a/Projects/MultiParticleRES/Inputs/BW.root b/Projects/MultiParticleRES/Inputs/BW.root
new file mode 100644
index 0000000000000000000000000000000000000000..2838ace698508d1f053ff30ce8ec6fab773189ca
Binary files /dev/null and b/Projects/MultiParticleRES/Inputs/BW.root differ
diff --git a/Projects/MultiParticleRES/Inputs/ExampleSequential15F.detector b/Projects/MultiParticleRES/Inputs/ExampleSequential15F.detector
new file mode 100755
index 0000000000000000000000000000000000000000..bfc19a7ede9a20ab79fb30003940ba1806e8bc8e
--- /dev/null
+++ b/Projects/MultiParticleRES/Inputs/ExampleSequential15F.detector
@@ -0,0 +1,75 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Target
+ THICKNESS= 106.5 micrometer
+ RADIUS= 30 mm
+ MATERIAL= CH2
+ ANGLE= 0 deg
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 cm
+
+% %%%%%%%%%%%%%%%%%%%%%%%%%
+% % Experiment thickness = 34um
+% beam_dump
+% R= 830 mm
+% THETA = 0
+% PHI = 0
+% Thickness= 45 micrometer
+%
+%%%%%%% Telescope 1 DSSD5 %%%%%%%
+M2Telescope
+ X128_Y128= -12.05 104.5 147.54
+ X128_Y1= -102.43 104.7 111.59
+ X1_Y1= -114.08 12.53 140.36
+ X1_Y128= -23.71 12.32 176.3
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+
+%%%%%%% Telescope 2 DSSD 4 %%%%%%%
+M2Telescope
+ X128_Y128= -114.49 -10.46 140.4
+ X128_Y1= -103.09 -102.67 111.64
+ X1_Y1= -12.56 -102.58 147.47
+ X1_Y128= -23.97 -10.38 176.22
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+
+%%%%%%% Telescope 3 DSSD 7 %%%%%%%
+M2Telescope
+ X128_Y128= 13.11 -103.21 147.01
+ X128_Y1= 103.38 -103.24 110.82
+ X1_Y1= 114.96 -10.99 139.34
+ X1_Y128= 24.66 -10.9 175.56
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+
+%%%%%%% Telescope 4 DSSD 11 %%%%%%%
+M2Telescope
+ X128_Y128= 115.09 11.03 140.54
+ X128_Y1= 103.67 103.32 111.26
+ X1_Y1= 13.52 103.53 147.05
+ X1_Y128= 24.94 11.24 176.34
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+
+% %%%%%%% Telescope 5 DSSD 2 %%%%%%%
+% M2Telescope
+% X128_Y128= -48.51 47.21 860.3
+% X128_Y1= -48.84 -50.07 860.45
+% X1_Y1= 48.45 -50.36 860.26
+% X1_Y128= 48.81 46.91 860.11
+% SI= 1
+% SILI= 0
+% CSI= 1
+% VIS= all
+
+
+
diff --git a/Projects/MultiParticleRES/Inputs/ExampleSequential15F.reaction b/Projects/MultiParticleRES/Inputs/ExampleSequential15F.reaction
new file mode 100644
index 0000000000000000000000000000000000000000..de9c8857e273013f85d3f491dae93affd010a83e
--- /dev/null
+++ b/Projects/MultiParticleRES/Inputs/ExampleSequential15F.reaction
@@ -0,0 +1,38 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 14O
+ Energy= 103.885 MeV
+ %Energy= 93.8 MeV
+ SigmaEnergy= 0.005 MeV
+ SigmaThetaX= 0. deg
+ SigmaPhiY= 0. deg
+ %SigmaX= 7 mm
+ %SigmaY= 7 mm
+ SigmaX= 8 mm
+ SigmaY= 8 mm
+ MeanThetaX= 0. deg
+ MeanPhiY= 0. deg
+ MeanX= -4. mm
+ MeanY= -1. mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ThreeBody
+ Type= sequential
+ Beam= 14O
+ Target= 1H
+ CompoundNucleusEnergy= 6.3 MeV
+ CompoundNucleusStateWidth= 0.02 MeV
+ IntermediaryNucleus= 14O
+ IntermediaryState= 0 MeV
+ IntermediaryStateWidth= 0.01 MeV
+ Light1= 1H
+ Heavy= 13N
+ Light2= 1H
+ ExcitationEnergyHeavy= 0. MeV
+ ShootLight1= 1
+ ShootLight2= 1
+ ShootHeavy= 1
+ UserEnergyCS= 1
+ EnergyCSPath= $NPTOOL/Projects/MultiParticleRES/Inputs/BW.root
+ EnergyCSName= BW
diff --git a/Projects/MultiParticleRES/Inputs/simu.mac b/Projects/MultiParticleRES/Inputs/simu.mac
new file mode 100644
index 0000000000000000000000000000000000000000..add2e059e4e5a5806132ca1713ec1c5a75dd1a49
--- /dev/null
+++ b/Projects/MultiParticleRES/Inputs/simu.mac
@@ -0,0 +1,2 @@
+/run/beamOn 100000
+
diff --git a/Projects/MultiParticleRES/simu.sh b/Projects/MultiParticleRES/simu.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b85e7e901680feaf0991af32d5e4456308ab9ea1
--- /dev/null
+++ b/Projects/MultiParticleRES/simu.sh
@@ -0,0 +1,2 @@
+npsimulation -D Inputs/ExampleSequential15F.detector -E Inputs/ExampleSequential15F.reaction -B Inputs/simu.mac -O testMultiParticuleRES
+
diff --git a/Projects/MultiParticleRES/utilities/BW.root b/Projects/MultiParticleRES/utilities/BW.root
new file mode 100644
index 0000000000000000000000000000000000000000..2838ace698508d1f053ff30ce8ec6fab773189ca
Binary files /dev/null and b/Projects/MultiParticleRES/utilities/BW.root differ
diff --git a/Projects/MultiParticleRES/utilities/ThreeBodyCS.cc b/Projects/MultiParticleRES/utilities/ThreeBodyCS.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5674aae83cd7c69573a1d313117b4f4acec1ed56
--- /dev/null
+++ b/Projects/MultiParticleRES/utilities/ThreeBodyCS.cc
@@ -0,0 +1,267 @@
+/*****************************************************************************
+ * 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: V. Alcindor mail: valcindor@ikp.tu-darmstadt.de *
+ * *
+ * Creation Date : July 2022 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This files is for the generation of 15F Breit Wigner shape in the case *
+ * of the sequential 2p emission 14O+p -> 15F -> 14O(1-) + p -> 13N + p *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ * *
+ *****************************************************************************/
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_sf_coulomb.h>
+
+// dirty way of avoiding a header or input file...
+namespace input_conditions {
+// Incoming beam, here 14O
+NPL::Particle Beam(8, 14);
+// Target particle, here protons
+NPL::Particle Target(1, 1);
+// Intermediate state in the sequential emission (14O(1-), Ex=5.173 MeV)
+NPL::Particle Intermediate(8, 14);
+double IntermediateState = 5.173;
+// Final particle in the exit channel, here 13N
+NPL::Particle Final(7, 13);
+} // namespace input_conditions
+
+using namespace input_conditions;
+
+double CoulombFunction(std::string type, double E, double R, double mu,
+ Int_t Z1, Int_t Z2, double L, Int_t k = 0) {
+ gsl_sf_result F_res;
+ gsl_sf_result Fp_res;
+ gsl_sf_result G_res;
+ gsl_sf_result Gp_res;
+ double expF;
+ double expG;
+ double eta, rho;
+ eta = 0.15748 * Z1 * Z2 * sqrt(mu / E);
+ rho = 0.218735 * sqrt(mu * E) * R;
+ gsl_sf_coulomb_wave_FG_e(eta, rho, L, k, &F_res, &Fp_res, &G_res, &Gp_res,
+ &expF, &expG);
+ if (type == "F")
+ return F_res.val * exp(expF);
+ else if (type == "G")
+ return G_res.val * exp(expG);
+ else
+ return 1;
+}
+
+double rho(double* E, double* par) {
+ double R = par[0];
+ double mu_uma = par[1];
+
+ return R * .218735 * sqrt(mu_uma * E[0]);
+} // computation of the reduced variable for the penetrability
+
+double P_l(double* E, double* par) {
+
+ gsl_sf_result F_res;
+ gsl_sf_result Fp_res;
+ gsl_sf_result G_res;
+ gsl_sf_result Gp_res;
+
+ double R = par[0];
+ double mu_uma = par[1];
+ double Z1 = par[2];
+ double Z2 = par[3];
+ double L = par[4];
+ double par_P[2] = {par[0], par[1]};
+ double coulF = CoulombFunction("F", E[0], R, mu_uma, Z1, Z2, L, 0);
+ double coulG = CoulombFunction("G", E[0], R, mu_uma, Z1, Z2, L, 0);
+ double Penetrability = rho(E, par_P) / (pow(coulF, 2.) + pow(coulG, 2.));
+
+ return Penetrability;
+} // computation of the penetrability cf Abrahamovitz & Stegun Ch14 et C.
+
+double f_theoretical_proton_width(double* x, double* par) {
+ double Er = par[0];
+ double l = par[1];
+ double z = par[2];
+ double a = par[3];
+ double Z = par[4];
+ double A = par[5];
+ double E = x[0];
+
+ double r0 = 1.2; // fm
+ double R = r0 * (pow(A, 1. / 3.) + pow(a, 1. / 3.)); // fm
+
+ double mBeam_uma = Beam.GetA();
+ +Beam.GetMassExcess() / 1000.;
+ double mTarget_uma = Target.GetA();
+ +Target.GetMassExcess() / 1000.;
+ double mCN_uma = (mBeam_uma * mTarget_uma) / (mBeam_uma + mTarget_uma);
+
+ double P_par[5] = {R, mCN_uma, z, Z, l};
+ double P = P_l(&E, P_par);
+
+ double h_bar_c = 197.3; // MeV.fm
+ double max_proton_width
+ = 3. * pow(h_bar_c, 2.) / (((A * a) / (A + a) * 931.5) * pow(R, 2.));
+
+ double proton_width = max_proton_width * P;
+
+ if (E > 3.339 && E < 3.99) {
+ double x1 = 3.339;
+ double x2 = 3.99;
+ P = P_l(&x1, P_par);
+ double y1 = max_proton_width * P;
+ P = P_l(&x2, P_par);
+ double y2 = max_proton_width * P;
+
+ return (y1 - y2) / (x1 - x2) * E + (y2 - x2 * (y1 - y2) / (x1 - x2));
+ } // trying to correct GSLCoulomb problem in this region
+
+ return proton_width;
+} // Computation of proton width using COULOMB function from GSL
+
+double f_proton_width(double* x, double* par) {
+ double Er = par[0];
+ double res_width = par[1];
+ double l = par[2];
+ double z = par[3];
+ double a = par[4];
+ double Z = par[5];
+ double A = par[6];
+ double E = x[0];
+
+ double r0 = 1.2; // fm
+ double R = r0 * (pow(A, 1. / 3.) + pow(a, 1. / 3.)); // fm
+
+ double mBeam_uma = Beam.GetA() + Beam.GetMassExcess() / 1000.;
+ double mTarget_uma = Target.GetA() + Target.GetMassExcess() / 1000.;
+ double mCN_uma = (mBeam_uma * mTarget_uma) / (mBeam_uma + mTarget_uma);
+
+ double P_par[5] = {R, mCN_uma, z, Z, l};
+ double P = P_l(&E, P_par);
+ double Pr = P_l(&Er, P_par);
+
+ double proton_width = (res_width * P / Pr);
+
+ if (E > 3.339 && E < 3.99) {
+ double x1 = 3.339;
+ double x2 = 3.99;
+ P = P_l(&x1, P_par);
+ double y1 = (res_width * P / Pr);
+ P = P_l(&x2, P_par);
+ double y2 = (res_width * P / Pr);
+
+ return (y1 - y2) / (x1 - x2) * E + (y2 - x2 * (y1 - y2) / (x1 - x2));
+ } // trying to correct GSLCoulomb problem in this region for 15F
+
+ return proton_width;
+} // Computation of proton width using COULOMB function from GSL
+
+double f_double_sequential_proton_cross_section(double* x, double* par) {
+ double E = x[0]; // MeV
+ if (E < 0.) {
+ std::cerr << "ENERGY IS NEGATIVE" << std::endl;
+ return 1.;
+ } else {
+
+ double Er = par[0]; // MeV
+ double res_width = par[1]; // MeV
+
+ double J1 = par[2];
+ double J2 = par[3];
+ double J = J1 + J2;
+ double w = (2. * J + 1.) / ((2. * J1 + 1.) * (2. * J2 + 1.));
+
+ double li = par[4];
+ double lf = par[5];
+
+ double p = std::sqrt(2. * 14. / 15. * E); // MeV/c
+ double h = 4.14e-15 * 1.0e-6; // MeV.s
+ double lambda = (h / p);
+
+ double EFinal = Beam.GetBindingEnergy() - Final.GetBindingEnergy();
+
+ // Width of the entrance channel
+ double x_Gp1[1] = {E};
+ double par_Gp1[7] = {Er,
+ res_width,
+ li,
+ (double)Target.GetZ(),
+ (double)Target.GetA(),
+ (double)Beam.GetZ(),
+ (double)Beam.GetA()};
+ double Gp1 = f_proton_width(x_Gp1, par_Gp1);
+
+ // Width of the direct 2p channel
+ double Gp2 = 0.;
+ if (E - EFinal > 0.) {
+ double x_Gp2[1] = {E - EFinal};
+ // 2*Target because I considered that 2p = 1H + 1H, so Z=2, A=2
+ double par_Gp2[6] = {Er - EFinal,
+ lf,
+ 2 * (double)Target.GetZ(),
+ 2 * (double)Target.GetA(),
+ (double)Final.GetZ(),
+ (double)Final.GetA()};
+ Gp2 = f_theoretical_proton_width(x_Gp2, par_Gp2);
+ }
+
+ // Width of the sequential 2p channel
+ double EIntermediate = Intermediate.GetExcitationEnergy();
+ Intermediate.SetExcitationEnergy(IntermediateState);
+ double Gp3 = 1.e-3;
+ if (E - EIntermediate < 0.) {
+ Gp3 = 0.;
+ } else if (E - EIntermediate > 0.) {
+ double x_Gp3[1] = {E - EIntermediate};
+ double par_Gp3[6] = {Er - EIntermediate, lf,
+ (double)Target.GetZ(), (double)Target.GetA(),
+ (double)Beam.GetZ(), (double)Beam.GetA()};
+ Gp3 = f_theoretical_proton_width(x_Gp3, par_Gp3);
+ }
+
+ double Gt = Gp1 + Gp2 + Gp3;
+
+ double mBeam_uma = Beam.GetA() + Beam.GetMassExcess() / 1000.;
+ double mTarget_uma = Target.GetA() + Target.GetMassExcess() / 1000.;
+ double mCN_uma = (mBeam_uma * mTarget_uma) / (mBeam_uma + mTarget_uma);
+
+ double CN_formation = 0.657 / (mCN_uma * E) * w * (Gp1 * Gt)
+ / (pow((E - Er), 2.) + 1. / 4. * pow(Gt, 2.));
+ double CN_decay = Gp3 / Gt;
+ double BW_cross_section = CN_formation * CN_decay;
+ return BW_cross_section;
+ }
+} // Breit-Wigner cross section of the reaction 14O(p,p)14O(p)13N
+
+// "Main of thsd
+void ThreeBodyCS(double Er = 6.3, double res_width = 0.02, int Npx = 1000) {
+
+ // Generating BW cross section for the 15F sequential emission, here 14O is
+ // taken as the reference, E(14Og.s.) = 0 MeV, as a consequence, the ground
+ // state of 13N is 4.628 MeV.
+
+ double J1 = 0.;
+ double J2 = 5. / 2.;
+ double li = 3.;
+ double lf = 2.;
+
+ TF1* g_double_cross_section
+ = new TF1("g_double_sequential_cross_section",
+ f_double_sequential_proton_cross_section, 0.1, 10.1, 6, 1);
+
+ g_double_cross_section->SetParameters(Er, res_width, J1, J2, li, lf);
+ g_double_cross_section->SetNpx(Npx);
+ g_double_cross_section->Draw();
+ g_double_cross_section->SetName("BW");
+ TFile* fout = new TFile("BW.root", "recreate");
+ g_double_cross_section->Write();
+}
diff --git a/Projects/MultiParticleRES/utilities/rootlogon.C b/Projects/MultiParticleRES/utilities/rootlogon.C
new file mode 100644
index 0000000000000000000000000000000000000000..9f03d73f93c0f366dcbc03c47ac09054622d9b99
--- /dev/null
+++ b/Projects/MultiParticleRES/utilities/rootlogon.C
@@ -0,0 +1 @@
+{ int status = gSystem->Load("libgsl"); }