diff --git a/.gitignore b/.gitignore
index 2c26578..3cdaf75 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,4 +9,6 @@ build/
 *_job.txt
 Config_EIC_*
 sjdkay_job.txt
-*#*
\ No newline at end of file
+*#*
+.cc*
+.h*
\ No newline at end of file
diff --git a/Batch_Submission_EIC.sh b/Batch_Submission_EIC.sh
index 4f54b3b..bec4b33 100755
--- a/Batch_Submission_EIC.sh
+++ b/Batch_Submission_EIC.sh
@@ -14,7 +14,7 @@ fi
 if [[ "$#" -ne 7 && "$#" -ne 8 ]]; then
     echo ""
     echo "!!! ERROR !!! - Expected 7 or 8 arguments - !!! ERROR !!!"
-    echo "Expect - NumFiles NumEvents EBeamE HBeamE OutputType InteractionPoint Particle Hadron(optional)"
+    echo "Expect - NumFiles NumEvents EBeamE HBeamE OutputType InteractionPoint Ejectile RecoilHadron(optional)"
     echo "See the Config_EIC.json file or the README for options and try again, exiting"
     echo "!!! ERROR !!! - Expected 7 or 8 arguments - !!! ERROR !!!"
     echo ""
@@ -28,35 +28,35 @@ EBeamE=$3
 HBeamE=$4
 OutputType=$5
 InteractionPoint=$6
-Particle=$7
+Ejectile=$7
 
 # If K+ specified, check the 8th argument, expect this to exist for K+, if it does NOT (case 1), set a default
-if [[ $Particle == "K+" && -z "$8" ]]; then
+if [[ $Ejectile == "K+" && -z "$8" ]]; then
     echo "!!! WARNING !!! - For K+ production expect a hadron specified, defaulting to Lambda - !!! WARNING !!!"
-    Hadron="Lambda"
-elif [[ $Particle == "K+" && ! -z "$8" ]]; then # If 8th argument is not a blank string (i.e. it exists), set the Hadron to this
-    Hadron=$8
-else # Any other case (non K+), set Hadron to be a blank string. We don't actually care for Pi+, Pi0 production etc.
-    Hadron=""
+    RecoilHadron="Lambda"
+elif [[ $Ejectile == "K+" && ! -z "$8" ]]; then # If 8th argument is not a blank string (i.e. it exists), set the RecoilHadron to this
+    RecoilHadron=$8
+else # Any other case (non K+), set RecoilHadron to be a blank string. We don't actually care for Pi+, Pi0 production etc.
+    RecoilHadron=""
 fi
 
 i=1
 while [[ $i -le $NumFiles ]]; do
     # This is the name of the job submission script the shell script creates
-    batch="${USER}_EICDempGen_${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}_Job.txt" # The name of the job submission script it'll create each time
+    batch="${USER}_EICDempGen_${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}_Job.txt" # The name of the job submission script it'll create each time
     echo "Running ${batch} for file ${i}"
     cp /dev/null ${batch}
     RandomSeed=$(od -An -N3 -i /dev/urandom)
     echo "#!/bin/csh" >> ${batch} # Tells your job which shell to run in
-    echo "#PBS -N DEMPGen_${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}" >> ${batch} # Name your job                     
+    echo "#PBS -N DEMPGen_${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}" >> ${batch} # Name your job                     
     echo "#PBS -m abe" >> ${batch} # Email you on job start, end or error
     #echo "#PBS -M ${USER}@jlab.org" >>${batch} # Your email address, change it to be what you like
     echo "#PBS -r n" >> ${batch} # Don't re-run if it crashes
-    echo "#PBS -o  /home/${USER}/trq_output/${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}.out" >> ${batch} # Output directory and file name, set to what you like
-    echo "#PBS -e  /home/${USER}/trq_output/${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}.err" >> ${batch} # Error output directory and file name
+    echo "#PBS -o  /home/${USER}/trq_output/${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}.out" >> ${batch} # Output directory and file name, set to what you like
+    echo "#PBS -e  /home/${USER}/trq_output/${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}.err" >> ${batch} # Error output directory and file name
     echo "date" >> ${batch} 
     echo "cd /home/apps/DEMPgen/" >> ${batch} # Tell your job to go to the directory with the script you want to run
-    echo "./Process_EIC.csh ${i} ${NumEvents} ${EBeamE} ${HBeamE} ${RandomSeed} ${OutputType} ${InteractionPoint} ${Particle} ${Hadron}" >> ${batch} # Run your script, change this to what you like
+    echo "./Process_EIC.csh ${i} ${NumEvents} ${EBeamE} ${HBeamE} ${RandomSeed} ${OutputType} ${InteractionPoint} ${Ejectile} ${RecoilHadron}" >> ${batch} # Run your script, change this to what you like
     echo "date">>${batch}
     echo "exit">>${batch} # End of your job script
     echo "Submitting batch"
diff --git a/Config_EIC.json b/Config_EIC.json
index 3777c99..6e7f2ea 100644
--- a/Config_EIC.json
+++ b/Config_EIC.json
@@ -23,13 +23,14 @@
     //**************************************
     /// This section if for EIC simulation only
     "Kinematics_type" : 1,         // Kinematics type (1->FF, 2->TSSA)
-    "particle": "Pion+", 		   // Choices: omega, pi+, pi0, K+
-    "hadron": "neutron",                // Choices: Neutron, proton, Lambda or Sigma0
+    "ejectile": "Pion+", 		   // Choices: omega, pi+, pi0, K+
+    "recoil_hadron": "neutron",                // Choices: Neutron, proton, Lambda or Sigma0
     "ebeam": 5,                         // Electron beam energy in GeV
     "hbeam": 100,                        // Hadron beam energy in GeV
     "hbeam_part":"Proton",           // Hadron beam particle, proton, deut or helium3, work in progress, will need to be added to shell script later on as an argument
     "OutputType": "HEPMC3",		   // choices: LUND (Docker), Pythia6 (ECCE Fun4All) or HEPMC3 (ePIC)   
     "det_location": "ip6",		   // choices: ip6 for STAR, ip8 for PHENIX
+    "calc_method": "Analytical",            // choices: Analytical or Solve, affects how the ejectile/recoil hadron 4-vectors are calculated. Default is analytical
     "Ee_Low": 0.5, // The minimum scattered electron energy that will be generated as a fraction of the electron beam energy
     "Ee_High": 2.5, // The maximum scattered electron energy that will be generated as a fraction of the electron beam energy
     "e_Theta_Low": 60.0, // The minimum scattered electron angle (theta) that will be generated in degrees
diff --git a/JLab_Batch_Submission.sh b/JLab_Batch_Submission.sh
index fa1dad0..2c8158c 100755
--- a/JLab_Batch_Submission.sh
+++ b/JLab_Batch_Submission.sh
@@ -11,7 +11,7 @@ echo "Running as ${USER}" # Checks who you're running this as
 if [[ "$#" -ne 7 && "$#" -ne 8 ]]; then
     echo ""
     echo "!!! ERROR !!! - Expected 7 or 8 arguments - !!! ERROR !!!"
-    echo "Expect - NumFiles NumEvents EBeamE HBeamE OutputType InteractionPoint Particle Hadron(optional)"
+    echo "Expect - NumFiles NumEvents EBeamE HBeamE OutputType InteractionPoint Ejectile RecoilHadron(optional)"
     echo "See the Config_EIC.json file or the README for options and try again, exiting"
     echo "!!! ERROR !!! - Expected 7 or 8 arguments - !!! ERROR !!!"
     echo ""
@@ -25,16 +25,16 @@ EBeamE=$3
 HBeamE=$4
 OutputType=$5
 InteractionPoint=$6
-Particle=$7
+Ejectile=$7
 
 # If K+ specified, check the 8th argument, expect this to exist for K+, if it does NOT (case 1), set a default
-if [[ $Particle == "K+" && -z "$8" ]]; then
+if [[ $Ejectile == "K+" && -z "$8" ]]; then
     echo "!!! WARNING !!! - For K+ production expect a hadron specified, defaulting to Lambda - !!! WARNING !!!"
-    Hadron="Lambda"
-elif [[ $Particle == "K+" && ! -z "$8" ]]; then # If 8th argument is not a blank string (i.e. it exists), set the Hadron to this
-    Hadron=$8
-else # Any other case (non K+), set Hadron to be a blank string. We don't actually care for Pi+, Pi0 production etc.
-    Hadron=""
+    RecoilHadron="Lambda"
+elif [[ $Ejectile == "K+" && ! -z "$8" ]]; then # If 8th argument is not a blank string (i.e. it exists), set the RecoilHadron to this
+    RecoilHadron=$8
+else # Any other case (non K+), set RecoilHadron to be a blank string. We don't actually care for Pi+, Pi0 production etc.
+    RecoilHadron=""
 fi
 
 Workflow="EIC_DEMPGen_${USER}" # Change this as desired
@@ -47,16 +47,16 @@ while true; do
             (
 		while [[ $i -le $NumFiles ]]; do
 		    # This is the name of the job submission script the shell script creates
-		    batch="${USER}_EICDempGen_${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}_Job.txt" # The name of the job submission script it'll create each time
+		    batch="${USER}_EICDempGen_${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}_Job.txt" # The name of the job submission script it'll create each time
 		    echo "Running ${batch} for file ${i}"
 		    cp /dev/null ${batch}
 		    RandomSeed=$(od -An -N3 -i /dev/urandom)
 		    echo "PROJECT: c-kaonlt"  >> ${batch} # Is eic a valid project?
 		    echo "TRACK: analysis" >> ${batch}
-		    echo "JOBNAME: DEMPGen_${EBeamE}on${HBeamE}_${Particle}${Hadron}_${InteractionPoint}_${NumEvents}_${i}" >> ${batch}
+		    echo "JOBNAME: DEMPGen_${EBeamE}on${HBeamE}_${Ejectile}${RecoilHadron}_${InteractionPoint}_${NumEvents}_${i}" >> ${batch}
                     echo "MEMORY: 2000 MB" >> ${batch} # Request 2GB RAM - probably too much
 		    echo "CPU: 1" >> ${batch} # Request 1 CPU core per job
-		    echo "COMMAND:/group/eic/users/${USER}/DEMPGen/Process_EIC_iFarm.csh ${i} ${NumEvents} ${EBeamE} ${HBeamE} ${RandomSeed} ${OutputType} ${InteractionPoint} ${Particle} ${Hadron}" >> ${batch}
+		    echo "COMMAND:/group/eic/users/${USER}/DEMPGen/Process_EIC_iFarm.csh ${i} ${NumEvents} ${EBeamE} ${HBeamE} ${RandomSeed} ${OutputType} ${InteractionPoint} ${Ejectile} ${RecoilHadron}" >> ${batch}
                     echo "MAIL: ${USER}@jlab.org" >> ${batch}
 		    echo "Submitting batch"
 		    eval "swif2 add-jsub ${Workflow} -script ${batch} 2>/dev/null" # Swif2 job submission, uses old jsub scripts
diff --git a/Process_EIC.csh b/Process_EIC.csh
index 647028f..b0f03e2 100755
--- a/Process_EIC.csh
+++ b/Process_EIC.csh
@@ -4,7 +4,7 @@
 
 echo""
 echo "This file is intended to be run as part of a batch job submission, however, you can also run it on its own."
-echo "Expected input is - FileNumber NumberOfEvents ElectronBeamEnergy HadronBeamEnergy RandomSeed OutputType IP Particle Hadron(Optional, for K+)"
+echo "Expected input is - FileNumber NumberOfEvents ElectronBeamEnergy HadronBeamEnergy RandomSeed OutputType IP Ejectile RecoilHadron(Optional, for K+)"
 echo "Please see the README for more info."
 echo ""
 
@@ -22,35 +22,35 @@ set HBeamE=$4
 set RandomSeed=$5
 set OutputType=$6
 set InteractionPoint=$7
-set Particle=$8
+set Ejectile=$8
 
-if ($Particle == "K+" && $#argv == 8 ) then
+if ($Ejectile == "K+" && $#argv == 8 ) then
     echo "! WARNING !"
     echo "! WARNING ! - For K+ production expect a hadron specified, defaulting to Lambda - ! WARNING !"
     echo "! WARNING !"
-    set Hadron="Lambda"
-else if ($Particle == "K+" && $#argv == 9 ) then
-    set Hadron=$9
+    set RecoilHadron="Lambda"
+else if ($Ejectile == "K+" && $#argv == 9 ) then
+    set RecoilHadron=$9
 else 
-    set Hadron=""
+    set RecoilHadron=""
 endif
 
-echo "Running target polarisation up, FF setting for file $FileNum with $NumEvents events per file for $EBeamE GeV e- on $HBeamE GeV p using random seed $RandomSeed, using $OutputType format output for $Particle $Hadron events."
+echo "Running target polarisation up, FF setting for file $FileNum with $NumEvents events per file for $EBeamE GeV e- on $HBeamE GeV p using random seed $RandomSeed, using $OutputType format output for $Ejectile $RecoilHadron events."
     
 # Set the config file name based upon inputs
-set ConfigFilename = 'Config_EIC_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.json'
+set ConfigFilename = 'Config_EIC_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.json'
 
 # Copy the default config file to our constructed filename
 cp Config_EIC.json $ConfigFilename
 
 # Use sed commands to change our config file based upon inputs
-sed -i 's/"file_name" \:.*/"file_name" \: "DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'",/' $ConfigFilename
+sed -i 's/"file_name" \:.*/"file_name" \: "DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'",/' $ConfigFilename
 sed -i 's/"n_events" \:.*/"n_events" \: '$NumEvents',/' $ConfigFilename
 sed -i 's/"generator_seed"\:.*/"generator_seed" \: '$RandomSeed',/' $ConfigFilename
 sed -i 's/"ebeam"\:.*/"ebeam" \: '$EBeamE',/' $ConfigFilename
 sed -i 's/"hbeam"\:.*/"hbeam" \: '$HBeamE',/' $ConfigFilename
-sed -i 's/"particle"\:.*/"particle" \: "'$Particle'",/' $ConfigFilename
-sed -i 's/"hadron"\:.*/"hadron" \: "'$Hadron'",/' $ConfigFilename
+sed -i 's/"ejectile"\:.*/"ejectile" \: "'$Ejectile'",/' $ConfigFilename
+sed -i 's/"recoil_hadron"\:.*/"recoil_hadron" \: "'$RecoilHadron'",/' $ConfigFilename
 sed -i 's/"det_location"\:.*/"det_location" \: "'$InteractionPoint'",/' $ConfigFilename
 sed -i 's/"OutputType"\:.*/"OutputType"\: "'$OutputType'",/' $ConfigFilename
 
@@ -60,8 +60,8 @@ cd data/
 sleep 5
 
 # Filename as it's created is a bit odd, so rename it
-set OriginalOutput = 'eic_input_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.dat'
-set RenamedOutput = 'eic_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.dat'
+set OriginalOutput = 'eic_input_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.dat'
+set RenamedOutput = 'eic_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.dat'
 mv "OutputFiles/"$OriginalOutput "OutputFiles/"$RenamedOutput
 
 rm -rf ../$ConfigFilename
diff --git a/Process_EIC_iFarm.csh b/Process_EIC_iFarm.csh
index 35cb6c3..bcea831 100755
--- a/Process_EIC_iFarm.csh
+++ b/Process_EIC_iFarm.csh
@@ -22,7 +22,7 @@ cd "$DEMPGENPath"
 
 echo""
 echo "This file is intended to be run as part of a batch job submission, however, you can also run it on its own."
-echo "Expected input is - FileNumber NumberOfEvents ElectronBeamEnergy HadronBeamEnergy RandomSeed OutputType IP Particle Hadron(Optional, for K+)"
+echo "Expected input is - FileNumber NumberOfEvents ElectronBeamEnergy HadronBeamEnergy RandomSeed OutputType IP Ejectile Recoil_Hadron(Optional, for K+)"
 echo "Please see the README for more info."
 echo ""
 
@@ -40,23 +40,23 @@ set HBeamE=$4
 set RandomSeed=$5
 set OutputType=$6
 set InteractionPoint=$7
-set Particle=$8
+set Ejectile=$8
 
-if ($Particle == "K+" && $#argv == 8 ) then
+if ($Ejectile == "K+" && $#argv == 8 ) then
     echo "! WARNING !"
-    echo "! WARNING ! - For K+ production expect a hadron specified, defaulting to Lambda - ! WARNING !"
+    echo "! WARNING ! - For K+ production expect a recoil hadron specified, defaulting to Lambda - ! WARNING !"
     echo "! WARNING !"
-    set Hadron="Lambda"
-else if ($Particle == "K+" && $#argv == 9 ) then
-    set Hadron=$9
+    set RecoilHadron="Lambda"
+else if ($Ejectile == "K+" && $#argv == 9 ) then
+    set RecoilHadron=$9
 else 
-    set Hadron=""
+    set RecoilHadron=""
 endif
 
-echo "Running target polarisation up, FF setting for file $FileNum with $NumEvents events per file for $EBeamE GeV e- on $HBeamE GeV p using random seed $RandomSeed, using $OutputType format output for $Particle $Hadron events."
+echo "Running target polarisation up, FF setting for file $FileNum with $NumEvents events per file for $EBeamE GeV e- on $HBeamE GeV p using random seed $RandomSeed, using $OutputType format output for $Ejectile $RecoilHadron events."
     
 # Set the config file name based upon inputs
-set ConfigFilename = $DEMPGENPath'/Config_EIC_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.json'
+set ConfigFilename = $DEMPGENPath'/Config_EIC_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.json'
 
 # Copy the default config file to our constructed filename
 cp "$DEMPGENPath/Config_EIC.json" $ConfigFilename
@@ -67,8 +67,8 @@ sed -i 's/"n_events" \:.*/"n_events" \: '$NumEvents',/' $ConfigFilename
 sed -i 's/"generator_seed"\:.*/"generator_seed" \: '$RandomSeed',/' $ConfigFilename
 sed -i 's/"ebeam"\:.*/"ebeam" \: '$EBeamE',/' $ConfigFilename
 sed -i 's/"hbeam"\:.*/"hbeam" \: '$HBeamE',/' $ConfigFilename
-sed -i 's/"particle"\:.*/"particle" \: "'$Particle'",/' $ConfigFilename
-sed -i 's/"hadron"\:.*/"hadron" \: "'$Hadron'",/' $ConfigFilename
+sed -i 's/"ejectile"\:.*/"ejectile" \: "'$Ejectile'",/' $ConfigFilename
+sed -i 's/"recoil_hadron"\:.*/"recoil_hadron" \: "'$RecoilHadron'",/' $ConfigFilename
 sed -i 's/"det_location"\:.*/"det_location" \: "'$InteractionPoint'",/' $ConfigFilename
 sed -i 's/"OutputType"\:.*/"OutputType"\: "'$OutputType'",/' $ConfigFilename
 
@@ -78,8 +78,8 @@ eval $DEMPGENPath/build/DEMPgen $ConfigFilename
 sleep 5
 
 # Filename as it's created is a bit odd, so rename it
-set OriginalOutput = $DEMPGENPath'/data/OutputFiles/eic_input_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.dat'
-set RenamedOutput =  $DEMPGENPath'/data/OutputFiles/eic_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Particle$Hadron'_'$NumEvents'_'$FileNum'.dat'
+set OriginalOutput = $DEMPGENPath'/data/OutputFiles/eic_input_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.dat'
+set RenamedOutput =  $DEMPGENPath'/data/OutputFiles/eic_DEMPGen_'$EBeamE'on'$HBeamE'_'$InteractionPoint'_'$Ejectile$RecoilHadron'_'$NumEvents'_'$FileNum'.dat'
 
 mv $OriginalOutput $RenamedOutput
 
diff --git a/README.md b/README.md
index cf82515..041c2c0 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,5 @@
 # DEMPgen
-Event generator for Deep Exclusive Meson Production
+Event generator for Deep Exclusive Meson Production.
 
 ## Building
 
@@ -9,7 +9,7 @@ Event generator for Deep Exclusive Meson Production
 
 - After downloading the source create a build directory and cd to it. Take note of the location of the source directory (where CMakeLists.txt should be stored) and run the commands:
 
-  - mdkir build.
+  - mdkir build
   - cd build
   - cmake ..
   - make -j8
@@ -25,7 +25,7 @@ Event generator for Deep Exclusive Meson Production
 
 ### Building on the iFarm
 
-- Building on the JLab iFarm requires you to set up some software versions beforehand, to build successfully, I did the following
+- Building on the JLab iFarm requires you to set up some software versions beforehand, to build successfully, I did the following - 
   - Comment out any CUE or other initialisation in your .login/.cshrc scripts
   - Login to an ifarm node
   - module load cmake/3.19.4
@@ -36,7 +36,22 @@ Event generator for Deep Exclusive Meson Production
 
 ## Configuration
 
-The file Config.json contains all the configuration options. Use this as a template for other configuration files, which may be given as an argument to the event generator.
+The file Config_EIC.json and Config_SoLID.json contain all the configuration options for EIC and SoLID event generation respectively. Use these as a template for other configuration files, which may be given as an argument to the event generator.
+- Note that Config_EIC.json is used in shell scripts (see below) that run large numbers of jobs on batch queueing systems. *Do not* modify this file unless you know what you're doing with it or you really need to edit something.
+  - Copying the template to a new file and editing that is strongly recommended in all cases.
+
+
+### Ejectile calculation methods
+
+- The EIC module of DEMPgen has two different calculation methods that may be used to calculate the ejectile properties.
+- These are referred to as the "Analytical" and "Solve" methods. Either can be used in a simulation run, just change the config file to switch between them.
+- In your config .json file, select between the two versions using the "calc_method" argument, this can be set to Analytical or Solve, e.g.
+  - "calc_method": "Analytical"
+  - OR
+  - "calc_method": "Solve"
+  - If the method is not entered exactly as above, it will not be recognised and the generator will default to the analytical method.
+
+- If you are using the shell scripts provided to run a lot of simulations, modify Config_EIC.json BEFORE running the shell scripts.
 
 ## Output
 
@@ -96,28 +111,29 @@ This project uses [JsonCpp](https://github.com/open-source-parsers/jsoncpp "Json
 
 ## Processing Scripts and json examples
 
-- There are several shell scripts in the main directory that can be used to run the generator and produce some output, information on these scripts is provided below
-- To remove clutter from the main directory, example .json scripts were moved to a new folder
+- There are several shell scripts in the main directory that can be used to run the generator and produce some output, information on these scripts is provided below.
+- To remove clutter from the main directory, example .json scripts were moved to a new folder -
   - json_examples
+- Remember to set the ejectile calculation method in the Config_EIC.json file *before* running these scripts if you do not want to use the default analytical method.
 
 ### Process_EIC.csh
 
 !!! NOTICE !!!  
-This script copies Config_EIC.json and formats a new file based upon this, DO NOT MODIFY Config_EIC.json if you want to use this script!  
+This script copies Config_EIC.json and formats a new file based upon this, DO NOT MODIFY Config_EIC.json (other than the calculation method) if you want to use this script!  
 !!! NOTICE !!!  
 
 - To facilitate the submission of batch jobs, I created a csh script to automatically construct .json config files and run them. This script can also be utilised to run the generator manually, without the need to go and edit a json file. 
 - The script requires 8 arguments (which is a lot, I know), but in the K+ case, it expects 9. They are as follows -  
 
-  - Arg 1 - FileNum -> For batch running, we typically run X files of Y events, this argument is just X, if you're running manually as a test, just input 1 or whatever you fancy  
-  - Arg 2 - NumEvents -> The number of events thrown for this file, set this to whatever you want to run. For reference, with the Pi+/K+ generator, 1B files takes ~1 hour  
-  - Arg 3 - EBeamE -> The electron beam energy, set this to whatever you want, typically, we use 5, 10 or 18 (the nominal max for the EIC)  
-  - Arg 4 - HBeamE -> The hadron beam energy, again, set this to whatevr you want. Typically we use 41, 100 or 275 (41 and 275 being the nominal min/max)  
-  - Arg 5 - RandomSeed -> The random seed, self explanatory. Set this however you like, the batch submission job randomly generates a random seed to feed in here  
-  - Arg 6 - OutputType -> The format of the output file, select from LUND, Pythia6 (for ECCE/Fun4All) or HEPMC3 (for ePIC), the default is HEPMC3 if your choice is invalid  
-  - Arg 7 - InteractionPoint -> The interaction point, choose from ip6 or ip8. The default is ip6 if your choice is invalid  
-  - Arg 8 - Particle -> The produced particle (meson) in the reaction, choose from omega, pi+, pi0 or K+  
-  - Arg 9 - Hadron -> OPTIONAL - This only matters if you select K+ as the particle, in this case, choose from Lambda or Sigma0 here. If your choice is invalid (or you don't specify arg9), the default is Lambda  
+  - Arg 1 - FileNum -> For batch running, we typically run X files of Y events, this argument is just X, if you're running manually as a test, just input 1 or whatever you fancy.
+  - Arg 2 - NumEvents -> The number of events thrown for this file, set this to whatever you want to run. For reference, with the Pi+/K+ generator, 1B files takes ~1 hour.
+  - Arg 3 - EBeamE -> The electron beam energy, set this to whatever you want, typically, we use 5, 10 or 18 (the nominal max for the EIC). 
+  - Arg 4 - HBeamE -> The hadron beam energy, again, set this to whatevr you want. Typically we use 41, 100 or 275 (41 and 275 being the nominal min/max).
+  - Arg 5 - RandomSeed -> The random seed, self explanatory. Set this however you like, the batch submission job randomly generates a random seed to feed in here.
+  - Arg 6 - OutputType -> The format of the output file, select from LUND, Pythia6 (for ECCE/Fun4All) or HEPMC3 (for ePIC), the default is HEPMC3 if your choice is invalid.
+  - Arg 7 - InteractionPoint -> The interaction point, choose from ip6 or ip8. The default is ip6 if your choice is invalid.
+  - Arg 8 - Ejectile -> The produced ejectile (meson) in the reaction, choose from omega, pi+, pi0 or K+.
+  - Arg 9 - RecoilHadron -> OPTIONAL - This only matters if you select K+ as the ejectile, in this case, choose from Lambda or Sigma0 here. If your choice is invalid (or you don't specify arg9), the default is Lambda.
   
 - So as an example if you executed the following -  
   
@@ -130,33 +146,33 @@ This script copies Config_EIC.json and formats a new file based upon this, DO NO
 - This script creates and submits batch jobs. It is designed for use with the torque queueing system on Lark at the University of Regina.
 The jobs the script creates and submits all execute the Process_EIC.csh script described above. This script requries a very similar set of arguments -  
 
-  - Arg 1 - NumFiles -> The batch script is designed to run X jobs of Y events, this number is just X, the number of files you want to run  
-  - Arg 2 - NumEvents -> The number of events thrown for this file, set this to whatever you want to run. For reference, with the Pi+/K+ generator, 1B files takes ~1 hour  
-  - Arg 3 - EBeamE -> The electron beam energy, set this to whatever you want, typically, we use 5, 10 or 18 (the nominal max for the EIC)  
-  - Arg 4 - HBeamE -> The hadron beam energy, again, set this to whatevr you want. Typically we use 41, 100 or 275 (41 and 275 being the nominal min/max)  
-  - Arg 5 - OutputType -> The format of the output file, select from LUND, Pythia6 (for ECCE/Fun4All) or HEPMC3 (for ePIC), the default is HEPMC3 if your choice is invalid  
-  - Arg 6 - InteractionPoint -> The interaction point, choose from ip6 or ip8. The default is ip6 if your choice is invalid  
-  - Arg 7 - Particle -> The produced particle (meson) in the reaction, choose from omega, pi+, pi0 or K+  
-  - Arg 8 - Hadron -> OPTIONAL - This only matters if you select K+ as the particle, in this case, choose from Lambda or Sigma0 here. If your choice is invalid (or you don't specify arg9), the default is Lambda  
+  - Arg 1 - NumFiles -> The batch script is designed to run X jobs of Y events, this number is just X, the number of files you want to run.
+  - Arg 2 - NumEvents -> The number of events thrown for this file, set this to whatever you want to run. For reference, with the Pi+/K+ generator, 1B files takes ~1 hour.
+  - Arg 3 - EBeamE -> The electron beam energy, set this to whatever you want, typically, we use 5, 10 or 18 (the nominal max for the EIC).
+  - Arg 4 - HBeamE -> The hadron beam energy, again, set this to whatevr you want. Typically we use 41, 100 or 275 (41 and 275 being the nominal min/max).
+  - Arg 5 - OutputType -> The format of the output file, select from LUND, Pythia6 (for ECCE/Fun4All) or HEPMC3 (for ePIC), the default is HEPMC3 if your choice is invalid. 
+  - Arg 6 - InteractionPoint -> The interaction point, choose from ip6 or ip8. The default is ip6 if your choice is invalif.
+  - Arg 7 - Ejectile -> The produced ejectile (meson) in the reaction, choose from omega, pi+, pi0 or K+.
+  - Arg 8 - RecoilHadron -> OPTIONAL - This only matters if you select K+ as the ejectile, in this case, choose from Lambda or Sigma0 here. If your choice is invalid (or you don't specify arg9), the default is Lambda.
 
-- The script automatically generates a random seed itself using the /dev/urandom function  
+- The script automatically generates a random seed itself using the /dev/urandom function.
 
 ### Process_EIC_iFarm.csh
 
-- This version is for use on the JLab iFarm/Farm
-- This script uses the same arguments as Process_EIC.csh
+- This version is for use on the JLab iFarm/Farm.
+- This script uses the same arguments as Process_EIC.csh.
 - If you are processing interactively (i.e. on the iFarm), you will need to make sure that you have executed -
   - module use /group/halla/modulefiles
   - module load root
-- You should also check the path set at the top looks OK
-    - By default, /eic/users/${USER}/DEMPGen is assumed
+- You should also check the path set at the top looks OK.
+    - By default, /eic/users/${USER}/DEMPGen is assumed.
 
 ### JLab_Batch_Submission.sh
 
-- This version should be used to submit jobs to the Farm batch queueing system (swif2)
-- It uses the same arguments as Batch_Submission_EIC.sh
-- You should also check the paths set throughout the script look ok (for example, in the COMMAND: ... line)
-    - By default, /eic/users/${USER}/DEMPGen is assumed
+- This version should be used to submit jobs to the Farm batch queueing system (swif2).
+- It uses the same arguments as Batch_Submission_EIC.sh.
+- You should also check the paths set throughout the script look ok (for example, in the COMMAND: ... line).
+    - By default, /eic/users/${USER}/DEMPGen is assumed.
 
 ### json_examples
 
diff --git a/src/eic_evgen/eic.cc b/src/eic_evgen/eic.cc
index 0d42170..f487d23 100644
--- a/src/eic_evgen/eic.cc
+++ b/src/eic_evgen/eic.cc
@@ -33,81 +33,81 @@ using namespace std;
 
 void eic() {
 
-    Int_t target_direction, kinematics_type;
-    Double_t EBeam, HBeam;
+  Int_t target_direction, kinematics_type;
+  Double_t EBeam, HBeam;
  
-   	cout << "Target Direction (1->Up, 2->Down): "; cin >> target_direction; cout << endl;
-   	cout << "Kinematics type (1->FF, 2->TSSA): ";  cin >> kinematics_type;  cout << endl;
-   	cout << "Enter the number of events: ";        cin >> fNEvents;         cout << endl;
-   	cout << "Enter the file number: ";             cin >> fNFile;           cout << endl;
-	cout << "Enter the electron beam energy: ";    cin >> EBeam;            cout << endl;
-	cout << "Enter the hadron beam energy: ";      cin >> HBeam;            cout << endl;
+  cout << "Target Direction (1->Up, 2->Down): "; cin >> target_direction; cout << endl;
+  cout << "Kinematics type (1->FF, 2->TSSA): ";  cin >> kinematics_type;  cout << endl;
+  cout << "Enter the number of events: ";        cin >> fNEvents;         cout << endl;
+  cout << "Enter the file number: ";             cin >> fNFile;           cout << endl;
+  cout << "Enter the electron beam energy: ";    cin >> EBeam;            cout << endl;
+  cout << "Enter the hadron beam energy: ";      cin >> HBeam;            cout << endl;
  
-//	eic(target_direction, kinematics_type, fNEvents);
+  //	eic(target_direction, kinematics_type, fNEvents);
 
 }
 
 /*--------------------------------------------------*/
 // 18/01/23 - SJDK- This function is never used since eic() is only called with a json object as the argument. Commented out for now, delete later?
-/* void eic(int event_number, int target_direction, int kinematics_type, TString file_name, int fEIC_seed, TString particle, TString hadron, TString det_location, TString OutputType, double EBeam, double HBeam) {
+/* void eic(int event_number, int target_direction, int kinematics_type, TString file_name, int fEIC_seed, TString Ejectile, TString RecoilHadron, TString det_location, TString OutputType, double EBeam, double HBeam) {
 
-   	TString targetname;
-	TString charge;
+   TString targetname;
+   TString charge;
 
-   	if( target_direction == 1 ) targetname = "up";
-  	if( target_direction == 2 ) targetname = "down";
+   if( target_direction == 1 ) targetname = "up";
+   if( target_direction == 2 ) targetname = "down";
 	
-	gKinematics_type = kinematics_type;
-	gfile_name = file_name;
-
-	fNFile = 1;
-
-	fNEvents = event_number;
-
-	fSeed = fEIC_seed;
-	cout << EBeam << " elec " << HBeam << " hadrons" << endl; 
-	fEBeam = EBeam;
-	fPBeam = HBeam;
-
-	pim* myPim = new pim(fSeed);
-  	myPim->Initilize();
-	// 09/02/22 - SJDK - Special case for the kaon, if hadron not specified, default to Lambda
-	if (particle == "K+"){
-	  if (hadron != "Lambda" && hadron != "Sigma0"){
-	    hadron = "Lambda";
-	  }
-	  else{
-	    hadron = ExtractParticle(hadron);
-	  }
-	  Reaction* r1 = new Reaction(particle, hadron);
-	  r1->process_reaction();
-	  delete r1;
-	}
-	else if (particle == "pi+" || particle == "Pion+" ||  particle == "Pi+"){
-	  hadron = "Neutron";
-	  particle = ExtractParticle(particle);
-	  charge = ExtractCharge(particle);
-	  Reaction* r1 = new Reaction(particle, hadron);
-	  r1->process_reaction();
-	  delete r1;
-	}
-	else if (particle == "pi0" || particle == "Pion0" || particle == "Pi0"){
-	  hadron = "Proton";
-	  particle = ExtractParticle(particle);
-	  charge = ExtractCharge(particle);
-	  //Reaction* r1 = new Reaction(particle);
-	  Reaction* r1 = new Reaction(particle, hadron);
-	  r1->process_reaction();
-	  delete r1;
-	}
-	else{
-	  particle = ExtractParticle(particle);
-	  charge = ExtractCharge(particle);
-	  Reaction* r1 = new Reaction(particle);
-	  r1->process_reaction();
-	  delete r1;
-	}
-}
+   gKinematics_type = kinematics_type;
+   gfile_name = file_name;
+
+   fNFile = 1;
+
+   fNEvents = event_number;
+
+   fSeed = fEIC_seed;
+   cout << EBeam << " elec " << HBeam << " RecoilHadrons" << endl; 
+   fEBeam = EBeam;
+   fPBeam = HBeam;
+
+   pim* myPim = new pim(fSeed);
+   myPim->Initilize();
+   // 09/02/22 - SJDK - Special case for the kaon, if RecoilHadron not specified, default to Lambda
+   if (Ejectile == "K+"){
+   if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){
+   RecoilHadron = "Lambda";
+   }
+   else{
+   RecoilHadron = ExtractParticle(RecoilHadron);
+   }
+   Reaction* r1 = new Reaction(Ejectile, RecoilHadron);
+   r1->process_reaction();
+   delete r1;
+   }
+   else if (Ejectile == "pi+" || Ejectile == "Pion+" ||  Ejectile == "Pi+"){
+   RecoilHadron = "Neutron";
+   Ejectile = ExtractParticle(particle);
+   charge = ExtractCharge(Ejectile);
+   Reaction* r1 = new Reaction(Ejectile, RecoilHadron);
+   r1->process_reaction();
+   delete r1;
+   }
+   else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){
+   RecoilHadron = "Proton";
+   Ejectile = ExtractParticle(Ejectile);
+   charge = ExtractCharge(Ejectile);
+   //Reaction* r1 = new Reaction(Ejectile);
+   Reaction* r1 = new Reaction(Ejectile, RecoilHadron);
+   r1->process_reaction();
+   delete r1;
+   }
+   else{
+   Ejectile = ExtractParticle(Ejectile);
+   charge = ExtractCharge(Ejectile);
+   Reaction* r1 = new Reaction(Ejectile);
+   r1->process_reaction();
+   delete r1;
+   }
+   }
 */
 
 /*--------------------------------------------------*/
@@ -115,223 +115,242 @@ void eic() {
 // SJDK 21/12/22 - Note that this is the one that actually gets used, reads in the .json file
 void eic(Json::Value obj) {
    	
-        TString targetname;  
- 	TString charge;
+  TString targetname;  
+  TString charge;
 
-	int target_direction = obj["Targ_dir"].asInt();
- 	gKinematics_type     = obj["Kinematics_type"].asInt();
+  int target_direction = obj["Targ_dir"].asInt();
+  gKinematics_type     = obj["Kinematics_type"].asInt();
 
-   	if( target_direction == 1 ) targetname = "up";
-   	if( target_direction == 2 ) targetname = "down";
+  if( target_direction == 1 ) targetname = "up";
+  if( target_direction == 2 ) targetname = "down";
  
- 	gfile_name = obj["file_name"].asString();
+  gfile_name = obj["file_name"].asString();
  
- 	gPi0_decay = obj["pi0_decay"].asBool();
+  gPi0_decay = obj["pi0_decay"].asBool();
 
- 	fNFile = 1;
- 	fNEvents = obj["n_events"].asUInt64();
+  fNFile = 1;
+  fNEvents = obj["n_events"].asUInt64();
 
- 	fSeed = obj["generator_seed"].asInt();
+  fSeed = obj["generator_seed"].asInt();
 
- 	pim* myPim = new pim(fSeed);
-   	myPim->Initilize();
+  pim* myPim = new pim(fSeed);
+  myPim->Initilize();
  
-//  	TDatime dsTime;
-//  	cout << "Start Time:   " << dsTime.GetHour() << ":" << dsTime.GetMinute() << endl;
-	// 21/12/22 - SJDK - Should do a check if these are defined or not, should crash if not defined or set defaults, see other quantities below
-	TString particle = obj["particle"].asString();
-	TString hadron = obj["hadron"].asString(); // 09/02/22 - SJDK - Added in hadron type argument for K+
-	// SJDK - 08/02/22 - This is terrible, need to change this, particle should just be K+
-	// Add a new flag which, hadron - where this is specified too, then add conditionals elsewhere based on this
-	// New conditional, special case for Kaon	
-	particle = ExtractParticle(particle);
-	charge = ExtractCharge(particle);
-	if(hadron == "Sigma" || hadron == "sigma"){ // SJDK - 31/01/23 - If hadron specified as Sigma, interpret this as Sigma0. Also correct for lower case
-	  hadron = "Sigma0";
-	}
-	if (hadron == "lambda"){ // SJDK - 31/01/23 - Make Lambda selection case insensitive
-	  hadron = "Lambda"; 
-	}
-	if (particle == "K+"){
-	  if (hadron != "Lambda" && hadron != "Sigma0"){
-	    hadron = "Lambda";
-	    cout << "! WARNING !" << endl;
-	    cout << "! WARNING !- K+ production specified but hadron not recognised, deaulting to Lambda - ! WARNING!" << endl;
-	    cout << "! WARNING !" << endl;
-	  }
-	  else{
-	    hadron = ExtractParticle(hadron);
-	  }
-	}
-	// SJDK - 19/12/22 - Specify hadron to neutron/proton for pi+/pi0 production, for pi0 production, may want to adjust? 
-	else if (particle == "pi+" || particle == "Pion+" || particle == "Pi+"){
-	  hadron = "Neutron";
-	}
-	else if (particle == "pi0" || particle == "Pion0" || particle == "Pi0"){
-	  hadron = "Proton";
-	}
-	else { // SJDK -09/02/22 - Note that in future this could be changed to get different hadrons in other reactions if desired
-	  hadron = "";
-	}
+  //  	TDatime dsTime;
+  //  	cout << "Start Time:   " << dsTime.GetHour() << ":" << dsTime.GetMinute() << endl;
+  // 21/12/22 - SJDK - Should do a check if these are defined or not, should crash if not defined or set defaults, see other quantities below
+  TString Ejectile = obj["ejectile"].asString();
+  TString RecoilHadron = obj["recoil_hadron"].asString(); // 09/02/22 - SJDK - Added in RecoilHadron type argument for K+
+  // SJDK - 08/02/22 - This is terrible, need to change this, Ejectile should just be K+
+  // Add a new flag which, RecoilHadron - where this is specified too, then add conditionals elsewhere based on this
+  // New conditional, special case for Kaon	
+  Ejectile = ExtractParticle(Ejectile);
+  charge = ExtractCharge(Ejectile);
+  if(RecoilHadron == "Sigma" || RecoilHadron == "sigma"){ // SJDK - 31/01/23 - If RecoilHadron specified as Sigma, interpret this as Sigma0. Also correct for lower case
+    RecoilHadron = "Sigma0";
+  }
+  if (RecoilHadron == "lambda"){ // SJDK - 31/01/23 - Make Lambda selection case insensitive
+    RecoilHadron = "Lambda"; 
+  }
+  if (Ejectile == "K+"){
+    if (RecoilHadron != "Lambda" && RecoilHadron != "Sigma0"){
+      RecoilHadron = "Lambda";
+      cout << "! WARNING !" << endl;
+      cout << "! WARNING !- K+ production specified but RecoilHadron not recognised, deaulting to Lambda - ! WARNING!" << endl;
+      cout << "! WARNING !" << endl;
+    }
+    else{
+      RecoilHadron = ExtractParticle(RecoilHadron);
+    }
+  }
+  // SJDK - 19/12/22 - Specify RecoilHadron to neutron/proton for pi+/pi0 production, for pi0 production, may want to adjust? 
+  else if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){
+    RecoilHadron = "Neutron";
+  }
+  else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){
+    RecoilHadron = "Proton";
+  }
+  else { // SJDK -09/02/22 - Note that in future this could be changed to get different RecoilHadrons in other reactions if desired
+    RecoilHadron = "";
+  }
 
-	// SJDK 03/04/23 - Change to how Qsq range is set/chosen, could add as an override variable later too
-	// Set min/max Qsq values depending upon particle type
-	if (particle == "pi+" || particle == "Pion+" || particle == "Pi+"){
-	  fQsq_Min = 3.0; fQsq_Max = 35.0;
-	  fW_Min = 2.0; fW_Max = 10.2;
-	}
-	else if (particle == "pi0" || particle == "Pion0" || particle == "Pi0"){
-	  fQsq_Min = 5.0; fQsq_Max = 1000.0;
-	  fW_Min = 2.0; fW_Max = 10.0;
-	}
-	else if (particle == "K+"){
-	  fQsq_Min = 1.0; fQsq_Max = 35.0;
-	  fW_Min = 2.0; fW_Max = 10.0;
-	}
-	else{
-	  fQsq_Min = 5.0; fQsq_Max = 35.0;
-	  fW_Min = 2.0; fW_Max = 10.0;
-	}
+  // SJDK 03/04/23 - Change to how Qsq range is set/chosen, could add as an override variable later too
+  // Set min/max Qsq values depending upon Ejectile type
+  if (Ejectile == "pi+" || Ejectile == "Pion+" || Ejectile == "Pi+"){
+    fQsq_Min = 3.0; fQsq_Max = 35.0;
+    fW_Min = 2.0; fW_Max = 10.2;
+    fT_Max = 1.3;
+  }
+  else if (Ejectile == "pi0" || Ejectile == "Pion0" || Ejectile == "Pi0"){
+    fQsq_Min = 5.0; fQsq_Max = 1000.0;
+    fW_Min = 2.0; fW_Max = 10.0;
+    fT_Max = 0.5;
+  }
+  else if (Ejectile == "K+"){
+    fQsq_Min = 1.0; fQsq_Max = 35.0;
+    fW_Min = 2.0; fW_Max = 10.0;
+    fT_Max = 2.0;
+  }
+  else{
+    fQsq_Min = 5.0; fQsq_Max = 35.0;
+    fW_Min = 2.0; fW_Max = 10.0;
+    fT_Max = 1.3;
+  }
 
-	// SJDK - 01/06/21
-	// Set beam energies from .json read in
-	if (obj.isMember("ebeam")){
-	  fEBeam = obj["ebeam"].asDouble();
-	}
-	else{
-	  fEBeam = 5;
-	  cout << "Electron beam energy not specified in .json file, defaulting to 5 GeV." << endl;
-	}
-	if (obj.isMember("hbeam")){
-	  fPBeam = obj["hbeam"].asDouble();
-	}
-	else{
-	  fPBeam = 100;
-	  cout << "Ion beam energy not specified in .json file, defaulting to 100 GeV." << endl;
-	}
+  // SJDK - 01/06/21
+  // Set beam energies from .json read in
+  if (obj.isMember("ebeam")){
+    fEBeam = obj["ebeam"].asDouble();
+  }
+  else{
+    fEBeam = 5;
+    cout << "Electron beam energy not specified in .json file, defaulting to 5 GeV." << endl;
+  }
+  if (obj.isMember("hbeam")){
+    fPBeam = obj["hbeam"].asDouble();
+    fHBeam = obj["hbeam"].asDouble();
+  }
+  else{
+    fPBeam = 100;
+    fHBeam = 100;
+    cout << "Ion beam energy not specified in .json file, defaulting to 100 GeV." << endl;
+  }
 
-	if (obj.isMember("hbeam_part")){
-	  gBeamPart = obj["hbeam_part"].asString();
-	} 
+  if (obj.isMember("hbeam_part")){
+    gBeamPart = obj["hbeam_part"].asString();
+  } 
+  else {
+    gBeamPart = "Proton";
+  }
+
+  // SJDK - 12/01/22
+  // Set output type as a .json read in
+  // Should be Pythia6, LUND or HEPMC3
+  if (obj.isMember("OutputType")){
+    gOutputType = obj["OutputType"].asString();
+    if (gOutputType == "Pythia6"){
+      cout << "Using Pythia6 output format for Fun4All" << endl;
+    }
+    else if (gOutputType == "LUND"){
+      cout << "Using LUND output format" << endl;
+    }
+    else if (gOutputType == "HEPMC3"){
+      cout << "Using HEPMC3 output format for ePIC" << endl;
+    }
+    else{
+      cout << "Output type not recognised!" << endl;
+      cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl;
+      gOutputType = "HEPMC3";
+    }
+  }
+  else{
+    cout << "Output type not specified in .json file!" << endl;
+    cout << "Setting output type to HEPMC3 (ePIC) by default!" << endl;
+    gOutputType = "HEPMC3";
+  }
+  ///*--------------------------------------------------*/
+  /// The detector selection is determined here
+  /// The incidence proton phi angle is 
+  if (obj.isMember("det_location")){
+    gDet_location = obj["det_location"].asString();
+    if (gDet_location == "ip8") {
+      fProton_incidence_phi = 0.0;
+    } 
+    else if (gDet_location == "ip6") {
+      fProton_incidence_phi = fPi;
+    }
     else {
-      gBeamPart = "Proton";
+      fProton_incidence_phi = 0.0;
+      cout << "The interaction point requested is not recognized!" << endl;
+      cout << "Therefore ip6 is used by default." << endl;
     }
+  }
+  else{ // 21/12/22 - This could probably be combined with the else statement above in some way
+    fProton_incidence_phi = 0.0;
+    cout << "The interaction points was not specified in the .json file!" << endl;
+    cout << "Therefore ip6 is used by default" << endl;
+  }
 
-	// SJDK - 12/01/22
-	// Set output type as a .json read in
-	// Should be Pythia6, LUND or HEPMC3
-	if (obj.isMember("OutputType")){
-	  gOutputType = obj["OutputType"].asString();
-	  if (gOutputType == "Pythia6"){
-	    cout << "Using Pythia6 output format for Fun4All" << endl;
-	  }
-	  else if (gOutputType == "LUND"){
-	    cout << "Using LUND output format" << endl;
-	  }
-	  else if (gOutputType == "HEPMC3"){
-	    cout << "Using HEPMC3 output format for EPIC" << endl;
-	  }
-	  else{
-	    cout << "Output type not recognised!" << endl;
-	    cout << "Setting output type to HEPMC3 by default!" << endl;
-	    gOutputType = "HEPMC3";
-	  }
-	}
-	else{
-	  cout << "Output type not specified in .json file!" << endl;
-	  cout << "Setting output type to HEPMC3 by default!" << endl;
-	  gOutputType = "HEPMC3";
-	}
-	///*--------------------------------------------------*/
-	/// The detector selection is determined here
-	/// The incidence proton phi angle is 
-	if (obj.isMember("det_location")){
-	  gDet_location = obj["det_location"].asString();
-	  if (gDet_location == "ip8") {
-	    fProton_incidence_phi = 0.0;
-	  } 
-	  else if (gDet_location == "ip6") {
-	    fProton_incidence_phi = fPi;
-	  }
-	  else {
-	    fProton_incidence_phi = 0.0;
-	    cout << "The interaction point requested is not recognized!" << endl;
-	    cout << "Therefore ip6 is used by default." << endl;
-	  }
-	}
-	else{ // 21/12/22 - This could probably be combined with the else statement above in some way
-	    fProton_incidence_phi = 0.0;
-	    cout << "The interaction points was not specified in the .json file!" << endl;
-	    cout << "Therefore ip6 is used by default" << endl;
-	}
-
-	if (obj.isMember("Ee_Low")){
-	  fScatElec_E_Lo = obj["Ee_Low"].asDouble();
-	}
-	else{
-	  fScatElec_E_Lo = 0.5;
-	  cout << "Minumum scattered electron energy not specified in .json file, defaulting to 0.5*EBeam." << endl;
-	}
+  if (obj.isMember("Ee_Low")){
+    fScatElec_E_Lo = obj["Ee_Low"].asDouble();
+  }
+  else{
+    fScatElec_E_Lo = 0.5;
+    cout << "Minumum scattered electron energy not specified in .json file, defaulting to 0.5*EBeam." << endl;
+  }
 
-	if (obj.isMember("Ee_High")){
-	  fScatElec_E_Hi = obj["Ee_High"].asDouble();
-	}
-	else{
-	  fScatElec_E_Hi = 2.5;
-	  cout << "Max scattered electron energy not specified in .json file, defaulting to 2.5*EBeam." << endl;
-	}
+  if (obj.isMember("Ee_High")){
+    fScatElec_E_Hi = obj["Ee_High"].asDouble();
+  }
+  else{
+    fScatElec_E_Hi = 2.5;
+    cout << "Max scattered electron energy not specified in .json file, defaulting to 2.5*EBeam." << endl;
+  }
 
-	if (obj.isMember("e_Theta_Low")){
-	  fScatElec_Theta_I = obj["e_Theta_Low"].asDouble() * fDEG2RAD;
-	}
-	else{
-	  fScatElec_Theta_I = 60.0 * fDEG2RAD;
-	  cout << "Min scattered electron theta not specified in .json file, defaulting to 60 degrees." << endl;
-	}
+  if (obj.isMember("e_Theta_Low")){
+    fScatElec_Theta_I = obj["e_Theta_Low"].asDouble() * fDEG2RAD;
+  }
+  else{
+    fScatElec_Theta_I = 60.0 * fDEG2RAD;
+    cout << "Min scattered electron theta not specified in .json file, defaulting to 60 degrees." << endl;
+  }
 
-	if (obj.isMember("e_Theta_High")){
-	  fScatElec_Theta_F = obj["e_Theta_High"].asDouble() * fDEG2RAD;
-	}
-	else{
-	  fScatElec_Theta_F = 175.0 * fDEG2RAD;
-	  cout << "Max scattered electron theta not specified in .json file, defaulting to 175 degrees." << endl;
-	}
+  if (obj.isMember("e_Theta_High")){
+    fScatElec_Theta_F = obj["e_Theta_High"].asDouble() * fDEG2RAD;
+  }
+  else{
+    fScatElec_Theta_F = 175.0 * fDEG2RAD;
+    cout << "Max scattered electron theta not specified in .json file, defaulting to 175 degrees." << endl;
+  }
 
-	if (obj.isMember("EjectileX_Theta_Low")){
-	  fEjectileX_Theta_I = obj["EjectileX_Theta_Low"].asDouble() * fDEG2RAD;
-	}
-	else{
-	  fEjectileX_Theta_I = 0.0 * fDEG2RAD;
-	  cout << "Min ejectile X theta not specified in .json file, defaulting to 0 degrees." << endl;
-	}
+  if (obj.isMember("EjectileX_Theta_Low")){
+    fEjectileX_Theta_I = obj["EjectileX_Theta_Low"].asDouble() * fDEG2RAD;
+  }
+  else{
+    fEjectileX_Theta_I = 0.0 * fDEG2RAD;
+    cout << "Min ejectile X theta not specified in .json file, defaulting to 0 degrees." << endl;
+  }
 
-	if (obj.isMember("EjectileX_Theta_High")){
-	  fEjectileX_Theta_F = obj["EjectileX_Theta_High"].asDouble() * fDEG2RAD;
-	}
-	else{
-	  fEjectileX_Theta_F = 60.0 * fDEG2RAD;
-	  cout << "Max ejectile X theta not specified in .json file, defaulting to 60 degrees." << endl;
-	}
-        
-	SigPar = ReadCrossSectionPar(particle, hadron);
+  if (obj.isMember("EjectileX_Theta_High")){
+    fEjectileX_Theta_F = obj["EjectileX_Theta_High"].asDouble() * fDEG2RAD;
+  }
+  else{
+    fEjectileX_Theta_F = 60.0 * fDEG2RAD;
+    cout << "Max ejectile X theta not specified in .json file, defaulting to 60 degrees." << endl;
+  }
 
-        if(particle != "pi0"){ // Default case now
-	  Reaction* r1 = new Reaction(particle, hadron);
-	  r1->process_reaction();
-	  delete r1;
-	}
-	else{  // Treat pi0 slightly differently for now
-	  Reaction* r1 = new Reaction(particle);
-	  r1->process_reaction();
-	  delete r1;
-	}
+  // 06/09/23 - SJDK - Added string to check method chosen
+  TString CalcMethod = obj["calc_method"].asString(); 
+  if(CalcMethod == "Analytical"){
+    UseSolve = false;
+  }
+  else if (CalcMethod == "Solve"){
+    UseSolve = true;
+  }
+  else{
+    cout << "! WARNING !" << endl  << "! WARNING !- Calculation method not specified or not recognised, defaulting to Analytical - ! WARNING!" << endl << "! WARNING !" << endl;
+    UseSolve = false;
+  }
+	
+  SigPar = ReadCrossSectionPar(Ejectile, RecoilHadron);
+	
+  if(Ejectile != "pi0"){ // Default case now
+    Reaction* r1 = new Reaction(Ejectile, RecoilHadron);
+    r1->process_reaction();
+    delete r1;
+  }
+  else{  // Treat pi0 slightly differently for now
+    Reaction* r1 = new Reaction(Ejectile);
+    r1->process_reaction();
+    delete r1;
+  }
 }
 
 /*--------------------------------------------------*/
 /*--------------------------------------------------*/
 
 void SetEICSeed (int seed) {
-	fSeed = seed;
+  fSeed = seed;
 }
 
 ///*--------------------------------------------------*/
@@ -344,27 +363,27 @@ void SetEICSeed (int seed) {
 
 TString ExtractParticle(TString particle) {
 
-	/// Make the input particle case insansitive
-	particle.ToLower();
-	if (particle.Contains("on")) {
-		particle.ReplaceAll("on", "");
-	};
+  /// Make the input particle case insansitive
+  particle.ToLower();
+  if (particle.Contains("on")) {
+    particle.ReplaceAll("on", "");
+  };
  	
-	if (particle.Contains("plus")) {
-		particle.ReplaceAll("plus", "+");
-	}
+  if (particle.Contains("plus")) {
+    particle.ReplaceAll("plus", "+");
+  }
 
-	if (particle.Contains("minus")) {
-		particle.ReplaceAll("minus", "-");
-	}
+  if (particle.Contains("minus")) {
+    particle.ReplaceAll("minus", "-");
+  }
 
-	if (particle.Contains("zero")) {
-		particle.ReplaceAll("zero", "0");
-	}
+  if (particle.Contains("zero")) {
+    particle.ReplaceAll("zero", "0");
+  }
 
-	particle[0] = toupper(particle[0]);
-	cout << "Particle: " << particle << endl;
-	return particle;
+  particle[0] = toupper(particle[0]);
+  cout << "Particle: " << particle << endl;
+  return particle;
 
 }
 
@@ -385,31 +404,30 @@ TString ExtractCharge(TString particle) {
   return charge;
 }
 
-vector<vector<vector<vector<double>>>> ReadCrossSectionPar(TString particle, TString hadron){
+vector<vector<vector<vector<double>>>> ReadCrossSectionPar(TString EjectileX, TString RecoilHad){
   
   string sigL_ParamFile, sigT_ParamFile;
  
-  if (particle == "Pi+" && hadron == "Neutron"){
-    cout << "Add Pi+/Neutron case here" << endl;
+  if (EjectileX == "Pi+" && RecoilHad == "Neutron"){
+    // When pion model parameterised in some way, add Pi+/Neutron case here - 
   }
-  else if (particle == "Pi-" && hadron == "Proton"){
-    cout << "Add Pi-/Proton case here" << endl;
+  else if (EjectileX == "Pi-" && RecoilHad == "Proton"){
+    // When pion model parameterised in some way, add Pi-/Proton case here - 
   }
-  else if (particle == "K+" && hadron == "Lambda"){
-    cout << "Add K+/Lambda case here" << endl;
+  else if (EjectileX == "K+" && RecoilHad == "Lambda"){
     sigL_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigL";
     sigT_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusLambda_Param_sigT"; // Shouldn't really have a relative path, should look at setting a DEMPGen variable and doing this in a better way later
   }
-  else if (particle == "K+" && hadron == "Sigma0"){
-    cout << "Add K+/Sigma case here" << endl;
+  else if (EjectileX == "K+" && RecoilHad == "Sigma0"){
     sigL_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigL";
     sigT_ParamFile = "../src/eic_evgen/CrossSection_Params/KPlusSigma_Param_sigT";
   }
-  else if (particle == "Pi0"){
-    cout << "Add Pi0 case here" << endl;
+  else if (EjectileX == "Pi0"){
+    // When pi0 model parameterised, add it here
   }
   else{
-    cout << "Throw some error" << endl;
+    cerr << " !!!!! " << endl << "Warning!" << endl << "Combination of specified ejectile and recoil hadron not found!" << "Cross section parameters cannot be read, check inputs!" << endl << "Warning!" << endl << " !!!!! " << endl;
+    exit(-1);
   }
  
   //....................................................................................................
diff --git a/src/eic_evgen/eic.h b/src/eic_evgen/eic.h
index 81b036d..c4847bf 100644
--- a/src/eic_evgen/eic.h
+++ b/src/eic_evgen/eic.h
@@ -58,7 +58,7 @@ void SetEICSeed(int);
 TString ExtractParticle(TString);
 TString ExtractCharge(TString);
 
-vector<vector<vector<vector<double>>>> ReadCrossSectionPar(TString particle, TString hadron);
+vector<vector<vector<vector<double>>>> ReadCrossSectionPar(TString EjectileX, TString RecoilHad);
 
 #endif
 
diff --git a/src/eic_evgen/eic_pim.cc b/src/eic_evgen/eic_pim.cc
index 14c527f..93346ff 100644
--- a/src/eic_evgen/eic_pim.cc
+++ b/src/eic_evgen/eic_pim.cc
@@ -24,6 +24,7 @@ TTree *t1;
 
 int gKinematics_type;
 bool gPi0_decay;
+bool UseSolve;
 string gDet_location;
 string gOutputType; // SJDK 12/01/22 - Added output type as a variable you can specify in the .json file
 string gBeamPart; // SJDK 12/01/22 - Added output type as a variable you can specify in the .json file
@@ -38,14 +39,15 @@ int fWLessShell, fWLess1P9, fSDiff;
 
 //long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNSigmaNeg, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile; 
 
-unsigned long long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNSigmaNeg, fNaN, fConserve, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile;
+unsigned long long int fNEvents, fNRecorded, fNGenerated, fWSqNeg, fNSigmaNeg, fNaN, fConserve, fNWeightUnphys, fNWeightReject, fLundRecorded, fNFile, fSolveEvents_0Sol, fSolveEvents_1Sol, fSolveEvents_2Sol;
 
-// SJDK 03/04/23 - Added in Qsq Min/Max and W Min/Max
-double fK, fm, fElectron_Kin_Col_GeV, fElectron_Kin_Col, fRand, fLumi, fuBcm2, fPI, fDEG2RAD, fRAD2DEG, fEBeam, fPBeam, fScatElec_Theta_I, fScatElec_Theta_F, fPion_Theta_I, fPion_Theta_F, fEjectileX_Theta_I, fEjectileX_Theta_F, fScatElec_E_Hi, fScatElec_E_Lo, fPSF, fQsq_Min, fQsq_Max, fW_Min, fW_Max; 
+// SJDK 03/04/23 - Added in Qsq Min/Max, W Min/Max and t max (06/09/23)
+// 13/09/23 - SJDK - New generic HBeam value (rather than proton beam)
+double fK, fm, fElectron_Kin_Col_GeV, fElectron_Kin_Col, fRand, fLumi, fuBcm2, fPI, fDEG2RAD, fRAD2DEG, fEBeam, fPBeam, fHBeam, fScatElec_Theta_I, fScatElec_Theta_F, fPion_Theta_I, fPion_Theta_F, fEjectileX_Theta_I, fEjectileX_Theta_F, fScatElec_E_Hi, fScatElec_E_Lo, fPSF, fQsq_Min, fQsq_Max, fW_Min, fW_Max, fT_Max; 
 
 double fOmega_Theta_I, fOmega_Theta_F, fOmega_Theta_Col, fOmega_Phi_Col;
 
-double fDiff_E, conserve, ene, mom, ene_mom, mom_px, mom_py, mom_pz, mom_pxpy, mom_pxpz, mom_pypz, mom_pxpypz;     // 18/06/21 AU -> New variables to count envents passing/not passing conservation laws
+double fDiff_E, conserve, ene, mom, ene_mom, mom_px, mom_py, mom_pz, mom_pxpy, mom_pxpz, mom_pypz, mom_pxpypz; // 18/06/21 AU -> New variables to count envents passing/not passing conservation laws
 
 double fMandSConserve, fTop_Pion_Mom, fBot_Pion_Mom, fPion_Mom_Same, fEnergyConserve, fXMomConserve, fYMomConserve, fZMomConserve, fXMomConserve_RF, fYMomConserve_RF, fZMomConserve_RF, fEnergyConserve_RF; 
 
@@ -239,7 +241,6 @@ void pim::Initilize() {
     // 18/01/23 - The luminosity below is some default assumtpion, more up to date values are set in DEMP prod and depend upon beam energy combinations if they are specified
     // See slide 11 in https://indico.cern.ch/event/1072579/contributions/4796856/attachments/2456676/4210776/CAP-EIC-June-7-2022-Seryi-r2.pdf for more info
     // fLumi                                     = 0.374e33; // Jlab design
-    //fLumi                                       = 1e34; // https://eic.jlab.org/wiki/index.php/EIC_luminosity - OUTDATED
     fLumi                                       = 1e33; // 18/01/23, this seems a better default based upon more up to date info, see link above
     fuBcm2                                      = 1.0e-30;
     fPI                                         = 3.1415926;
@@ -278,14 +279,12 @@ void pim::Initilize() {
     fRecoilProton_Mass_GeV                      = fRecoilProton_Mass/1000.0;
     fPion_Mass                                  = 139.57039 ;
     fPion_Mass_GeV                              = fPion_Mass/1000.0;
-
     fKaon_Mass                                  = 493.677;
     fKaon_Mass_GeV                              = fKaon_Mass/1000.0;
     fLambda_Mass                                = 1115.683;
     fLambda_Mass_GeV                            = fLambda_Mass/1000.0;
     fSigma_Mass                                 = 1192.642;
     fSigma_Mass_GeV                             = fSigma_Mass/1000.0;
-
     fOmega_Mass                                 = 782.66;
     fOmega_Mass_GeV                             = fOmega_Mass/1000.0;
 
@@ -332,6 +331,9 @@ void pim::Initilize() {
     fConserve                                   = 0;
     fNWeightUnphys                              = 0;
     fNWeightReject                              = 0;
+    fSolveEvents_0Sol                           = 0;
+    fSolveEvents_1Sol                           = 0;
+    fSolveEvents_2Sol                           = 0;
     fSDiff                                      = 0;
     fScatElecEnergyLess                         = 0;
     fScatElecThetaLess                          = 0;
diff --git a/src/eic_evgen/eic_pim.h b/src/eic_evgen/eic_pim.h
index 99f8fc1..06d43e7 100644
--- a/src/eic_evgen/eic_pim.h
+++ b/src/eic_evgen/eic_pim.h
@@ -50,6 +50,7 @@ extern TString gfile_name;
 extern TString gParticle;
 extern TString gHadron;
 extern bool gPi0_decay;
+extern bool UseSolve;
 extern std::string gDet_location;
 extern std::string gOutputType;
 extern std::string gBeamPart;
@@ -115,6 +116,10 @@ extern unsigned long long int fConserve;
 extern unsigned long long int fNWeightUnphys;
 extern unsigned long long int fNWeightReject;
 
+extern unsigned long long int fSolveEvents_0Sol;
+extern unsigned long long int fSolveEvents_1Sol;
+extern unsigned long long int fSolveEvents_2Sol;
+
 extern unsigned long long int fNFile;
 
 extern double fK;
@@ -129,6 +134,8 @@ extern double fDEG2RAD;
 extern double fRAD2DEG;
 extern double fEBeam;
 extern double fPBeam;
+// 13/09/23 - SJDK - New generic HBeam value (rather than proton beam)
+extern double fHBeam;
 extern double fScatElec_Theta_I;
 extern double fScatElec_Theta_F;
 extern double fPion_Theta_I; // SJDK 19/12/22 - These should be removed in future, specific to pion reaction cases. Should be generic MesonX
@@ -143,6 +150,7 @@ extern double fQsq_Min;
 extern double fQsq_Max;
 extern double fW_Min;
 extern double fW_Max;
+extern double fT_Max;
 
 extern double fMandSConserve;
 extern double fTop_Pion_Mom;
diff --git a/src/eic_evgen/process_routine/DEMP_Reaction.cc b/src/eic_evgen/process_routine/DEMP_Reaction.cc
index 595eab8..9446a10 100644
--- a/src/eic_evgen/process_routine/DEMP_Reaction.cc
+++ b/src/eic_evgen/process_routine/DEMP_Reaction.cc
@@ -54,7 +54,7 @@ void DEMP_Reaction::process_reaction() {
     rNEvent_itt = i;
     fNGenerated ++;
  
-    Progress_Report();  // This is happens at each 10% of the total event is processed
+    Progress_Report();  // This happens at each 10% of the total event is processed
     Processing_Event();
   }
  
@@ -93,13 +93,12 @@ void DEMP_Reaction::Init() {
  
   dRootTree->Branch("EventWeight", &fEventWeight, "fEventWeight/D");
 
-
   /*--------------------------------------------------*/ 
 	
   qsq_ev = 0, t_ev = 0, w_neg_ev = 0, w_ev = 0;
   rNEvents = fNEvents;
   rNEvent_itt = 0;
-
+  
   // 02/06/21 SJDK 
   // Set these values once the beam energies are read in
   fPSF = ( fEBeam * ( fScatElec_E_Hi - fScatElec_E_Lo ) *( sin( fScatElec_Theta_F ) - sin( fScatElec_Theta_I ) ) * 2 * fPI *( sin( fEjectileX_Theta_F ) - sin( fEjectileX_Theta_I ) ) * 2 * fPI );
@@ -107,12 +106,14 @@ void DEMP_Reaction::Init() {
   fElectron_Kin_Col = fElectron_Kin_Col_GeV * 1000.0;
 
   // cout << rNEvents << "    " << fNEvents << endl;
-	
-  rFermiMomentum = pd->fermiMomentum();
+
+  // 08/09/23 - SJDK - Fermi momentum commented out for now, this is not fully implemented yet
+  // In future, this will be enabled/disabled automatically depending upon the specified hadron beam
+  //rFermiMomentum = pd->fermiMomentum();
 
   // ----------------------------------------------------
   // Proton in collider (lab) frame
-
+  // 08/09/23 - SJDK - The naming needs to be adjusted to be the incoming hadron beam, not the proton. Make this generic.
   r_lproton = GetProtonVector_lab();
   r_lprotong = GetProtonVector_lab() * fm;
 
@@ -128,13 +129,14 @@ void DEMP_Reaction::Init() {
   // ----------------------------------------------------
   // Electron in collider (lab) frame
 
-  cout << "Fermi momentum: " << rFermiMomentum << endl;
+  // 06/09/23 - SJDK - Commenting out for now, should be disabled for e/p collisions
+  //cout << "Fermi momentum: " << rFermiMomentum << endl;
 
   r_lelectron	 = GetElectronVector_lab();
   r_lelectrong = r_lelectron * fm;
 
   ///*--------------------------------------------------*/
-  /// Getting the particle mass from the data base
+  /// Getting the ejectile (produced meson) particle mass from the data base
  
   produced_X = ParticleEnum(rEjectile);
   f_Ejectile_Mass = ParticleMass(produced_X)*1000; //MeV
@@ -176,33 +178,75 @@ void DEMP_Reaction::Init() {
   f_Ejectile_Theta_F = fEjectileX_Theta_F;
 
   cout << "Produced particle in exclusive production: " << rEjectile << ";  with mass: " << f_Ejectile_Mass << " MeV "<< endl;
-  cout << fEBeam << " GeV electrons on " << fPBeam << " GeV ions" << endl;
-  
+  cout << fEBeam << " GeV electrons on " << fHBeam << " GeV ions" << endl;
+  if(UseSolve == true){
+    cout << rEjectile << " and " << rEjectile_scat_hadron << " 4-vectors calculated using Solve function" << endl;
+  }
+  else if(UseSolve == false){
+    cout << rEjectile << " and " << rEjectile_scat_hadron << " 4-vectors calculated using analytical solution" << endl;
+  }
   // Set luminosity value based upon beam energy combination, note that if no case matches, a default of 1e33 is assumed. Cases are a set of nominal planned beam energy combinations for the EIC (and EICC)
   // See slide 11 in https://indico.cern.ch/event/1072579/contributions/4796856/attachments/2456676/4210776/CAP-EIC-June-7-2022-Seryi-r2.pdf
   // If available in the future, this could be replaced by some fixed function
-  if ((fEBeam == 5.0 ) && (fPBeam == 41.0) ){
+  if ((fEBeam == 5.0 ) && (fHBeam == 41.0) ){
     fLumi = 0.44e33;
   }
-  else if ((fEBeam == 5.0 ) && (fPBeam == 100.0) ){
+  else if ((fEBeam == 5.0 ) && (fHBeam == 100.0) ){
     fLumi = 3.68e33;
   }
-  else if ((fEBeam == 10.0 ) && (fPBeam == 100.0) ){
+  else if ((fEBeam == 10.0 ) && (fHBeam == 100.0) ){
     fLumi = 4.48e33;
   }
-  else if ((fEBeam == 18.0 ) && (fPBeam == 275.0) ){
+  else if ((fEBeam == 18.0 ) && (fHBeam == 275.0) ){
     fLumi = 1.54e33;
   }
-  else if ((fEBeam == 3.5 ) && (fPBeam == 20) ){ // EICC optimal beam energy combination
+  else if ((fEBeam == 3.5 ) && (fHBeam == 20) ){ // EICC optimal beam energy combination
     fLumi = 2e33;
   }
-  else if ((fEBeam == 2.8 ) && (fPBeam == 13) ){ // EICC lowest beam energy combination
+  else if ((fEBeam == 2.8 ) && (fHBeam == 13) ){ // EICC lowest beam energy combination
     fLumi = 0.7e33;
   }
   else{
     cout << "!!! Notice !!! The beam energy combination simulated does not match an expected case, a default luminosity value of - " << fLumi << " cm^2s^-1 has been assumed. !!! Notice !!!" << endl;
   }
- 
+
+  if(UseSolve == true){
+    /*--------------------------------------------------*/ 
+    // SJDK 03/04/22 -  New set of initialisation stuff for the solve function from Ishan and Bill
+    
+    CoinToss = new TRandom3();
+
+    F = new TF1("F",
+		"[6]-sqrt([7]**2+x**2)-sqrt([8]**2+([3]-[0]*x)**2+([4]-[1]*x)**2+([5]-[2]*x)**2)",
+		0, r_lproton.E());
+
+    char AngleGenName[100] = "AngleGen";
+    double dummy[2] = {0,1};
+
+    f_Ejectile_Theta_I = fEjectileX_Theta_I;
+    f_Ejectile_Theta_F = fEjectileX_Theta_F;
+
+    double ThetaRange[2] = {f_Ejectile_Theta_I, f_Ejectile_Theta_F};
+    double PhiRange[2] = {0, 360*TMath::DegToRad()};
+
+    AngleGen = new CustomRand(AngleGenName, dummy,
+			      ThetaRange, PhiRange);
+
+    UnitVect = new TVector3(0,0,1);
+
+    //  ///*--------------------------------------------------*/ 
+    //  // Produced hadron and recoilded hadron from the solve function 
+
+    VertBeamElec = new TLorentzVector();
+    VertScatElec = new TLorentzVector();
+
+    Initial      = new TLorentzVector();
+    Target       = new TLorentzVector();
+    Photon       = new TLorentzVector();
+    Interaction  = new TLorentzVector();
+    Final        = new TLorentzVector();
+  }
+  
 }
 
 void DEMP_Reaction::Processing_Event() {
@@ -210,10 +254,10 @@ void DEMP_Reaction::Processing_Event() {
   // ----------------------------------------------------
   // Considering Fermi momentum for the proton
   // ----------------------------------------------------
-  // SJDK - 31/01/23 - This doesn't seem to do anything?
-  if( kCalcFermi ) {
-    Consider_Proton_Fermi_Momentum(); 
-  }
+  // SJDK - 06/06/23 - Commenting out for now, this increases the number of recorded events - Should have no fermi momentum for e/p collisions
+  // if( kCalcFermi ) {
+  //   Consider_Proton_Fermi_Momentum(); 
+  // }
 
   // ----------------------------------------------------
   // Boost vector from collider (lab) frame to protons rest frame (Fix target)
@@ -230,10 +274,9 @@ void DEMP_Reaction::Processing_Event() {
   fScatElec_Energy_Col = fRandom->Uniform( fScatElec_E_Lo * fElectron_Energy_Col , fScatElec_E_Hi * fElectron_Energy_Col );
 
   // ----------------------------------------------------
-  // Produced Particle X in Collider frame
+  // Produced ejectile in Collider frame
   // ----------------------------------------------------  
 
-  /// The generic produced particle in the exclusive reaction is labelled as X 
   f_Ejectile_Theta_Col      = acos( fRandom->Uniform( cos(f_Ejectile_Theta_I), cos(f_Ejectile_Theta_F ) ) ); 
   f_Ejectile_Phi_Col        = fRandom->Uniform( 0 , 2.0 * fPi );
     	
@@ -283,62 +326,72 @@ void DEMP_Reaction::Processing_Event() {
     return;
   }
 
-  // ---------------------------------------------------------
-  // Pion momentum in collider frame, analytic solution starts
-  // ---------------------------------------------------------
+  // SJDK - 06/09/23 - Check UseSolved boolean, process through relevant loop
+  if (UseSolve == false){
+    // ---------------------------------------------------------
+    // Pion momentum in collider frame, analytic solution starts
+    // ---------------------------------------------------------
  
-  double fupx = sin( f_Ejectile_Theta_Col ) * cos( f_Ejectile_Phi_Col );
-  double fupy = sin( f_Ejectile_Theta_Col ) * sin( f_Ejectile_Phi_Col );
-  double fupz = cos( f_Ejectile_Theta_Col );
+    double fupx = sin( f_Ejectile_Theta_Col ) * cos( f_Ejectile_Phi_Col );
+    double fupy = sin( f_Ejectile_Theta_Col ) * sin( f_Ejectile_Phi_Col );
+    double fupz = cos( f_Ejectile_Theta_Col );
  
-  double fuqx = sin( r_lphoton.Theta() ) * cos( r_lphoton.Phi() );
-  double fuqy = sin( r_lphoton.Theta() ) * sin( r_lphoton.Phi() );
-  double fuqz = cos( r_lphoton.Theta() );
+    double fuqx = sin( r_lphoton.Theta() ) * cos( r_lphoton.Phi() );
+    double fuqy = sin( r_lphoton.Theta() ) * sin( r_lphoton.Phi() );
+    double fuqz = cos( r_lphoton.Theta() );
  
-  double fa = -(r_lphoton.Vect()).Mag() * ( fupx * fuqx +  fupy * fuqy +  fupz * fuqz );
-  double fb = pow ( (r_lphoton.Vect()).Mag() , 2 );
-  double fc = r_lphoton.E() + fProton_Mass;
+    double fa = -(r_lphoton.Vect()).Mag() * ( fupx * fuqx +  fupy * fuqy +  fupz * fuqz );
+    double fb = pow ( (r_lphoton.Vect()).Mag() , 2 );
+    double fc = r_lphoton.E() + fProton_Mass;
  
-  fa = ( fa - std::abs( (r_lproton.Vect()).Mag() ) * ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fupx ) + 
-  						       ( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fupy ) + 
-  						       ( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fupz ) ) );
+    fa = ( fa - std::abs( (r_lproton.Vect()).Mag() ) * ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fupx ) + 
+							 ( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fupy ) + 
+							 ( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fupz ) ) );
      
-  double factor = ( pow( (r_lproton.Vect()).Mag() , 2 ) + 2.0 * (r_lphoton.Vect()).Mag() * (r_lproton.Vect()).Mag() *  
-  		    ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fuqx ) + 
-  		      ( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fuqy ) + 
-  		      ( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fuqz ) ) );
+    double factor = ( pow( (r_lproton.Vect()).Mag() , 2 ) + 2.0 * (r_lphoton.Vect()).Mag() * (r_lproton.Vect()).Mag() *  
+		      ( ( ( r_lproton.X() / (r_lproton.Vect()).Mag() ) * fuqx ) + 
+			( ( r_lproton.Y() / (r_lproton.Vect()).Mag() ) * fuqy ) + 
+			( ( r_lproton.Z() / (r_lproton.Vect()).Mag() ) * fuqz ) ) );
      
-  fb =  fb + factor;  
-  fc = r_lphoton.E() + r_lproton.E();
+    fb =  fb + factor;  
+    fc = r_lphoton.E() + r_lproton.E();
      
-  double ft = fc * fc - fb + f_Ejectile_Mass * f_Ejectile_Mass - f_Recoil_Mass * f_Recoil_Mass;
+    double ft = fc * fc - fb + f_Ejectile_Mass * f_Ejectile_Mass - f_Recoil_Mass * f_Recoil_Mass;
      
-  double fQA = 4.0 * ( fa * fa - fc * fc );
-  double fQB = 4.0 * fc * ft;
+    double fQA = 4.0 * ( fa * fa - fc * fc );
+    double fQB = 4.0 * fc * ft;
 
-  double fQC = -4.0 * fa * fa * f_Ejectile_Mass * f_Ejectile_Mass - ft * ft;    
+    double fQC = -4.0 * fa * fa * f_Ejectile_Mass * f_Ejectile_Mass - ft * ft;    
  
-  fradical = fQB * fQB - 4.0 * fQA * fQC;
+    fradical = fQB * fQB - 4.0 * fQA * fQC;
  
-  fepi1 = ( -fQB - sqrt( fradical ) ) / ( 2.0 * fQA );
-  fepi2 = ( -fQB + sqrt( fradical ) ) / ( 2.0 * fQA );
-
-  ///---------------------------------------------------------
-  /// Particle X momentum in collider frame, analytic solution
-  /// And obtain recoiled proton in collider (lab) frame
-  ///---------------------------------------------------------
-
-   r_l_Ejectile.SetPxPyPzE( (sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * cos(f_Ejectile_Phi_Col),
-  		   ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * sin(f_Ejectile_Phi_Col),
-  		   ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * cos(f_Ejectile_Theta_Col),
-  		   fepi1 );
+    fepi1 = ( -fQB - sqrt( fradical ) ) / ( 2.0 * fQA );
+    fepi2 = ( -fQB + sqrt( fradical ) ) / ( 2.0 * fQA );
+
+    ///---------------------------------------------------------
+    /// Particle X momentum in collider frame, analytic solution
+    /// And obtain recoiled proton in collider (lab) frame
+    ///---------------------------------------------------------
+
+    r_l_Ejectile.SetPxPyPzE( (sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * cos(f_Ejectile_Phi_Col),
+			     ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * sin(f_Ejectile_Theta_Col) * sin(f_Ejectile_Phi_Col),
+			     ( sqrt( pow( fepi1 , 2) - pow(f_Ejectile_Mass , 2) ) ) * cos(f_Ejectile_Theta_Col),
+			     fepi1 );
+  
+    l_Recoil.SetPxPyPzE( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile).X(),
+			 ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Y(),
+			 ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Z(),
+			 sqrt( pow( ( ( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Vect() ).Mag()),2) +
+			       pow( f_Recoil_Mass , 2) ) );
+  }
+   else if (UseSolve == true){
+     if(!Solve()){
+       return;
+     }  
+     r_l_Ejectile = r_l_Ejectile_solved;
+     l_Recoil = r_l_Recoil_solved;
+   }
   
-   l_Recoil.SetPxPyPzE( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile).X(),
-  			       ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Y(),
-  			       ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Z(),
-  			       sqrt( pow( ( ( ( r_lproton + r_lelectron - r_lscatelec - r_l_Ejectile ).Vect() ).Mag()),2) +
-  				     pow( f_Recoil_Mass , 2) ) );
-
   ///--------------------------------------------------
   
   r_l_Ejectile_g = r_l_Ejectile * fm;
@@ -351,41 +404,24 @@ void DEMP_Reaction::Processing_Event() {
   r_lw = r_lproton + r_lphoton;
   fW = r_lw.Mag();
 
-  // ----------------------------------------------------------------------------------------------
-  // Calculate w prime w' = (proton + photon - pion)^2                                             
-  // ----------------------------------------------------------------------------------------------
- 
-  lwp = r_lprotong + r_lphotong - r_l_Ejectile_g;
-  fW_Prime_GeV = lwp.Mag();    
-
-  fsini = r_lelectron + r_lproton;
-  fsfin = r_lscatelec + r_l_Ejectile + l_Recoil;
-     
-  fsinig = fsini * fm;
-  fsfing = fsfin * fm; 
-  // SJDK 15/06/21 - Mandlestam S conservation check - doesn't actually seem to be utilised?
-  fMandSConserve = std::abs( fsinig.Mag() - fsfing.Mag() );
-
   // SJDK 15/06/21 - Added integer counters for conservation law check and for NaN check
   if (r_l_Ejectile.E() != r_l_Ejectile.E()){ // SJDK 15/06/21 - If the energy of the produced meson is not a number, return and add to counter
     fNaN++;
     return;
   }
-  kSConserve = false;
-  if( std::abs( fsinig.Mag() - fsfing.Mag() ) < fDiff ) {
-    kSConserve = true;
-  }
-
-///*--------------------------------------------------*/ 
-//-> 10/05/23 - Love added a slimmed down, simpler to read version of the CheckLaws fn
-// 
-// To check the conservation of the energy and momentum, there two methods avalaible:
-// Method 1: Give the four-vectors of the initial and final states partciles, 
-//           tolerance factor will be defaulted 1e-6 MeV
-// Method 2: Give the four-vectors of the initial and final states partciles, 
-//           and the prefered tolerance factor.
-//
 
+  ///*--------------------------------------------------*/ 
+  //-> 10/05/23 - Love added a slimmed down, simpler to read version of the CheckLaws fn
+  // 
+  // To check the conservation of the energy and momentum, there two methods avalaible:
+  // Method 1: Give the four-vectors of the initial and final states partciles, 
+  //           tolerance factor will be defaulted 1e-6 MeV
+  //           CheckLaws(e_beam, h_beam, scatt_e, ejectile, recoil) <- input 4 vectors
+  // Method 2: Give the four-vectors of the initial and final states partciles, 
+  //           and the prefered tolerance factor.
+  //           CheckLaws(e_beam, h_beam, scatt_e, ejectile, recoil, tolerance) <- input 4 vectors and tolerance value in GeV
+  // Both functions return 1 if conservations laws are satisified
+  
    if( pd->CheckLaws(r_lelectron, r_lproton, r_lscatelec, r_l_Ejectile, l_Recoil) !=1 ){
      fConserve++;
      return;
@@ -461,16 +497,9 @@ void DEMP_Reaction::Processing_Event() {
   }
   */ 
 
-  // 31/01/23 SJDK - New limit on t, remove only events outside the parameterisation range, limits depend upon particle typ
-  if (rEjectile == "Pi+" && fT_GeV > 1.3 ) {
-    t_ev++;
-    return;
-  }
-  else if (rEjectile == "K+" && fT_GeV > 2.0) {
-    t_ev++;
-    return;
-  }
-  else if (rEjectile == "Pi0+" && fT_GeV > 0.5){ // 03/02/23 - SJDK - Not sure what range is used for pi0, assume < 0.5 for now, would be u in this case anyway?
+  // 31/01/23 SJDK - New limit on t, remove only events outside the parameterisation range
+  // 06/09/23 SJDK - fT_Max set in eic.cc depending upon ejectile type
+  if (fT_GeV > fT_Max ) {
     t_ev++;
     return;
   }
@@ -637,7 +666,7 @@ TLorentzVector DEMP_Reaction::GetProtonVector_lab() {
   //     fProton_Phi_Col   = fPi; 
   fProton_Phi_Col   = fProton_incidence_phi; 
 
-  fProton_Mom_Col   = fPBeam * 1e3; 
+  fProton_Mom_Col   = fHBeam * 1e3; 
   fVertex_X         = 0.; 
   fVertex_Y         = 0.; 
   fVertex_Z         = 0.; 
@@ -754,32 +783,67 @@ Double_t DEMP_Reaction::Get_Total_Cross_Section() {
 
 /*--------------------------------------------------*/
 /// Output generator detail
-
+// 06/09/23 SJDK - Formatting of these is all messed up annoyingly, would be nice to get the final numbers to align. They don't currently.
+// Cuts are now ordered as they are applied in the generator
 void DEMP_Reaction::Detail_Output() {
 
-  DEMPDetails << "Total events tried                                           " << setw(20) << fNGenerated   << endl;
-  DEMPDetails << "Total events recorded                                        " << setw(20) << fNRecorded    << endl;
-  DEMPDetails << "Number of events with wsq negative                           " << setw(20) << w_neg_ev      << endl;
-  DEMPDetails << "Number of events with " << fW_Min << " < w < " << fW_Max       << "                           " << setw(20) << w_ev << endl;
-  DEMPDetails << "Number of events with " << fQsq_Min << " < qsq < " << fQsq_Max << "                           " << setw(20) << qsq_ev << endl;
-  DEMPDetails << "Number of events with Meson (X) energy NaN                   " << setw(20) << fNaN          << endl;
-  DEMPDetails << "Total events passing conservation law check with tolerance " << fDiff << setw(17) << conserve   << endl;
-  DEMPDetails << "Total events failing conservation law checks                 " << setw(20) << fConserve     << endl;
-  DEMPDetails << "Total events failing energy conservation check ONLY          " << setw(20) << ene   << endl; 
-  DEMPDetails << "Total events failing momentum conservation check ONLY        " << setw(20) << mom   << endl;
-  DEMPDetails << "Total events failing energy AND momentum conservation checks " << setw(20) << ene_mom   << endl;
-  DEMPDetails << "Total events failing px conservation law check               " << setw(20) << mom_px   << endl;
-  DEMPDetails << "Total events failing py conservation law check               " << setw(20) << mom_py   << endl;
-  DEMPDetails << "Total events failing pz conservation law check               " << setw(20) << mom_pz   << endl;
-  DEMPDetails << "Total events failing px and py conservation law checks       " << setw(20) << mom_pxpy   << endl;
-  DEMPDetails << "Total events failing px and pz conservation law checks       " << setw(20) << mom_pxpz   << endl;
-  DEMPDetails << "Total events failing py and pz conservation law checks       " << setw(20) << mom_pypz   << endl;
-  DEMPDetails << "Total events failing px, py and pz conservation law checks   " << setw(20) << mom_pxpypz   << endl;
-  DEMPDetails << "Number of events with -t > 2 (K+) or -t > 1.3 (Pi+) GeV      " << setw(20) << t_ev          << endl;
-  DEMPDetails << "Number of events with w less than threshold                  " << setw(20) << fWSqNeg       << endl;
-  DEMPDetails << "Number of events with Sigma negative                         " << setw(20) << fNSigmaNeg    << endl;
   DEMPDetails << "Seed used for the Random Number Generator                    " << setw(20) << fSeed         << endl;
+  DEMPDetails << endl;
+  DEMPDetails << "Total events tried                                           " << setw(20) << fNGenerated   << endl;
+  if(UseSolve == true){
+    DEMPDetails << "Total events cut                                             " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fSolveEvents_0Sol) << setw(20) << ((double) (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fSolveEvents_0Sol)/(double)fNGenerated)*100 << " %" << endl;
+    DEMPDetails << "Total events recorded                                        " << setw(20) << fNRecorded << setw(20) << ((double)fNRecorded/(double)fNGenerated)*100 << " %" << endl;
+    if (fNGenerated != (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fNRecorded + fSolveEvents_0Sol)){
+      DEMPDetails << "Total events cut + recorded = events tried?                " << setw(20) << "NO! ERROR!" << endl;
+    }
+    else{
+      DEMPDetails << "Total events cut + recorded = events tried?                " << setw(22) << "Yes! :)" << endl;
+    }
 
+  }
+  else{
+    DEMPDetails << "Total events cut                                             " << setw(20) << (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg) << setw(20) << ((double) (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg)/(double)fNGenerated)*100 << " %" << endl;
+    DEMPDetails << "Total events recorded                                        " << setw(20) << fNRecorded << setw(20) << ((double)fNRecorded/(double)fNGenerated)*100 << " %" << endl;
+    if (fNGenerated != (qsq_ev + w_ev + w_neg_ev + fNaN + fConserve + t_ev + fNSigmaNeg + fNRecorded)){
+      DEMPDetails << "Total events cut + recorded = events tried?                " << setw(20) << "NO! ERROR!" << endl;
+    }
+    else{
+      DEMPDetails << "Total events cut + recorded = events tried?                " << setw(22) << "Yes! :)" << endl;
+    }
+  }
+  
+  DEMPDetails << endl << "Cut details -" << endl;
+  DEMPDetails << "Events cut due to qsq < " << fQsq_Min << " or qsq > "<< fQsq_Max << "                        " << setw(20) << qsq_ev << setw(20) << ((double)qsq_ev/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to negative Wsq value                         " << setw(20) << w_neg_ev << setw(20) << ((double)w_neg_ev/(double)fNGenerated)*100 << " %" << endl;  
+  DEMPDetails << "Events cut due to W < " << fW_Min << " or W > " << fW_Max << "                          " << setw(20) << w_ev << setw(20) << ((double)w_ev/(double)fNGenerated)*100 << " %" << endl;
+  if(UseSolve == true){
+    DEMPDetails << "Events cut due to solve function finding 0 solutions         " << setw(20) << fSolveEvents_0Sol << setw(20) << ((double)fSolveEvents_0Sol/(double)fNGenerated)*100 << " %" << endl;
+  }
+  DEMPDetails << "Events cut due to ejectile (X) energy NaN                    " << setw(20) << fNaN << setw(20) << ((double)fNaN/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to conservation law check failure             " << setw(20) << fConserve << setw(20) << ((double)fConserve/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to -t > " << fT_Max << "GeV                      " << setw(30) << t_ev << setw(20) << ((double)t_ev/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to -ve cross section value                    " << setw(20) << fNSigmaNeg << setw(20) << ((double)fNSigmaNeg/(double)fNGenerated)*100 << " %" << endl;
+
+  DEMPDetails << endl << "Conservation law checks details -" << endl;
+  DEMPDetails << "Total events PASSING conservation law check with tolerance " << fDiff << setw(17) << conserve << endl;
+  DEMPDetails << "Events cut due to energy conservation check ONLY             " << setw(20) << ene << setw(20) << ((double)ene/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to momentum conservation check ONLY           " << setw(20) << mom << setw(20) << ((double)mom/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to energy AND momentum conservation checks    " << setw(20) << ene_mom << setw(20) << ((double)ene_mom/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to px conservation law check                  " << setw(20) << mom_px << setw(20) << ((double)mom_px/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to py conservation law check                  " << setw(20) << mom_py << setw(20) << ((double)mom_py/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to pz conservation law check                  " << setw(20) << mom_pz << setw(20) << ((double)mom_pz/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to px and py conservation law checks          " << setw(20) << mom_pxpy << setw(20) << ((double)mom_pxpy/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to px and pz conservation law checks          " << setw(20) << mom_pxpz << setw(20) << ((double)mom_pxpz/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to py and pz conservation law checks          " << setw(20) << mom_pypz << setw(20) << ((double)mom_pypz/(double)fNGenerated)*100 << " %" << endl;
+  DEMPDetails << "Events cut due to px, py and pz conservation law checks      " << setw(20) << mom_pxpypz << setw(20) << ((double)mom_pxpypz/(double)fNGenerated)*100 << " %" << endl;
+
+  if(UseSolve == true){
+    DEMPDetails << endl << "Solve function, addtional info -" << endl;
+    DEMPDetails << "Number of events with 0 Solution                             " << setw(20) << fSolveEvents_0Sol << endl;
+    DEMPDetails << "Number of events with 1 Solution                             " << setw(20) << fSolveEvents_1Sol << endl;
+    DEMPDetails << "Number of events with 2 Solution                             " << setw(20) << fSolveEvents_2Sol << endl;
+  }
+  
 }
 
 ////*--------------------------------------------------
@@ -1002,7 +1066,7 @@ void DEMP_Reaction::DEMPReact_HEPMC3_Out_Init() {
 
 void DEMP_Reaction::DEMPReact_HEPMC3_Output() {
   
-  // HEPMC3 output for Athena/EPIC simulations
+  // HEPMC3 output for Athena/ePIC simulations
 
   // First line - E - Event# - #Vertices - #Particles
   DEMPOut << std::scientific << std::setprecision(15) << "E" << " "  << print_itt <<  " " << "1" << " " << 5 << endl;
@@ -1025,3 +1089,167 @@ void DEMP_Reaction::DEMPReact_HEPMC3_Output() {
   DEMPOut << "P" << " " << "5" << " " << "-1" << " " << PDGtype(recoil_hadron) << " " << l_Recoil_g.X() << " "  << l_Recoil_g.Y() << " "  << l_Recoil_g.Z() << " " << l_Recoil_g.E() << " " <<  l_Recoil_g.M() << " " << "1" << endl;
 
 }
+
+/*--------------------------------------------------*/ 
+
+bool DEMP_Reaction::SolnCheck()
+{
+  //
+  //  // Double Checking for solution viability
+  //  if (TMath::Abs(f_Scat_hadron_Mass-r_l_scat_hadron_solved->M())>1){
+  //    //cerr << "Mass Missmatch" << endl;
+  //    //cerr << TMath::Abs(proton_mass_mev-Proton->M()) << endl;
+  //    return false;
+  //  }
+  //  if (TMath::Abs(W_in()-W_out())>1){
+  //    //cerr << "W Missmatch" << endl;
+  //    //cerr << TMath::Abs(W_in()-W_out()) << endl;
+  //    return false;
+  //  }
+  //  *Final = *r_l_scat_hadron_solved + *r_lX_solved;
+  //
+  //  if (TMath::Abs(Initial->Px()-Final->Px())>1){
+  //    //cerr << "Px Missmatch" << endl;
+  //    //cerr << TMath::Abs(Initial->Px()-Final->Px()) << endl;
+  //    return false;
+  //  }
+  //
+  //  if (TMath::Abs(Initial->Py()-Final->Py())>1){
+  //    //cerr << "Py Missmatch" << endl;
+  //    //cerr << TMath::Abs(Initial->Py()-Final->Py()) << endl;
+  //    return false;
+  //  }
+  //
+  //  if (TMath::Abs(Initial->Pz()-Final->Pz())>1){
+  //    //cerr << "Pz Missmatch" << endl;
+  //    //cerr << TMath::Abs(Initial->Pz()-Final->Pz()) << endl;
+  //    return false;
+  //  }
+  //
+  //  if (TMath::Abs(Initial->E()-Final->E())>1){
+  //    return false;
+  //  }
+  return true;
+}
+
+/*--------------------------------------------------*/ 
+double DEMP_Reaction::W_in()
+{
+  return 0;
+}
+
+/*--------------------------------------------------*/ 
+double DEMP_Reaction::W_out()
+{
+  return 0;
+}
+
+/*--------------------------------------------------*/ 
+
+int DEMP_Reaction::Solve()
+{
+
+
+  VertBeamElec->SetPxPyPzE(r_lelectron.Px(), r_lelectron.Py(), r_lelectron.Pz(), r_lelectron.E());
+  VertScatElec->SetPxPyPzE(r_lscatelec.Px(), r_lscatelec.Py(), r_lscatelec.Pz(), r_lscatelec.E());
+  Target->SetPxPyPzE(r_lproton.Px(), r_lproton.Py(), r_lproton.Pz(), r_lproton.E());
+  *Photon = *VertBeamElec - *VertScatElec;
+  *Interaction = *Photon;
+
+  *Initial = *Interaction+*Target;
+
+  theta =  f_Ejectile_Theta_Col;
+  phi   =  f_Ejectile_Phi_Col;  
+
+  return this->Solve(theta, phi);
+}
+
+
+int DEMP_Reaction::Solve(double theta, double phi)
+{
+
+  W_in_val = W_in();
+
+  if (W_in_val<0){
+    return 0;
+  }
+
+  UnitVect->SetTheta(theta);
+  UnitVect->SetPhi(phi);
+  UnitVect->SetMag(1);
+
+  double* pars = new double[9];
+
+  pars[0] = UnitVect->X();
+  pars[1] = UnitVect->Y();
+  pars[2] = UnitVect->Z();
+  pars[3] = Initial->Px();
+  pars[4] = Initial->Py();
+  pars[5] = Initial->Pz();
+  pars[6] = Initial->E();
+  pars[7] = f_Ejectile_Mass;
+  pars[8] = f_Recoil_Mass;
+
+  F->SetParameters(pars);
+
+
+  ///*--------------------------------------------------*/ 
+  // Looking for the 1st Solution:
+  //    If a solution found, then this will be the fist solution. Then we proceed to look for the 2nd solution. 
+  //    If no soluion found, then exit solve function
+
+  P = F->GetX(0, 0, pars[6], 0.0001, 10000);
+
+  if (TMath::Abs(F->Eval(P)) < 1){
+    fSolveEvents_1Sol++;
+  } else {
+    fSolveEvents_0Sol++; 
+    return 0;
+  }
+
+  TLorentzVector * r_l_Ejectile_solved_1_temp = new TLorentzVector();
+  TLorentzVector * r_l_Ejectile_solved_2_temp = new TLorentzVector();
+
+  Float_t r_l_Ejectile_E = sqrt( pow(P*pars[0],2) + pow(P*pars[1],2) + pow(P*pars[2],2) + pow(f_Ejectile_Mass,2) );
+  r_l_Ejectile_solved_1_temp->SetPxPyPzE(P*pars[0], P*pars[1],  P*pars[2], r_l_Ejectile_E);
+
+  ///*--------------------------------------------------*/ 
+  // Looking for the 2nd Solution 
+
+  P2 = F->GetX(0, P+100, pars[6], 0.0001, 10000);
+  Float_t r_l_Ejectile_E_2 = sqrt( pow(P2 * pars[0],2) + pow(P2 * pars[1],2) + pow(P2 * pars[2],2) + pow(f_Ejectile_Mass,2) );
+  r_l_Ejectile_solved_2_temp->SetPxPyPzE(P2 * pars[0], P2 * pars[1],  P2 * pars[2], r_l_Ejectile_E_2);
+
+  ///*--------------------------------------------------*/ 
+  // If a valid 2nd solution is found, then we are certian that there are two solutions.
+  //   - We then increament the counter for 2nd solution scenario
+  //   - We then decreament the counter for the 1st solution scenario 
+
+  if (TMath::Abs(F->Eval(P2)) < 1){
+    fSolveEvents_2Sol++;
+    fSolveEvents_1Sol--;
+    if ( Int_t(CoinToss->Uniform(0,100)) < 50) {
+      r_l_Ejectile_solved.SetPxPyPzE(r_l_Ejectile_solved_1_temp->X(), r_l_Ejectile_solved_1_temp->Y(), r_l_Ejectile_solved_1_temp->Z(), r_l_Ejectile_solved_1_temp->E());
+    } else {
+      r_l_Ejectile_solved.SetPxPyPzE(r_l_Ejectile_solved_2_temp->X(), r_l_Ejectile_solved_2_temp->Y(), r_l_Ejectile_solved_2_temp->Z(), r_l_Ejectile_solved_2_temp->E());
+    }
+  }
+  else {
+    r_l_Ejectile_solved.SetPxPyPzE(r_l_Ejectile_solved_1_temp->X(), r_l_Ejectile_solved_1_temp->Y(), r_l_Ejectile_solved_1_temp->Z(), r_l_Ejectile_solved_1_temp->E());
+  }
+
+  ///*--------------------------------------------------*/ 
+  /// Solve for the recoil information with the "solved" Ejectile informaiton
+  TLorentzVector * r_l_hadron_temp= new TLorentzVector();
+  *r_l_hadron_temp = *Initial- r_l_Ejectile_solved;
+  r_l_Recoil_solved.SetPxPyPzE(r_l_hadron_temp->Px(), r_l_hadron_temp->Py(), r_l_hadron_temp->Pz(), r_l_hadron_temp->E());
+
+  delete r_l_Ejectile_solved_1_temp;
+  delete r_l_Ejectile_solved_2_temp;
+
+  delete r_l_hadron_temp;
+  delete[] pars;
+ 
+  return 1;
+
+}
diff --git a/src/eic_evgen/reaction_routine.cc b/src/eic_evgen/reaction_routine.cc
index 90ece05..0bec8e6 100644
--- a/src/eic_evgen/reaction_routine.cc
+++ b/src/eic_evgen/reaction_routine.cc
@@ -19,12 +19,11 @@ using namespace std;
 
 /*--------------------------------------------------*/
 /// Reaction 
+Reaction::Reaction(TString ejectile_str) { 
 
-Reaction::Reaction(TString particle_str) { 
-
- 	rEjectile = particle_str;
-	cout << "Produced particle is: " << GetParticle() << endl; 
-	cout << "Generated process: e + p -> e' + p' + " << GetParticle() << endl; 
+ 	rEjectile = ejectile_str;
+	cout << "Produced ejectile is: " << GetEjectile() << endl; 
+	cout << "Generated process: e + p -> e' + p' + " << GetEjectile() << endl; 
     	tTime.Start(); 
  
  	cout << "/*--------------------------------------------------*/" << endl;
@@ -37,13 +36,13 @@ Reaction::Reaction(TString particle_str) {
 }
 
 // SJDK 09/02/22 - New reaction where the particle and hadron are specified
-Reaction::Reaction(TString particle_str, TString hadron_str) { 
+Reaction::Reaction(TString ejectile_str, TString recoil_hadron_str) { 
 
- 	rEjectile = particle_str;
-	rRecoil = hadron_str;
-	cout << "Produced particle is: " << GetParticle() << endl; 
-	cout << "Produced hadron is: " << GetHadron() << endl;
-	cout << "Generated process: e + p -> e'+ " << GetHadron() << " + " << GetParticle() << endl; 
+ 	rEjectile = ejectile_str;
+	rRecoil = recoil_hadron_str;
+	cout << "Produced ejectile is: " << GetEjectile() << endl; 
+	cout << "Produced recoil hadron is: " << GetRecoilHadron() << endl;
+	cout << "Generated process: e + p -> e'+ " << GetRecoilHadron() << " + " << GetEjectile() << endl; 
     	tTime.Start(); 
  
  	cout << "/*--------------------------------------------------*/" << endl;
@@ -78,17 +77,10 @@ Reaction::~Reaction() {
 ///
 
 void Reaction::process_reaction() {
-//  if (rEjectile == "Pi0") {
-//    //		Pi0_Production* r1 = new Pi0_Production("Eta");
-//    Pi0_Production* rr1 = new Pi0_Production(rEjectile);
-//    rr1->process_reaction();
-//    delete rr1;
-//  }
-//  else{
+
   // SJDK - 19/12/22 - New generic DEMP reaction class, the intention is that this should be able to handle any case
     DEMP_Reaction* rr1 = new DEMP_Reaction(rEjectile, rRecoil);
     rr1->process_reaction();
     delete rr1;
-//  }
 
 }
diff --git a/src/eic_evgen/reaction_routine.h b/src/eic_evgen/reaction_routine.h
index c8e388e..8b3c7b0 100644
--- a/src/eic_evgen/reaction_routine.h
+++ b/src/eic_evgen/reaction_routine.h
@@ -29,8 +29,8 @@ class Reaction{
   ~Reaction();
 
   void process_reaction();		
-  TString GetParticle() {return rEjectile;};		
-  TString GetHadron() {return rRecoil;};
+  TString GetEjectile() {return rEjectile;};		
+  TString GetRecoilHadron() {return rRecoil;};
   
  protected:
   TStopwatch tTime;
@@ -48,8 +48,8 @@ class DEMP_Reaction {
   ~DEMP_Reaction();
 
   void process_reaction();		
-  TString GetParticle() {return rEjectile;};		
-  TString GetHadron() {return rRecoil;};
+  TString GetEjectile() {return rEjectile;};		
+  TString GetRecoilHadron() {return rRecoil;};
 
  protected:
 
@@ -74,7 +74,7 @@ class DEMP_Reaction {
 
   Double_t Get_Total_Cross_Section(); 
 
-//  Double_t GetPi0_CrossSection();
+  //  Double_t GetPi0_CrossSection();
 
   /*--------------------------------------------------*/
   // Parameters
@@ -87,7 +87,7 @@ class DEMP_Reaction {
 
   std::string sTFile;   /// Generator output files. For documentation and monitoring purposes 
   std::string sLFile;   /// Lund input file into the EIC simulation
-  std::string sDFile;   /// Root dianostic plot in root file format
+  std::string sDFile;   /// Root diagnostic plot in root file format
 
   std::ofstream DEMPOut;     
   std::ofstream DEMPDetails;
@@ -137,11 +137,8 @@ class DEMP_Reaction {
   TLorentzVector r_l_Ejectile;
   TLorentzVector r_l_Ejectile_g;
 
-//  Particle* r_l_Ejectile_solved;
-//  Particle* l_Recoil_solved;
-
   TLorentzVector r_l_Ejectile_solved;
-  TLorentzVector l_Recoil_solved;
+  TLorentzVector r_l_Recoil_solved;
 
   double f_Ejectile_Mass;
   double f_Ejectile_Mass_GeV;
@@ -220,28 +217,19 @@ class DEMP_Reaction {
 
   ///*--------------------------------------------------*/ 
   // Rory Check algorithm
+
+  TLorentzVector* Interaction_Solve;
+  TLorentzVector* Target_Solve;
+
+  TLorentzVector* VertBeamElec;
+  TLorentzVector* VertScatElec;
   
-  //  Particle* Interaction;
-  //  Particle* Target;
-  //
-  //  Particle* Initial;
-  //  Particle* Final;
-  //
-  //  Particle* VertBeamElec;
-  //  Particle* VertScatElec;
-  //  Particle* Photon;
-
-
-  TLorentzVector Interaction_Solve;
-  TLorentzVector Target_Solve;
-  
-  TLorentzVector Initial;
-  TLorentzVector Final;
+  TLorentzVector* Initial;
+  TLorentzVector* Target;
+  TLorentzVector* Photon;
+  TLorentzVector* Interaction;
+  TLorentzVector* Final;
   
-  TLorentzVector VertBeamElec;
-  TLorentzVector VertScatElec;
-  TLorentzVector Photon;
-
   bool SolnCheck();
   double W_in_Solve(); 
   double W_out_Solve();
@@ -256,6 +244,9 @@ class DEMP_Reaction {
   int Solve();
   int Solve(double theta, double phi);
 
+  double W_in();
+  double W_out();
+  
   ///*--------------------------------------------------*/ 
   // Needed for the Solve function