* Added (Phi, LAB_Theta, LAB_Phi) angle branches filling, and simplification...

* Added (Phi, LAB_Theta, LAB_Phi) angle branches filling, and simplification of change frame for emitted particles momentum

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();

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];}//!

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);