From ab28cb07f281bc3da8dd19ea245e533b4c381eff Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Val=C3=A9rian=20Alcindor?= <valerian.alcindor@hotmail.fr>
Date: Mon, 15 Aug 2022 09:49:54 +0200
Subject: [PATCH] Adding a three body decay process for sequential emission
from resonant states
---
NPSimulation/Process/ThreeBody.cc | 532 ++++++++++++++++++
NPSimulation/Process/ThreeBody.hh | 107 ++++
Projects/MultiParticleRES/.DS_Store | Bin 0 -> 8196 bytes
Projects/MultiParticleRES/Inputs/.DS_Store | Bin 0 -> 6148 bytes
Projects/MultiParticleRES/Inputs/BW.root | Bin 0 -> 6883 bytes
.../Inputs/ExampleSequential15F.detector | 75 +++
.../Inputs/ExampleSequential15F.reaction | 38 ++
Projects/MultiParticleRES/Inputs/simu.mac | 2 +
Projects/MultiParticleRES/simu.sh | 2 +
Projects/MultiParticleRES/utilities/BW.root | Bin 0 -> 6883 bytes
.../MultiParticleRES/utilities/ThreeBodyCS.cc | 267 +++++++++
.../MultiParticleRES/utilities/rootlogon.C | 1 +
12 files changed, 1024 insertions(+)
create mode 100644 NPSimulation/Process/ThreeBody.cc
create mode 100644 NPSimulation/Process/ThreeBody.hh
create mode 100644 Projects/MultiParticleRES/.DS_Store
create mode 100644 Projects/MultiParticleRES/Inputs/.DS_Store
create mode 100644 Projects/MultiParticleRES/Inputs/BW.root
create mode 100755 Projects/MultiParticleRES/Inputs/ExampleSequential15F.detector
create mode 100644 Projects/MultiParticleRES/Inputs/ExampleSequential15F.reaction
create mode 100644 Projects/MultiParticleRES/Inputs/simu.mac
create mode 100755 Projects/MultiParticleRES/simu.sh
create mode 100644 Projects/MultiParticleRES/utilities/BW.root
create mode 100644 Projects/MultiParticleRES/utilities/ThreeBodyCS.cc
create mode 100644 Projects/MultiParticleRES/utilities/rootlogon.C
diff --git a/NPSimulation/Process/ThreeBody.cc b/NPSimulation/Process/ThreeBody.cc
new file mode 100644
index 000000000..f8415e749
--- /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 000000000..91a3ec1be
--- /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
GIT binary patch
literal 8196
zcmeHLU2GIp6h3EK;LenvOQp0pE?uleD20|ETLF=!KSBj9vi+yJ%<haZGM%YAv%4T^
ztZ#}0G`<=C;)@y+i4P?BWa3X^G~!Phe-eB!Mtw5z#Rp#0bLY+yyQL34ASTRB?mctP
zJ?GvzXTCFc&Hw;*70k^54FJHX3(9L~xJDE6qR(kcB&jBfWDgL4pDp>5S;x!H(T)Xi
z1mXz95r`uYM<9;C{~-c&W{bw=IrrskJjM}-BXCnjK)fH4)CEmNbXwBz)j>m80uU7?
zK-lP->VVKD8PQZkrzH(dX-!!@U}TC>iGk9b=95C5Xey%9l1g(xX$~0ej8TPx;qGLY
z4AlWslEz~kfj9zlBOuaeA@qR@lTfIr-!En-9orAu+rL4js(Qh~MQjnPVGpN|d8Kq<
z1O>m#$R8Hmfaf{c^0kp2vaM2vujuyNz_48_9~#)EMJnUNwrhH&F0bI4erTJP6HH}l
zhBur#wRQ96)?{aA$5d<b)b=gyt;ua$JEo>oR@c<F^HAovGvRutr5hN%3SdcvThCVd
z<oLyI`A?!sDf{$Y*(X)2(FccywG6K<(v9~|2Daz=Lss4=n&HcQa<DJwxd-!>m><r0
zPGQ3JGkkf@aB>C52&_GhleV9*{E@&q5s0fX-^&Nxq2h>5wnCQSyx?1hEHbc*Ru7r^
z9MyLg*ClRWyP^HA9rt$cJ#cVl$<k$9tJPKv7j56pI@VyxHiKiMhM%)s({_)JQMlZ|
ze$=+|>T+yV2V7&qn&B(!8&=(_>y_H&RYnR~`|*Muoa9Z6$vqnNH7dWr#!b4B2Ru`(
zXofR-qu$KeVIl9%TQo|swoY$Rxxof}@<lGOQEOAVuMQliXDUs%XxkV!)j`L|T8<QT
zYTFs#uMQVwK&7x&+r`w8oHsE+&iGQ6(z_V{ggqiQK4ADk|D@&8-omw>MqM8j@;#<$
zndCvk%6hK%>H8V~nvMC#jJ(xfBA}J!{fvbVj^tvcUegY|tCMLl$!)i26Nc##ODMCZ
z)H(^HFrrPVcFC!9F0F!1bPDZ*aj+;$PQn>D3+Lb!cmvME1$ZAmh0ov;d;v4?1N;QP
z!5{D^0wO9{gNtz~F2mcf8Q0>SxE1fgUAP-l*ozNhKR$$qa2SuEfjKnsG4wIOB9`z;
zd>Ws@)A$0ugfHW(_!^$ai}(S4h#%qScnQD6ukbSdjKAQoctuf^Dy3SfR~nQx%37sO
z-ndfEj#0^%M@2X`E8{6UqCa-Z=(VyV(bIe14a$!17mIvYUVF=my46kVHg2ITxr!PI
zvNEpziPsRKvw$5!G<>j?@REt6B?bw@)`|J@hc6Pj#%6?eV^T}-1%z_R-qB3#LPEG?
z*S91RT%p90?54Ja#;Yi8)X-LKNodh!N2jLoYDynDvr{J)B4&4r77@w7{yi*ThqvKf
zcn>}%Sbhta;Yau#uE1ZofM8gUs|bn-Tuo41j~lQRH)03gjXQ8B0kaqPVjtd*`*8q=
za0EwjEQF<nE>2=T1m_e%a~hu_Xr95dcn+V#=kY}X=$jRw+gFC0r1(}P=!)sQ>v--_
zl9z!x+a%_$djzWF4h^FI@16hu|8>8Qc!M|saRhE$1hBd%-P1+;&-W%4wRV#F5$d7|
yy=h596B@#DoM>2%6OGF8P&Y}c!X_2bX-P>!?LYq!F#Gc|KL6wMKg8k9HvA1LJqxD*
literal 0
HcmV?d00001
diff --git a/Projects/MultiParticleRES/Inputs/.DS_Store b/Projects/MultiParticleRES/Inputs/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6
GIT binary patch
literal 6148
zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3
zem<@ulZcFPQ@L2!n>{z**<q8>++&mCkOWA81W14cNZ<zv;LbK1Poaz?KmsK2CSc!(
z0ynLxE!0092;Krf2c+FF_Fe*7ECH>lEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ
zLs35+`xjp>T0<F0fCPF1$Cyrb|F7^5{eNG?83~ZUUlGt@xh*qZDeu<Z%US-OSsOPv
j)R!Z4KLME7ReXlK;d!wEw5GODWMKRea10D2@KpjYNUI8I
literal 0
HcmV?d00001
diff --git a/Projects/MultiParticleRES/Inputs/BW.root b/Projects/MultiParticleRES/Inputs/BW.root
new file mode 100644
index 0000000000000000000000000000000000000000..2838ace698508d1f053ff30ce8ec6fab773189ca
GIT binary patch
literal 6883
zcmbt(Wl&r}x9#8-AUMGZ1cJMU;O+zn4uiWp1c%@*!3n`(Ah<gOouGp|gAA_0<&xa0
zS6{vN<K3!zs#bT`KHYomUtPP`a&dHY0{|{#004j)06+qLLG{a;>jiKx;6!{G%mDx#
zasU901_0O0oLr0$g>w`MLm)&)?>u?=`#<^|0Dlun*{{4I1~B}~{K5tRAZbWj+nJ+?
zYq7fgt>gbJ07!p#zJTOqlhX@0U-0k5OMRmMzcOAc{AHu&{KqEok8SGjc?1n<4tQ}b
zYD*(CM|TrDb0b%CXLoZ4H)~@%BU2YgSJ#(-shhQ<gRUa#Gd6&e8sG(Yl<Cc%y0kRd
z#8Ov#?VtOxIyO<(ET;crC8bPuI3-k;TD`C|xfJ5SpPD4O?UqD~$(bq{nZh7B5G7qA
z7n1a;qs`U7>%i~X>-qNaPiH{$oO`AFPhaf=f8l4)flrm5=rsHHVGJv7n#^>o2+3^e
z<l-eCYtXq74gnD{#UUkWI~^;(i#_bz>G^wl>AZF#wUxb<M&Wz8Olcxa4ibYxg(&V|
zX*}!NfTjKf$%2iMM582z3C>$RtFBF}E>>L12TOwA&k|K7cTXMn*NopMEvC#Qi;Z-i
z6H%X3IA<RGHJP49#M7TIf<m7MpG&vS)}7&@hy!Rbc*OMV%E0Vy{-FO1tUnXLxi;dS
zO#dPMnpWErERf%@WN9nNZ|I|OYIPMG8kL8`^a14N*PvtM#su0VR62^mS8QnY@{}}P
z8xXDVzs23_&rzGjWk5k0UubWbZM^-m(yX$&cvUufj?idZ`OQINeGKc~xPwjvN^r)!
z=Cf1<x0mqKc{aRLV03nfqE7hF*ttv7tpCAg^3uwxeSob2dC5&7s1b5^oVbnE$#$i8
zJp_*T%dj(lte4+0X<B7G@I1H1gD5ZQ1!PR_us4H3aBa~8GNw55k_2cbkv?EV@%+;8
z7-n&ncb?6-V&bS{h!x0G7WAZ0Ro`u>7;Mw{#6`O}yVPP|tk_rY<|Po4vMt8-gss9h
zBF$o}O~a8;d(h@&HSbNt(-o-cBQlntyA|Lw`8Z?A5I>F9KWXMUQnJ^}+!jBg(LGSz
z+*sWKLAO6TF#6G@P93S%Hx-7Y(SO#iyV8dX$!u&P!b_W1b6(ZtvvV_4Lv!}M<3GnP
zlRPM&ZpkpvIG<K(sY;C1KtI{&qZUJY2Rg}k@GVC8)3z+BRj}t;_LyfL`bi?WG=%Tr
z;$!a{jXRlwI6vn&4U2)*UzS$3rO>(Is!nS-NThAyUSt?oS@tP+?aq?UzS|cfN9zQU
z)MpE&(bG2B{V$fHE6_}=^3!oV8H7c>GLzqe7M+(0<(H!)j>X?l%fBps<BI3BDYbR4
z@~gXwR-0P`EhgJpORm)Go2?6Nu$5osYRlZmBtwo!J+u~gi$2A^3P1Xc9$nfN+e8pv
z9HnM|866D8^DVntPq20zCMi!Y2?%Ww$lag(+(?NAE2bdJDZP1$Ae}M7%bMnVm*n%H
zZf7p<it)!zW;BJ*hOe}zTHblIlW|{j7_+O-2OXjs?76nmc3*m};`v-mv@eZ$zenpk
z#ZrcIA-Qn=E)6ob^&PF~V7inNqby5>(pwD7)|o?<MD>*Dpg)|6Jj|vOZttT`)z+~(
zG9>gWxvHGsaC#i^qGwFRcS}a}iuLy_e?3mSyDqv47LVAi|5-M+q3MJw*zd}~2}Req
zY#gDV{CizM&w&sF1FUI^ZEUxwciveivht1>+^(&4OJjxvy3#QK<q;$4-KW+RAQw2n
zSim?RR^jG|+KF_h>3X2a^Os6rnbt*&Jh_w{2tNO|(>W&fd}X_O=rZAkH|WS^`&*a%
zuBOs_tmmm%5EWs@{)+giztJIVRLSPoyhf*p9hxP>Y_Zd2fd?^_xbI*Bh@Gq_>tc>a
z5|IWdBHQXf_61&A9LuY%2U;zq>XLxIJOspi+zHp7XtE4eTVn!O<TGs-iv)%^xwxMx
zv4Q=vXt((@u$OjP1%FPM%nuDdPO54mDel?1X;6gCl;wNTU%#nd#)S+;ss{ebXq3pW
z!v_kNi3zYXqgb>wUClbz14&3W>^<yp6$v?YdHMLd=lkRi#tFuO^^NQ0X5+%dJCtR%
z^<#Z8V(&~Gp-9AJS?dH&W9WBkMd5d(9ULos6%ktJ03W!LR03j`^fcg%Qu{sJ3fYXN
z!qn3D3#+o}1$0svY;2|0;r_Ft68hz2S)s)J@V;93lEs7}ZA91S06cXE?9yLHkYD&E
zry*;GP2v?wWsf~gx>iVIt)55XZ@FbkipEcDgDhH&j$|?MGV6rurtse@QVi`79w`oP
zb+C}``U%Dj2Y<}_Hs~!1z1Xp5S{;#q{uoclWy2)FF`yKj_RTkTGFMkm=7NZ5DCA9J
z1mC6FUUtor2j0zwN~!B%<{$f%?DnmYl~xohGLSxCHOO%0#b?M!u>4wxp&<@ljPXP8
zi5XP&OHhOSNk^LTLq2kS_2avs7SE0)M6)y|LuSUZJxD$y+704Dt+5Hd;;FntC#bvI
zqYYJkIu~MkZ%UoQB;kL0eVX6C^8TAH7hpRq*~lAPQdCDLgu+q;^Whu#_#&RCyK9P%
zGk=rfrcLYKNk(A(J_l<)wGHu0@)|m0;~J7tvF#^+50MmPWtq>fCi0n8;yqh+QJ#`b
zF5edg8L8EHMLN^ogm0Uik9!7Fx=sMKYYDz+wRN75Ms?5aqQ>(S1+!gmv;8`KRQ2Im
z-xC=K{A8*`u;h-HGG=MkM(K-Qq|wAyP93JIdp||(bcmXg=0!_h>F^b|L1y!)^6h);
zqU=&HCv16_G&w2q->MHzF0H?~KD@uuSFSnz8M^s<o^e-H;5{Mteda8^I7f$%jm39t
z{Nx<P=73oVjxRcT=OqZLmv@;sIxXZ;uBppkDF#F*Pirq{Ygtq^te?zy87QM%gg+Qc
zeh%W>-&?!cgY#lz+?|Y;?{4cN!jKsY)Y2p8SMx(`#~tS=daHG+EZ~qdjJtqM9^={3
ziq?KqOkcNp9o%GCB5~?R>P3xnITii5&T<Ra7pc&bwPqlEZ3^c#Ynr+n8y{@h{3BuJ
zfka3tOm2LZCuoo&R=Y42H8)1Oy-bp(*>q`D#Xbq;I%nYa@^EQg_%>zX9*nU(MkzJf
z$hgismiBpfF2wia&t;ufq-EST4roDi&u>f49a*c8ES4i;HbM-CMefM+XpEkfmPT-k
zCdP9A?#Hw#+irubo9*KxQKvqO-HYV3AkuD?SWpV$oVxkKdKTBZV1XegC%XR97W1^i
z<ht&wZ0*WN;^jB~MNKM40M8NG)Qm`i<`K(_!bRzxcetMovbQm+K+o3%cd6@9w^B`^
z7J(P<8XIbQ2($C;FgAMu{Nt3G@hmX--TXbAi~PKGH!^a?BOQ2UylWlfIb|ms7ZPOf
zR!>>=k*#_HQVGmYbf|*UY0Pj4=c`Vj>v>a+Onv0KnYY%rTtT=oF}-WOQJe0CkI_tt
zg@Cco73qbWw(EA;YsSs@H2gKK?^uLzv#ZaGnO3p8g&}IZm1P8n9ZLg?B)LHyeLII7
zWv)Q_VF9wTd32A5z-u`Rig{ayuf`O3ht=XO_VvWoA!FPDn4eLWu9w^GIo<(%^|!~=
zys)l<8)uaY3Ffi-+R4?SJ&kHBdat2yjoUvQD8^sgC9l<0nAr8?uHLnchH3(BE1E}a
z8G<Sab7nbCKFqfx4>ZXk0qcqfv%c8`vwTa-CUC*9L%X}clsKje0AXL2goluM2Uu%A
zat>YFVg-=9A>N)}e2nIEr^=%{2pwKSxJ{@5Yj+WuS&lJ>#gbF$ResVh45d0~+ZJhv
zeBeQoaa<}525At?Ju(rSEwCVWT@eiR@`S%)yb#UR#Ex4#!sIOZEpapFmoVrFzRH`z
z<_z&U5V<5Y+b`{V`2F_&_>;T(B54AH^Rev_v#3L^_s8o@)}fo9jWQTnVgnimXQj+Y
zZ#qqThgj7*n26+O4OfsjMfXCJK64Mf8y`(mx<+bg)_#_m$4Ow-GnDSn3UYI58zlLR
zY+c|zw|6d=paB}boW>B-4V!U1n*dms8L&Rmq8ECvuIzd=P$`LU`!gyaSTnz|(=J1p
zTF(dsufTtDe1?S{Zd18$LC&wRf+6SCPVF_Mg{^gBX>9_s{@b%vLJSG|f-Dyfyt0HU
zta_X9Ly>MZeEV;^AegQVT4a;Zz1uxd4R%CL=5L#_S}N(Stv+u$s?`IYidSaAUOa#P
z)QDIFdP(C|5k_oMq^0?8cc4u2JOh6TV85=M(~L%Jpdy*~$10c;gYdKQ^JOD!>YmnF
zDJi`!5(DF9t5dOpD-)dYf`^hTT=}Yy&jtww4z<;uaK1(K68KPkI>8^QCkmW5PpD8+
z3Y@Qvl(qInG?9?I-QXG&?RWDv*qjMMC4HD)7<ePOu1@5yjEI^?C&y*#`r3hh?%;kE
znS-r=Y}w=ks-AhVeg4}B9nG_-4ne9I!aCltf#8X~7`$dMeEzWsescjPVpFrh_$4N4
zm>iqk@ozZlmdlmOS+(j;%2)&P^Vgt_w5HyJACDE2;asM7-9Eh?B(I%FK)KdyKH*<`
zY<SEcJ$4Yao0NXI?H~qf7G>Qed7^?oDF@LG;(~%5{awm+fH{X`eu6hxAaHOHADl|J
z{L1TVzg?8-P~QIHZ*4vDx>{V<f!Nw;7crKEUC2ydMNHVr;|Z>G-Ca`N3zxR5YJY%Z
zULz;C8f*;Ou9>&qLBp;lEJfVdBd9*it6L$p0&eKfXL{lSb^4{r)MEgnt)Zt7d4R2j
z$t2RF0r0m|g2=W$cG?)N0VUg;YRJ;C3YmK~K9|oEsx~^|D4+Sa1)l0H2!v9k3|XPQ
z3p5S{fNSt5yoV)#{<^F+VqEy(CL4RbA4ovq0M%<l9b~Xf#))v1FS^CIZ=RFe=D<eW
z+SRH(;HtGnl%77IPPe&p=7KR@%RZ;|gds%d@+i=s69pW+S#*7*-Xla{=w)<60fy=>
z2o9bEW~^GOT6+?L{ZO>%85|MS`PsUq3b4Rt9LIAR%0Rjlay>T}mUIrBHSLy9z^d8I
zM+Yzl_<L@3y_YAu-*t%vnLP4z-_t%tnn*gUDdp<pF<P--V&yizEjx~`o0b%oB8EF(
zS)~t<#6vk}bqTxp*zh29N1UC0$(+eUFtF@$Mcpz5j!q5Nu-*YT{T%0~LMSUx-s<{7
zGYt6U#7H}J8*5-!aWJCq9QdlGNeDcQV5+(S)m)Wgzat%qU*<w7(BluGb3~WU${BZw
z#RhUmELD8H4R*d5Xd+Ll0&37Y9x-143i>T3YVBx<cYnq1W3j@~H1WDs6v7_tjK3Pq
zv9NPyA6OW8BGaVaLK9<Q**OMrpY{&A=ac+qZAB1bVyUI!^jHIHW?^TPcXYte_?dS(
zQvjNqtMAv@ZtM?J2IAAG?CqZ)yf~W>V!|x&?z}sI{Bl*g3mEKnt(d*bTEWKN4!7h%
zFvH|ZPK-mTHGvh&yRdkFpdv?S)AN)Xy9T{HEiVW@gv-m)Yk`*?-(k<|b}oe4&t~81
zik!Xwldji{=XdtyP2HMtBJ^>t-ijgF_imL5=^dU6&r>3&YtL|q+`3yTiOox-Ecz66
zLT|n7ofHp$A|pD!ys4kGpn;-id9w(G_QX#7Dh>~xoxom^Owl`|gDx_7lGtmNQMiwb
z2^98{@w~+{SrY|Hdy^Jj;!hY{;(fLMG>;c4;Z5B?DdI(V_-`cyV1}fjXzl9uO5M%H
z+}PgSMb5#(5m)lBsNhXyXYF8XZuT;{y1nQKi>QFQH!nKEZ=DrKGA#m;4~?dS9<pRf
zQ*@O++l{o1oA0GqaFv#0p#8rpO|=$&%b?k(uXLEIaDJ%E$oPZOut)eFK1x+jQ#_am
z5)dfkCHa6C=(UbRk#N^bioM}Eww#gUZf~g+ysmrE;(2`HIWBO$>>6k9v(LvZ`cAy=
zCe&l1c<O}WyXJ^oVb0gM9Rxnw_Q#lj0#g_(<>khg@3j#W?so@{MqPo93a|!-HCD}&
zG0DWbp0$;tPKOBLfgaUL)^r7B(CrcVe&4O&ALzpjSxf{n6S3x8636AG&Hm!lC#kS=
zWD&#O$r)k6+xxQ3=qClPZ6+_cxD+=y(_F)9@WdJ!R;};yC?rrb5BCC;#1FAE9
z6&XqH5U3nkTca2uEx_)_Nqw0c#hjkm#up_g4V>LTy)#^+c#`j-lvb3)>HV-zSuozJ
zpOLKt=NtJ|og>T9_#dy=Rft{j4RU4XjMnciJurtu?S5QOY+yH6;T#lFhZ_BsDO?87
zF)`_C=BH`57#ETeRO=8Ql}<mnNe|o6olY{xRat=P1GFHiOUexGH293-8)0w5@<$Yt
zn7x98!Svg-2_6lpFn-D^r;EMy)Z1PfbJ;W{x1Tw0XAt*T#4;b?r?PkN2X6|=nB&8{
zKi*ty(YN)i@HlA(wY_t~lWI23o!6vqEh(tEJbzsuX({|L=wNl?p3x5aJ|eA7sGq}}
z1tX_s*Q8}%U$yb-i41d6k|`%?4E8BsMI4FC24*O^4v{q<E@MIZKF`eYCax7f)iobS
z?t=1~SM7ue+cSIGln7PT!(GaK{iBuW1smsWG~kF^(N(G;k}`H54vX%b4;Sy+scIt0
z5H1LHych!}MK^@@4WTf7rqXwY#P5Y*R0HIHu+AkAYUR?D>KUoln8%9N4D6}`I@nWR
zADYZW)x7eSs~t<DDQVf9Av2O{eX~bw2sW06e2CBbI+%2tPDEkEdQR;03*__01#uOY
zs#OI(LrI7Dan;jH%$dQyMN43tL{){$9SJG0CCL4dVOy8>W<gkdJz);_D2;>f(jnwo
z9-=)rnmEg$7ag7DGq3uL<J1SRh3lm0wv6uXJ3fpeT663f6-vy)qARgkuvK)WtYG5~
z9EFFWuD85hPw<_=XqT>V?$w_uhc?w~t4w1jhw=>%V(o!S688QiQU}az?~J4(t<DBW
zw)<9&rm$J4H?n^-zE+}d#g#I~&r5O&n<?MUIPv=3y^~!!%8~LE$#%Vb)>R+%{#^`f
zVeW~h+u*8#1n=aay`UG_TpntHXV_%X+X68Q#9NhZw=ZX%opZp^XQL#S*$;)fn?Ihj
z3Rn4(7GKM);K}(_zAuAoIt1-ox(c*JbL6IGdFyZ(-kiy(&86vZvDALASCyTAgSJX3
zz5I*ahQIULz(6E*B}zrpSUvmEdhX1$Uz)VM%AjPMo(w?#tB0<}_TcB&d`G-T$m<(n
z9j`@vct?CRGF<ozG5Dsfo88*1WmQgV4IW;&i3#m~#vy^}&e;n~gJ>qqfwe6^<REs~
zV<2;OQ|Q)#RFcFo*=xPdi=@5EBvnoD2ai0S&?g&vDU^02DQ1Sat?qP@x5EVeGaDYe
zyeY93G@(mD+ltj{RCl&=@cs;v6;eNH@3N)8G<0)acjT;7#oPXACAuho^?>P|ES9eL
zX10ZT0CdX4Kd*JNh-J*yKf*T|j)^^4Jc&AYJ=`ns2fyXSWSZ37vQL^T<|A(%5qhX?
z&fNZ`$?ctH*inFo^%*i=kG3z~8-&CR8+~HZwG`gPB=wYI1?4kIRxV;^3jgHDN|wt!
znE+bCqKu<RBrZ$bp-@^vW(#rXZ8ePR(1lxs^v;-@+1W4F@hrp~OW(r*;Y^;`88}NL
zEwFIv9))wkBJFob%-7*ln`8nyvp*qbezaUFyX>lsz7!y(PPWnJo!o_C?@g~YD&*~k
z!O4;Cbtw_%I*f?JqGG}gDAWrBnv)L>jgUUVN1~&P-a`;ZVIdxj1uTh^$t?^$;uzle
z8Te?^tf^;6OEA)<qLc-z71=Uv0vMxgOslmnoeK#M!{SpIgexwBp>Lri#_mUX#(2|E
zuX?7Or>1~3_{+%rm^9eUWEn4CR??SV0-=gh5RH5CSxxU{faZW^`28|-Z^$z9C1CO9
zBR;UM1_SI(_M_+Mw}|j?kxRM|6mjd!1et^O^0?+Y;l|lG43VU@VvZvk>9>@P>H#{-
zrbMnZC_$#jF(<p}&N}pwRe2c~Ip(W&w#EsP`9FoU1WNU}FLJ&{6Emib5f|<@Z>iJ7
z=S2#2(qqx8&}maFDEqj$s<Qx-w*m_hn6~qCY0vfhtPR<b>bWAg$IM7O6x|uh4<ewp
zW78$}KRiCLR97a4mB@M-M?$e<=iDix6pN)=l<B}hiHHx@ko-%*k|8>!yq*%@uMc5w
z*ILyylV7gJDyL@Et%O+wI#x?_08bITSrVPvb&hnyG5$I=eg2s1azlpO6+*!`mPrMY
z>dqA7Izwr5J@>Z}Y>_I<2|{7+2+)}!*IHN${R$a}bsWQ}vPI;WR2D4%N^rAyolKzX
zi>s8xWuk5fOI0XK*O>En`S4vVvK~pFH$A}68GLUnbb02q)k6a1NW#U<6j&{P5i~=#
z`D;#DTxp6>m>5eIEW@KUgktSDS(Zd`pXahU_CV^*shM;4v3A_kyyYXAQR>9)d*MwM
zORxKC5J@rG5{}3;{-fqEkMm3_Pm+y?3CT#AbMpsn{a5nT!7Hy1>)$3iddB!A)U;qx
zP##mcS~t<#4kvADM{Q(ODw@s~$wkgDD4kAE=RC{Kb8kYAl6)&u8Og*81LPddD05mL
z*Ids+c8qW#$;E0}l@<JUQam3IJ~#_~jO+H$f2UjzO5V)U)VO~wrV_vAf3~~?-ND;i
zAxBmzH|h2($&g&Inuh)|8rZNq@bP(Asd}<(Sx(EwLYTYwvtSZHVZABj5<woIHqB^t
zD}0x&r3a0ItwAwB&;R_Hmlg>i^RGqv=k2^WG%s$;|K;txgkHL*|Ft#ql7{P_mxTTw
Te19FQ|A+5|jRfk+lLGiJiN-M7
literal 0
HcmV?d00001
diff --git a/Projects/MultiParticleRES/Inputs/ExampleSequential15F.detector b/Projects/MultiParticleRES/Inputs/ExampleSequential15F.detector
new file mode 100755
index 000000000..bfc19a7ed
--- /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 000000000..de9c8857e
--- /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 000000000..add2e059e
--- /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 000000000..b85e7e901
--- /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
GIT binary patch
literal 6883
zcmbt(Wl&r}x9#8-AUMGZ1cJMU;O+zn4uiWp1c%@*!3n`(Ah<gOouGp|gAA_0<&xa0
zS6{vN<K3!zs#bT`KHYomUtPP`a&dHY0{|{#004j)06+qLLG{a;>jiKx;6!{G%mDx#
zasU901_0O0oLr0$g>w`MLm)&)?>u?=`#<^|0Dlun*{{4I1~B}~{K5tRAZbWj+nJ+?
zYq7fgt>gbJ07!p#zJTOqlhX@0U-0k5OMRmMzcOAc{AHu&{KqEok8SGjc?1n<4tQ}b
zYD*(CM|TrDb0b%CXLoZ4H)~@%BU2YgSJ#(-shhQ<gRUa#Gd6&e8sG(Yl<Cc%y0kRd
z#8Ov#?VtOxIyO<(ET;crC8bPuI3-k;TD`C|xfJ5SpPD4O?UqD~$(bq{nZh7B5G7qA
z7n1a;qs`U7>%i~X>-qNaPiH{$oO`AFPhaf=f8l4)flrm5=rsHHVGJv7n#^>o2+3^e
z<l-eCYtXq74gnD{#UUkWI~^;(i#_bz>G^wl>AZF#wUxb<M&Wz8Olcxa4ibYxg(&V|
zX*}!NfTjKf$%2iMM582z3C>$RtFBF}E>>L12TOwA&k|K7cTXMn*NopMEvC#Qi;Z-i
z6H%X3IA<RGHJP49#M7TIf<m7MpG&vS)}7&@hy!Rbc*OMV%E0Vy{-FO1tUnXLxi;dS
zO#dPMnpWErERf%@WN9nNZ|I|OYIPMG8kL8`^a14N*PvtM#su0VR62^mS8QnY@{}}P
z8xXDVzs23_&rzGjWk5k0UubWbZM^-m(yX$&cvUufj?idZ`OQINeGKc~xPwjvN^r)!
z=Cf1<x0mqKc{aRLV03nfqE7hF*ttv7tpCAg^3uwxeSob2dC5&7s1b5^oVbnE$#$i8
zJp_*T%dj(lte4+0X<B7G@I1H1gD5ZQ1!PR_us4H3aBa~8GNw55k_2cbkv?EV@%+;8
z7-n&ncb?6-V&bS{h!x0G7WAZ0Ro`u>7;Mw{#6`O}yVPP|tk_rY<|Po4vMt8-gss9h
zBF$o}O~a8;d(h@&HSbNt(-o-cBQlntyA|Lw`8Z?A5I>F9KWXMUQnJ^}+!jBg(LGSz
z+*sWKLAO6TF#6G@P93S%Hx-7Y(SO#iyV8dX$!u&P!b_W1b6(ZtvvV_4Lv!}M<3GnP
zlRPM&ZpkpvIG<K(sY;C1KtI{&qZUJY2Rg}k@GVC8)3z+BRj}t;_LyfL`bi?WG=%Tr
z;$!a{jXRlwI6vn&4U2)*UzS$3rO>(Is!nS-NThAyUSt?oS@tP+?aq?UzS|cfN9zQU
z)MpE&(bG2B{V$fHE6_}=^3!oV8H7c>GLzqe7M+(0<(H!)j>X?l%fBps<BI3BDYbR4
z@~gXwR-0P`EhgJpORm)Go2?6Nu$5osYRlZmBtwo!J+u~gi$2A^3P1Xc9$nfN+e8pv
z9HnM|866D8^DVntPq20zCMi!Y2?%Ww$lag(+(?NAE2bdJDZP1$Ae}M7%bMnVm*n%H
zZf7p<it)!zW;BJ*hOe}zTHblIlW|{j7_+O-2OXjs?76nmc3*m};`v-mv@eZ$zenpk
z#ZrcIA-Qn=E)6ob^&PF~V7inNqby5>(pwD7)|o?<MD>*Dpg)|6Jj|vOZttT`)z+~(
zG9>gWxvHGsaC#i^qGwFRcS}a}iuLy_e?3mSyDqv47LVAi|5-M+q3MJw*zd}~2}Req
zY#gDV{CizM&w&sF1FUI^ZEUxwciveivht1>+^(&4OJjxvy3#QK<q;$4-KW+RAQw2n
zSim?RR^jG|+KF_h>3X2a^Os6rnbt*&Jh_w{2tNO|(>W&fd}X_O=rZAkH|WS^`&*a%
zuBOs_tmmm%5EWs@{)+giztJIVRLSPoyhf*p9hxP>Y_Zd2fd?^_xbI*Bh@Gq_>tc>a
z5|IWdBHQXf_61&A9LuY%2U;zq>XLxIJOspi+zHp7XtE4eTVn!O<TGs-iv)%^xwxMx
zv4Q=vXt((@u$OjP1%FPM%nuDdPO54mDel?1X;6gCl;wNTU%#nd#)S+;ss{ebXq3pW
z!v_kNi3zYXqgb>wUClbz14&3W>^<yp6$v?YdHMLd=lkRi#tFuO^^NQ0X5+%dJCtR%
z^<#Z8V(&~Gp-9AJS?dH&W9WBkMd5d(9ULos6%ktJ03W!LR03j`^fcg%Qu{sJ3fYXN
z!qn3D3#+o}1$0svY;2|0;r_Ft68hz2S)s)J@V;93lEs7}ZA91S06cXE?9yLHkYD&E
zry*;GP2v?wWsf~gx>iVIt)55XZ@FbkipEcDgDhH&j$|?MGV6rurtse@QVi`79w`oP
zb+C}``U%Dj2Y<}_Hs~!1z1Xp5S{;#q{uoclWy2)FF`yKj_RTkTGFMkm=7NZ5DCA9J
z1mC6FUUtor2j0zwN~!B%<{$f%?DnmYl~xohGLSxCHOO%0#b?M!u>4wxp&<@ljPXP8
zi5XP&OHhOSNk^LTLq2kS_2avs7SE0)M6)y|LuSUZJxD$y+704Dt+5Hd;;FntC#bvI
zqYYJkIu~MkZ%UoQB;kL0eVX6C^8TAH7hpRq*~lAPQdCDLgu+q;^Whu#_#&RCyK9P%
zGk=rfrcLYKNk(A(J_l<)wGHu0@)|m0;~J7tvF#^+50MmPWtq>fCi0n8;yqh+QJ#`b
zF5edg8L8EHMLN^ogm0Uik9!7Fx=sMKYYDz+wRN75Ms?5aqQ>(S1+!gmv;8`KRQ2Im
z-xC=K{A8*`u;h-HGG=MkM(K-Qq|wAyP93JIdp||(bcmXg=0!_h>F^b|L1y!)^6h);
zqU=&HCv16_G&w2q->MHzF0H?~KD@uuSFSnz8M^s<o^e-H;5{Mteda8^I7f$%jm39t
z{Nx<P=73oVjxRcT=OqZLmv@;sIxXZ;uBppkDF#F*Pirq{Ygtq^te?zy87QM%gg+Qc
zeh%W>-&?!cgY#lz+?|Y;?{4cN!jKsY)Y2p8SMx(`#~tS=daHG+EZ~qdjJtqM9^={3
ziq?KqOkcNp9o%GCB5~?R>P3xnITii5&T<Ra7pc&bwPqlEZ3^c#Ynr+n8y{@h{3BuJ
zfka3tOm2LZCuoo&R=Y42H8)1Oy-bp(*>q`D#Xbq;I%nYa@^EQg_%>zX9*nU(MkzJf
z$hgismiBpfF2wia&t;ufq-EST4roDi&u>f49a*c8ES4i;HbM-CMefM+XpEkfmPT-k
zCdP9A?#Hw#+irubo9*KxQKvqO-HYV3AkuD?SWpV$oVxkKdKTBZV1XegC%XR97W1^i
z<ht&wZ0*WN;^jB~MNKM40M8NG)Qm`i<`K(_!bRzxcetMovbQm+K+o3%cd6@9w^B`^
z7J(P<8XIbQ2($C;FgAMu{Nt3G@hmX--TXbAi~PKGH!^a?BOQ2UylWlfIb|ms7ZPOf
zR!>>=k*#_HQVGmYbf|*UY0Pj4=c`Vj>v>a+Onv0KnYY%rTtT=oF}-WOQJe0CkI_tt
zg@Cco73qbWw(EA;YsSs@H2gKK?^uLzv#ZaGnO3p8g&}IZm1P8n9ZLg?B)LHyeLII7
zWv)Q_VF9wTd32A5z-u`Rig{ayuf`O3ht=XO_VvWoA!FPDn4eLWu9w^GIo<(%^|!~=
zys)l<8)uaY3Ffi-+R4?SJ&kHBdat2yjoUvQD8^sgC9l<0nAr8?uHLnchH3(BE1E}a
z8G<Sab7nbCKFqfx4>ZXk0qcqfv%c8`vwTa-CUC*9L%X}clsKje0AXL2goluM2Uu%A
zat>YFVg-=9A>N)}e2nIEr^=%{2pwKSxJ{@5Yj+WuS&lJ>#gbF$ResVh45d0~+ZJhv
zeBeQoaa<}525At?Ju(rSEwCVWT@eiR@`S%)yb#UR#Ex4#!sIOZEpapFmoVrFzRH`z
z<_z&U5V<5Y+b`{V`2F_&_>;T(B54AH^Rev_v#3L^_s8o@)}fo9jWQTnVgnimXQj+Y
zZ#qqThgj7*n26+O4OfsjMfXCJK64Mf8y`(mx<+bg)_#_m$4Ow-GnDSn3UYI58zlLR
zY+c|zw|6d=paB}boW>B-4V!U1n*dms8L&Rmq8ECvuIzd=P$`LU`!gyaSTnz|(=J1p
zTF(dsufTtDe1?S{Zd18$LC&wRf+6SCPVF_Mg{^gBX>9_s{@b%vLJSG|f-Dyfyt0HU
zta_X9Ly>MZeEV;^AegQVT4a;Zz1uxd4R%CL=5L#_S}N(Stv+u$s?`IYidSaAUOa#P
z)QDIFdP(C|5k_oMq^0?8cc4u2JOh6TV85=M(~L%Jpdy*~$10c;gYdKQ^JOD!>YmnF
zDJi`!5(DF9t5dOpD-)dYf`^hTT=}Yy&jtww4z<;uaK1(K68KPkI>8^QCkmW5PpD8+
z3Y@Qvl(qInG?9?I-QXG&?RWDv*qjMMC4HD)7<ePOu1@5yjEI^?C&y*#`r3hh?%;kE
znS-r=Y}w=ks-AhVeg4}B9nG_-4ne9I!aCltf#8X~7`$dMeEzWsescjPVpFrh_$4N4
zm>iqk@ozZlmdlmOS+(j;%2)&P^Vgt_w5HyJACDE2;asM7-9Eh?B(I%FK)KdyKH*<`
zY<SEcJ$4Yao0NXI?H~qf7G>Qed7^?oDF@LG;(~%5{awm+fH{X`eu6hxAaHOHADl|J
z{L1TVzg?8-P~QIHZ*4vDx>{V<f!Nw;7crKEUC2ydMNHVr;|Z>G-Ca`N3zxR5YJY%Z
zULz;C8f*;Ou9>&qLBp;lEJfVdBd9*it6L$p0&eKfXL{lSb^4{r)MEgnt)Zt7d4R2j
z$t2RF0r0m|g2=W$cG?)N0VUg;YRJ;C3YmK~K9|oEsx~^|D4+Sa1)l0H2!v9k3|XPQ
z3p5S{fNSt5yoV)#{<^F+VqEy(CL4RbA4ovq0M%<l9b~Xf#))v1FS^CIZ=RFe=D<eW
z+SRH(;HtGnl%77IPPe&p=7KR@%RZ;|gds%d@+i=s69pW+S#*7*-Xla{=w)<60fy=>
z2o9bEW~^GOT6+?L{ZO>%85|MS`PsUq3b4Rt9LIAR%0Rjlay>T}mUIrBHSLy9z^d8I
zM+Yzl_<L@3y_YAu-*t%vnLP4z-_t%tnn*gUDdp<pF<P--V&yizEjx~`o0b%oB8EF(
zS)~t<#6vk}bqTxp*zh29N1UC0$(+eUFtF@$Mcpz5j!q5Nu-*YT{T%0~LMSUx-s<{7
zGYt6U#7H}J8*5-!aWJCq9QdlGNeDcQV5+(S)m)Wgzat%qU*<w7(BluGb3~WU${BZw
z#RhUmELD8H4R*d5Xd+Ll0&37Y9x-143i>T3YVBx<cYnq1W3j@~H1WDs6v7_tjK3Pq
zv9NPyA6OW8BGaVaLK9<Q**OMrpY{&A=ac+qZAB1bVyUI!^jHIHW?^TPcXYte_?dS(
zQvjNqtMAv@ZtM?J2IAAG?CqZ)yf~W>V!|x&?z}sI{Bl*g3mEKnt(d*bTEWKN4!7h%
zFvH|ZPK-mTHGvh&yRdkFpdv?S)AN)Xy9T{HEiVW@gv-m)Yk`*?-(k<|b}oe4&t~81
zik!Xwldji{=XdtyP2HMtBJ^>t-ijgF_imL5=^dU6&r>3&YtL|q+`3yTiOox-Ecz66
zLT|n7ofHp$A|pD!ys4kGpn;-id9w(G_QX#7Dh>~xoxom^Owl`|gDx_7lGtmNQMiwb
z2^98{@w~+{SrY|Hdy^Jj;!hY{;(fLMG>;c4;Z5B?DdI(V_-`cyV1}fjXzl9uO5M%H
z+}PgSMb5#(5m)lBsNhXyXYF8XZuT;{y1nQKi>QFQH!nKEZ=DrKGA#m;4~?dS9<pRf
zQ*@O++l{o1oA0GqaFv#0p#8rpO|=$&%b?k(uXLEIaDJ%E$oPZOut)eFK1x+jQ#_am
z5)dfkCHa6C=(UbRk#N^bioM}Eww#gUZf~g+ysmrE;(2`HIWBO$>>6k9v(LvZ`cAy=
zCe&l1c<O}WyXJ^oVb0gM9Rxnw_Q#lj0#g_(<>khg@3j#W?so@{MqPo93a|!-HCD}&
zG0DWbp0$;tPKOBLfgaUL)^r7B(CrcVe&4O&ALzpjSxf{n6S3x8636AG&Hm!lC#kS=
zWD&#O$r)k6+xxQ3=qClPZ6+_cxD+=y(_F)9@WdJ!R;};yC?rrb5BCC;#1FAE9
z6&XqH5U3nkTca2uEx_)_Nqw0c#hjkm#up_g4V>LTy)#^+c#`j-lvb3)>HV-zSuozJ
zpOLKt=NtJ|og>T9_#dy=Rft{j4RU4XjMnciJurtu?S5QOY+yH6;T#lFhZ_BsDO?87
zF)`_C=BH`57#ETeRO=8Ql}<mnNe|o6olY{xRat=P1GFHiOUexGH293-8)0w5@<$Yt
zn7x98!Svg-2_6lpFn-D^r;EMy)Z1PfbJ;W{x1Tw0XAt*T#4;b?r?PkN2X6|=nB&8{
zKi*ty(YN)i@HlA(wY_t~lWI23o!6vqEh(tEJbzsuX({|L=wNl?p3x5aJ|eA7sGq}}
z1tX_s*Q8}%U$yb-i41d6k|`%?4E8BsMI4FC24*O^4v{q<E@MIZKF`eYCax7f)iobS
z?t=1~SM7ue+cSIGln7PT!(GaK{iBuW1smsWG~kF^(N(G;k}`H54vX%b4;Sy+scIt0
z5H1LHych!}MK^@@4WTf7rqXwY#P5Y*R0HIHu+AkAYUR?D>KUoln8%9N4D6}`I@nWR
zADYZW)x7eSs~t<DDQVf9Av2O{eX~bw2sW06e2CBbI+%2tPDEkEdQR;03*__01#uOY
zs#OI(LrI7Dan;jH%$dQyMN43tL{){$9SJG0CCL4dVOy8>W<gkdJz);_D2;>f(jnwo
z9-=)rnmEg$7ag7DGq3uL<J1SRh3lm0wv6uXJ3fpeT663f6-vy)qARgkuvK)WtYG5~
z9EFFWuD85hPw<_=XqT>V?$w_uhc?w~t4w1jhw=>%V(o!S688QiQU}az?~J4(t<DBW
zw)<9&rm$J4H?n^-zE+}d#g#I~&r5O&n<?MUIPv=3y^~!!%8~LE$#%Vb)>R+%{#^`f
zVeW~h+u*8#1n=aay`UG_TpntHXV_%X+X68Q#9NhZw=ZX%opZp^XQL#S*$;)fn?Ihj
z3Rn4(7GKM);K}(_zAuAoIt1-ox(c*JbL6IGdFyZ(-kiy(&86vZvDALASCyTAgSJX3
zz5I*ahQIULz(6E*B}zrpSUvmEdhX1$Uz)VM%AjPMo(w?#tB0<}_TcB&d`G-T$m<(n
z9j`@vct?CRGF<ozG5Dsfo88*1WmQgV4IW;&i3#m~#vy^}&e;n~gJ>qqfwe6^<REs~
zV<2;OQ|Q)#RFcFo*=xPdi=@5EBvnoD2ai0S&?g&vDU^02DQ1Sat?qP@x5EVeGaDYe
zyeY93G@(mD+ltj{RCl&=@cs;v6;eNH@3N)8G<0)acjT;7#oPXACAuho^?>P|ES9eL
zX10ZT0CdX4Kd*JNh-J*yKf*T|j)^^4Jc&AYJ=`ns2fyXSWSZ37vQL^T<|A(%5qhX?
z&fNZ`$?ctH*inFo^%*i=kG3z~8-&CR8+~HZwG`gPB=wYI1?4kIRxV;^3jgHDN|wt!
znE+bCqKu<RBrZ$bp-@^vW(#rXZ8ePR(1lxs^v;-@+1W4F@hrp~OW(r*;Y^;`88}NL
zEwFIv9))wkBJFob%-7*ln`8nyvp*qbezaUFyX>lsz7!y(PPWnJo!o_C?@g~YD&*~k
z!O4;Cbtw_%I*f?JqGG}gDAWrBnv)L>jgUUVN1~&P-a`;ZVIdxj1uTh^$t?^$;uzle
z8Te?^tf^;6OEA)<qLc-z71=Uv0vMxgOslmnoeK#M!{SpIgexwBp>Lri#_mUX#(2|E
zuX?7Or>1~3_{+%rm^9eUWEn4CR??SV0-=gh5RH5CSxxU{faZW^`28|-Z^$z9C1CO9
zBR;UM1_SI(_M_+Mw}|j?kxRM|6mjd!1et^O^0?+Y;l|lG43VU@VvZvk>9>@P>H#{-
zrbMnZC_$#jF(<p}&N}pwRe2c~Ip(W&w#EsP`9FoU1WNU}FLJ&{6Emib5f|<@Z>iJ7
z=S2#2(qqx8&}maFDEqj$s<Qx-w*m_hn6~qCY0vfhtPR<b>bWAg$IM7O6x|uh4<ewp
zW78$}KRiCLR97a4mB@M-M?$e<=iDix6pN)=l<B}hiHHx@ko-%*k|8>!yq*%@uMc5w
z*ILyylV7gJDyL@Et%O+wI#x?_08bITSrVPvb&hnyG5$I=eg2s1azlpO6+*!`mPrMY
z>dqA7Izwr5J@>Z}Y>_I<2|{7+2+)}!*IHN${R$a}bsWQ}vPI;WR2D4%N^rAyolKzX
zi>s8xWuk5fOI0XK*O>En`S4vVvK~pFH$A}68GLUnbb02q)k6a1NW#U<6j&{P5i~=#
z`D;#DTxp6>m>5eIEW@KUgktSDS(Zd`pXahU_CV^*shM;4v3A_kyyYXAQR>9)d*MwM
zORxKC5J@rG5{}3;{-fqEkMm3_Pm+y?3CT#AbMpsn{a5nT!7Hy1>)$3iddB!A)U;qx
zP##mcS~t<#4kvADM{Q(ODw@s~$wkgDD4kAE=RC{Kb8kYAl6)&u8Og*81LPddD05mL
z*Ids+c8qW#$;E0}l@<JUQam3IJ~#_~jO+H$f2UjzO5V)U)VO~wrV_vAf3~~?-ND;i
zAxBmzH|h2($&g&Inuh)|8rZNq@bP(Asd}<(Sx(EwLYTYwvtSZHVZABj5<woIHqB^t
zD}0x&r3a0ItwAwB&;R_Hmlg>i^RGqv=k2^WG%s$;|K;txgkHL*|Ft#ql7{P_mxTTw
Te19FQ|A+5|jRfk+lLGiJiN-M7
literal 0
HcmV?d00001
diff --git a/Projects/MultiParticleRES/utilities/ThreeBodyCS.cc b/Projects/MultiParticleRES/utilities/ThreeBodyCS.cc
new file mode 100644
index 000000000..5674aae83
--- /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 000000000..9f03d73f9
--- /dev/null
+++ b/Projects/MultiParticleRES/utilities/rootlogon.C
@@ -0,0 +1 @@
+{ int status = gSystem->Load("libgsl"); }
--
GitLab