diff --git a/NPLib/Physics/TReactionConditions.cxx b/NPLib/Physics/TReactionConditions.cxx
index bd7476763a14478031e1ad82dd5a03e3b1aab34c..8f58de2feb9fdfa68f007917b515324fe56ca2e2 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 d0632b05c6469d32e9275f234d6ca119bb86df7a..a3898a8bc599f554ca88fe7ba6eaab649ccd8a1b 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);}//!
@@ -124,6 +128,8 @@ public:
string GetParticleName (const int &i) const {return fRC_Particle_Name[i];}//!
double GetTheta (const int &i) const {return fRC_Theta[i];}//!
double GetPhi (const int &i) const {return fRC_Phi[i];}//!
+ double GetLABTheta (const int &i) const {return fRC_LAB_Theta[i];}//!
+ double GetLABPhi (const int &i) const {return fRC_LAB_Phi[i];}//!
double GetKineticEnergy (const int &i) const {return fRC_Kinetic_Energy[i];}//!
double GetMomentumDirectionX (const int &i) const {return fRC_Momentum_Direction_X[i];}//!
double GetMomentumDirectionY (const int &i) const {return fRC_Momentum_Direction_Y[i];}//!
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);