diff --git a/NPLib/Physics/TReactionConditions.cxx b/NPLib/Physics/TReactionConditions.cxx
index d5969c63a8e5eed0dcce556c884e52d2b1c92b16..a4bd2d87427d2fc0317548cd42cb9fbdfce0ea4b 100644
--- a/NPLib/Physics/TReactionConditions.cxx
+++ b/NPLib/Physics/TReactionConditions.cxx
@@ -54,6 +54,8 @@ void TReactionConditions::Clear(){
fRC_Particle_Name.clear();
fRC_Theta.clear();
fRC_Phi.clear();
+ fRC_LAB_Theta.clear();
+ fRC_LAB_Phi.clear();
fRC_Kinetic_Energy.clear();
fRC_Momentum_Direction_X.clear();
fRC_Momentum_Direction_Y.clear();
diff --git a/NPLib/Physics/TReactionConditions.h b/NPLib/Physics/TReactionConditions.h
index 0fd75266dcb6c4aa59bdd442bdbf8b35039eb069..f54017e894f74a41a4cddd376d733a2cde7da3fb 100644
--- a/NPLib/Physics/TReactionConditions.h
+++ b/NPLib/Physics/TReactionConditions.h
@@ -63,6 +63,8 @@ private:
vector<string> fRC_Particle_Name;
vector<double> fRC_Theta;
vector<double> fRC_Phi;
+ vector<double> fRC_LAB_Theta;
+ vector<double> fRC_LAB_Phi;
vector<double> fRC_Kinetic_Energy;
vector<double> fRC_Momentum_Direction_X;
vector<double> fRC_Momentum_Direction_Y;
@@ -97,6 +99,8 @@ public:
void SetParticleName (const string & Particle_Name) {fRC_Particle_Name.push_back(Particle_Name);}//!
void SetTheta (const double & Angle) {fRC_Theta.push_back(Angle);}//!
void SetPhi (const double & AnglePhi) {fRC_Phi.push_back(AnglePhi);}//!
+ void SetLABTheta (const double & LABAngle) {fRC_LAB_Theta.push_back(LABAngle);}//!
+ void SetLABPhi (const double & LABAnglePhi) {fRC_LAB_Phi.push_back(LABAnglePhi);}//!
void SetKineticEnergy (const double & Kinetic_Energy) {fRC_Kinetic_Energy.push_back(Kinetic_Energy);}//!
void SetMomentumDirectionX (const double & Momentum_Direction_X) {fRC_Momentum_Direction_X.push_back(Momentum_Direction_X);}//!
void SetMomentumDirectionY (const double & Momentum_Direction_Y) {fRC_Momentum_Direction_Y.push_back(Momentum_Direction_Y);}//!
diff --git a/NPSimulation/Detectors/Dali/Dali.cc b/NPSimulation/Detectors/Dali/Dali.cc
index 9b0bbb4aeaf62851734706300cc004ba5851b14e..bbf2fa3c32561259e383a8472cd7c41818f2d491 100644
--- a/NPSimulation/Detectors/Dali/Dali.cc
+++ b/NPSimulation/Detectors/Dali/Dali.cc
@@ -73,7 +73,7 @@ namespace Dali_NS{
// Energy and time Resolution
const double EnergyThreshold = 0*MeV;
const double ResoTime = 0.0*ns; //4.5*ns ;
- const double ResoEnergy = 0.001*MeV ;
+ const double ResoEnergy = 1.36*MeV ; // mean Resolution(FWHM) 1.7% of 80MeV from slides 20170214-SAMURAI34-setup-DALI.pdf // 0.001*MeV ;
const double Radius = 50*mm ;
const double Width = 49.76*mm ;
const double Hight = 84.81*mm ;
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 5f4d79daa05116898f4b60d4a5ac54a9b98ad9fa..1cb6d413a91b25920cfb6c9c86f52936491a2f5f 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -229,7 +229,7 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep)
G4ThreeVector U(1,0,0);
G4ThreeVector V(0,1,0);
-
+ G4ThreeVector ZZ(0,0,1);
m_ReactionConditions->SetBeamEmittanceTheta(PrimaryTrack->GetMomentumDirection().theta()/deg);
m_ReactionConditions->SetBeamEmittancePhi(PrimaryTrack->GetMomentumDirection().phi()/deg);
m_ReactionConditions->SetBeamEmittanceThetaX(PrimaryTrack->GetMomentumDirection().angle(U)/deg);
@@ -238,17 +238,24 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep)
///// Build rotation matrix to go from the incident //////
///// beam frame to the "world" frame //////
//////////////////////////////////////////////////////////
+
+ /*
G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
cos(Beam_theta) * sin(Beam_phi),
-sin(Beam_theta));
- G4ThreeVector col2(-sin(Beam_phi),
- cos(Beam_phi),
+ G4ThreeVector col2(-sin(-Beam_phi),
+ cos(-Beam_phi),
0);
G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
sin(Beam_theta) * sin(Beam_phi),
cos(Beam_theta));
+
G4RotationMatrix BeamToWorld(col1, col2, col3);
+
+ */
+
+
/////////////////////////////////////////////////////////////////
///// Angles for emitted particles following Cross Section //////
///// Angles are in the beam frame //////
@@ -257,7 +264,7 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep)
// Angles
// Shoot and Set a Random ThetaCM
m_Reaction.ShootRandomThetaCM();
- double phi = RandFlat::shoot() * 2. * pi;
+ double phi = RandFlat::shoot() * 2. * pi;
//////////////////////////////////////////////////
///// Momentum and angles from kinematics /////
@@ -272,18 +279,20 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep)
G4ThreeVector momentum_kine3_beam(sin(Theta3) * cos(phi),
sin(Theta3) * sin(phi),
cos(Theta3));
- // Momentum in World frame
- G4ThreeVector momentum_kine3_world = BeamToWorld * momentum_kine3_beam;
-
+ // Momentum in World frame //to go from the incident beam frame to the "world" frame
+ G4ThreeVector momentum_kine3_world = /*BeamToWorld */ momentum_kine3_beam ;
+ momentum_kine3_world.rotate(Beam_theta, V ); // rotation of Beam_theta on Y axis
+ momentum_kine3_world.rotate(Beam_phi, ZZ ); // rotation of Beam_phi on Z axis
// Momentum in beam frame for heavy particle
G4ThreeVector momentum_kine4_beam(sin(Theta4) * cos(phi+pi),
sin(Theta4) * sin(phi+pi),
cos(Theta4));
// Momentum in World frame
- G4ThreeVector momentum_kine4_world = BeamToWorld * momentum_kine4_beam;
-
-
+ G4ThreeVector momentum_kine4_world = /*BeamToWorld */ momentum_kine4_beam;
+ momentum_kine4_world.rotate(Beam_theta, V ); // rotation of Beam_theta on Y axis
+ momentum_kine4_world.rotate(Beam_phi, ZZ ); // rotation of Beam_phi on Z axis
+
// Emitt secondary
if(m_Reaction.GetShoot3()){
G4DynamicParticle particle3(LightName,momentum_kine3_world,Energy3);
@@ -312,6 +321,14 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep)
m_ReactionConditions->SetPhi(phi/deg);
if((phi+pi)/deg > 360 ) m_ReactionConditions->SetPhi((phi-pi)/deg);
else m_ReactionConditions->SetPhi((phi+pi)/deg);
+
+ // Angle 3 and 4 in LAB ZAXIS FRAME //
+ m_ReactionConditions->SetLABTheta(momentum_kine3_world.theta()/deg);
+ m_ReactionConditions->SetLABTheta(momentum_kine4_world.theta()/deg);
+
+ m_ReactionConditions->SetLABPhi(momentum_kine3_world.phi()/deg);
+ m_ReactionConditions->SetLABPhi(momentum_kine4_world.phi()/deg);
+
// Energy 3 and 4 //
m_ReactionConditions->SetKineticEnergy(Energy3);
m_ReactionConditions->SetKineticEnergy(Energy4);