diff --git a/NPLib/Detectors/MUST2/TMust2Physics.cxx b/NPLib/Detectors/MUST2/TMust2Physics.cxx
index b82b7d03f89db45b4267076896e023fe54da6fa8..e03492f3b77aadbd32a8b2e6646b802599b3c6b1 100644
--- a/NPLib/Detectors/MUST2/TMust2Physics.cxx
+++ b/NPLib/Detectors/MUST2/TMust2Physics.cxx
@@ -1232,11 +1232,12 @@ void TMust2Physics::AddParameterToCalibrationManager() {
vector<double> standardCsI;
for (int i = 0; i < m_NumberOfTelescope; ++i) {
- if (m_CsIOffset[i] == 1) {
+ if (m_CsIOffset[i+1] == 1) {
standardCsI = {0, 500. / 16384.};
}
- else
+ else{
standardCsI = {-250, 250. / 8192.};
+ }
for (int j = 0; j < 128; ++j) {
Cal->AddParameter("MUST2", "T" + NPL::itoa(i + 1) + "_Si_X" + NPL::itoa(j + 1) + "_E",
diff --git a/NPSimulation/Detectors/MUST2/MUST2Array.cc b/NPSimulation/Detectors/MUST2/MUST2Array.cc
index ec47e723258b3c572ef3c83058cbdc7108e8f6a3..831edab0fd9535d954542613cfdf110a5642eb6a 100644
--- a/NPSimulation/Detectors/MUST2/MUST2Array.cc
+++ b/NPSimulation/Detectors/MUST2/MUST2Array.cc
@@ -933,7 +933,7 @@ void MUST2Array::ReadSensitive(const G4Event*) {
vector<unsigned int> level = CsIScorer->GetLevel(i);
if (ECsI > ThresholdCsI) {
if (m_CsIOffset[level[0]] == 1)
- m_Event->SetCsIE(level[0], level[1], NPL::EnergyToADC(ECsI, 0, 0, 0, 16384));
+ m_Event->SetCsIE(level[0], level[1], NPL::EnergyToADC(ECsI, 0, 500, 0, 16384));
else
m_Event->SetCsIE(level[0], level[1], NPL::EnergyToADC(ECsI, 0, 250, 8192, 16384));
double timeCsI = RandGauss::shoot(CsIScorer->GetTime(i), ResoTimeMust);
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index d1fd91c7c81adb7f00a804ef6356ce82dc078dab..2ac1d414d37b4272fc65f3a2861088d6eaebcf8f 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -115,6 +115,7 @@ void NPS::BeamReaction::ReadConfiguration() {
else if (m_ReactionType == Fusion) {
vector<InputBlock*> blocks = input.GetAllBlocksWithToken("FusionReaction");
m_BeamName = NPL::ChangeNameToG4Standard(blocks[0]->GetString("Beam"));
+ m_BeamNameNPL = blocks[0]->GetString("Beam");
m_TargetNuclei = blocks[0]->GetString("Target");
m_FusionProduct = blocks[0]->GetString("Product");
m_FusionExcitation = blocks[0]->GetDouble("ExcitationEnergy", "MeV");
@@ -676,6 +677,8 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
NPL::Particle N(m_FusionProduct);
N.SetExcitationEnergy(m_FusionExcitation);
NPL::Particle T(m_TargetNuclei);
+ NPL::Particle B(m_BeamNameNPL);
+
int PZ = N.GetZ();
int PA = N.GetA();
Product = IonTable->GetIon(PZ, PA, m_FusionExcitation * MeV);
@@ -684,7 +687,7 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
TVector3 BeamP = NPS::ConvertVector(PrimaryTrack->GetMomentum());
TLorentzVector BeamLV;
- BeamLV.SetVectM(BeamP, N.Mass() * MeV);
+ BeamLV.SetVectM(BeamP, B.Mass() * MeV);
TLorentzVector TargetLV;
TargetLV.SetVectM(TVector3(0, 0, 0), T.Mass() * MeV);
TLorentzVector TotalLV = BeamLV + TargetLV;
@@ -692,6 +695,13 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
// energy lost in the fusion process to be removed to the total energy
// Total Available Ek = Initial Ek + (InitialMass-FinalMass)
double KineAvailable = TotalLV.Et() + (TotalLV.Mag() - N.Mass());
+ // beam energy too low to allow for fusion to occur
+ if (KineAvailable < 0) {
+ // FIXME
+ // avoid the creation of secondary but the primary is still killed
+ return;
+ }
+
G4ThreeVector momentum_dir = NPS::ConvertVector(TotalLV.Vect().Unit());
//////FIXME Unsure of this part
diff --git a/NPSimulation/Process/BeamReaction.hh b/NPSimulation/Process/BeamReaction.hh
index de59daf728e9277ea2de1bc0fabbc71fc7d0a3ad..9d1bf2d57693aa7087f611893e4fc0ca2f122594 100644
--- a/NPSimulation/Process/BeamReaction.hh
+++ b/NPSimulation/Process/BeamReaction.hh
@@ -24,68 +24,64 @@
#ifndef BeamReaction_h
#define BeamReaction_h
-#include "G4VFastSimulationModel.hh"
#include "G4Abla.hh"
#include "G4AblaInterface.hh"
#include "G4Fragment.hh"
-#include "PhysicsList.hh"
-#include "NPReaction.h"
-#include "NPQFS.h"
+#include "G4VFastSimulationModel.hh"
#include "NPPhaseSpace.h"
+#include "NPQFS.h"
+#include "NPReaction.h"
+#include "PhysicsList.hh"
#include "TReactionConditions.h"
class G4VPhysicalVolume;
-namespace NPS{
- enum ReactionType{
- TwoBody,
- QFS,
- PhaseSpace,
- Fusion
- };
+namespace NPS {
+ enum ReactionType { TwoBody, QFS, PhaseSpace, Fusion };
- class BeamReaction : public G4VFastSimulationModel{
- public:
- BeamReaction (G4String, G4Region*);
- BeamReaction (G4String);
- ~BeamReaction ();
+ class BeamReaction : public G4VFastSimulationModel {
+ public:
+ BeamReaction(G4String, G4Region*);
+ BeamReaction(G4String);
+ ~BeamReaction();
+
+ public:
+ void ReadConfiguration();
+ G4bool IsApplicable(const G4ParticleDefinition&);
+ G4bool ModelTrigger(const G4FastTrack&);
+ void DoIt(const G4FastTrack&, G4FastStep&);
+
+ private:
+ NPL::Reaction m_Reaction;
+ NPL::QFS m_QFS;
+ NPL::PhaseSpace m_PhaseSpace;
+ string m_BeamName;
+ string m_BeamNameNPL;
+ int m_ReactionType;
+ G4AblaInterface* ABLA;
+
+ bool m_active; // is the process active
+ bool m_shoot;
+ double m_StepSize;
+ double m_Z;
+ double m_S;
+ double m_rand;
+ double m_length;
+ int m_Parent_ID;
+ double m_event_weight;
+ double SlowDownBeam(const G4ParticleDefinition* Beam, double IncidentEnergy, double Thickness,
+ G4Material* Material);
- public:
- void ReadConfiguration();
- G4bool IsApplicable(const G4ParticleDefinition&);
- G4bool ModelTrigger(const G4FastTrack &);
- void DoIt(const G4FastTrack&, G4FastStep&);
-
- private:
- NPL::Reaction m_Reaction;
- NPL::QFS m_QFS;
- NPL::PhaseSpace m_PhaseSpace;
- string m_BeamName;
- int m_ReactionType;
- G4AblaInterface* ABLA;
+ private: // specific for the simple case of fusion
+ string m_TargetNuclei;
+ string m_FusionProduct;
+ double m_FusionExcitation;
- bool m_active;// is the process active
- bool m_shoot;
- double m_StepSize;
- double m_Z;
- double m_S;
- double m_rand;
- double m_length;
- int m_Parent_ID;
- double m_event_weight;
- double SlowDownBeam(const G4ParticleDefinition* Beam, double IncidentEnergy, double Thickness,G4Material* Material);
-
- private:// specific for the simple case of fusion
- string m_TargetNuclei;
- string m_FusionProduct;
- double m_FusionExcitation;
-
private:
- TReactionConditions* m_ReactionConditions;
-
+ TReactionConditions* m_ReactionConditions;
+
public:
void AttachReactionConditions();
- void SetStepSize(double step){m_StepSize=step;};
+ void SetStepSize(double step) { m_StepSize = step; };
};
-}
-
+} // namespace NPS
-#endif
+#endif
diff --git a/Projects/e870/DetectorConfiguration/MUST2_E805.detector b/Projects/e870/DetectorConfiguration/MUST2_E805.detector
new file mode 100644
index 0000000000000000000000000000000000000000..4c78f7563ad7c4d4783c632bee6c0a61c96cc24f
--- /dev/null
+++ b/Projects/e870/DetectorConfiguration/MUST2_E805.detector
@@ -0,0 +1,60 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Target
+ %THICKNESS= 26 micrometer
+ THICKNESS= 0.000001 micrometer
+ ANGLE= 0 deg
+ RADIUS= 15 mm
+ MATERIAL= CH2
+ ANGLE= 0 deg
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+%%%%%%% Telescope 1 %%%%%%%
+M2Telescope
+ X1_Y1= -13.57 -104.78 299.83 mm
+ X1_Y128= -25.1 -12.62 328.63 mm
+ X128_Y1= -104.01 -104.85 263.88 mm
+ X128_Y128= -115.53 -12.71 292.67 mm
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+ CsIOffset= 1
+
+%%%%%%% Telescope 2 %%%%%%%
+M2Telescope
+ X1_Y1= -114.9 9.68 291.84 mm
+ X1_Y128= -24.56 9.8 327.88 mm
+ X128_Y1= -103.58 101.8 262.73 mm
+ X128_Y128= -13.24 102.09 298.75 mm
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+ CsIOffset= 1
+
+%%%%%%% Telescope 3 %%%%%%%
+M2Telescope
+ X1_Y1= 11.84 101.72 299.05 mm
+ X1_Y128= 23.44 9.68 328.38 mm
+ X128_Y1= 102.16 101.57 263.01 mm
+ X128_Y128= 113.81 9.57 292.3 mm
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+ CsIOffset= 1
+
+
+%%%%%%% Telescope 4 %%%%%%%
+M2Telescope
+ X1_Y1= 113.56 -13.18 292.11 mm
+ X1_Y128= 23.23 -13.37 328.15 mm
+ X128_Y1= 102.39 -105.49 263.59 mm
+ X128_Y128= 12.04 -105.69 299.63 mm
+ SI= 1
+ SILI= 0
+ CSI= 1
+ VIS= all
+ CsIOffset= 1
+