From f1c10d8a8b28e569d5bd3e54f11fbaaf9bba4d06 Mon Sep 17 00:00:00 2001
From: fysquare <freddy.flavigny@laposte.net>
Date: Tue, 5 May 2020 17:41:36 +0200
Subject: [PATCH] Debuging QFS, heavy residue properties

---
 NPLib/Physics/NPQFS.cxx              |  3 ++-
 NPSimulation/Process/BeamReaction.cc | 18 +++++++++---------
 2 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/NPLib/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index 6063f494d..3d3850e62 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -420,6 +420,7 @@ TGraph* QFS::GetTheta2VsTheta1(double AngleStep_CM){
   for (double angle=0 ; angle <= 180 ; angle+=AngleStep_CM){
     SetThetaCM(angle*TMath::Pi()/180.);
     KineRelativistic(theta1, phi1, E1, theta2, phi2, E2);
+    Dump();
 
     vx.push_back(theta1*180./M_PI);
     vy.push_back(theta2*180./M_PI);
@@ -438,7 +439,7 @@ TGraph* QFS::GetPhi2VsPhi1(double AngleStep_CM){
   double theta1,phi1,E1,theta2,phi2,E2;
 
   for (double theta=0 ; theta <= 180 ; theta+=AngleStep_CM){
-      for (double angle=0 ; angle <= 180 ; angle+=AngleStep_CM){
+      for (double angle=-180 ; angle <= 180 ; angle+=AngleStep_CM){
           SetThetaCM(theta*TMath::Pi()/180.);
           SetPhiCM(angle*TMath::Pi()/180.);
           KineRelativistic(theta1, phi1, E1, theta2, phi2, E2);
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index b077304f5..446ffc9d6 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -463,10 +463,10 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
         /////  Momentum and angles from  kinematics  /////
         //////////////////////////////////////////////////
         // Variable where to store results
-        double Theta1, Phi1, Energy1, Theta2, Phi2, Energy2, ThetaB, PhiB, EnergyB;
+        double Theta1, Phi1, TKE1, Theta2, Phi2, TKE2, ThetaB, PhiB, TKEB;
 
         // Compute Kinematic using previously defined ThetaCM
-        m_QFS.KineRelativistic(Theta1, Phi1, Energy1, Theta2, Phi2, Energy2);
+        m_QFS.KineRelativistic(Theta1, Phi1, TKE1, Theta2, Phi2, TKE2);
 
         G4ThreeVector U(1, 0, 0);
         G4ThreeVector V(0, 1, 0);
@@ -499,7 +499,7 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
         
         G4ThreeVector momentum_kineB_beam( P_B->Px(), P_B->Py(), P_B->Pz() );
         momentum_kineB_beam = momentum_kineB_beam.unit();
-        EnergyB = m_QFS.GetEnergyImpulsionLab_B()->Energy();
+        TKEB = m_QFS.GetEnergyImpulsionLab_B()->Energy() - m_QFS.GetNucleusB()->Mass();
         G4ThreeVector momentum_kineB_world =  momentum_kineB_beam;
         momentum_kineB_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
         momentum_kineB_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
@@ -512,16 +512,16 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
   
         // Emitt secondary
         if (m_QFS.GetShoot1()) {
-            G4DynamicParticle particle1(Light1Name, momentum_kine1_world, Energy1);
+            G4DynamicParticle particle1(Light1Name, momentum_kine1_world, TKE1);
             fastStep.CreateSecondaryTrack(particle1, localPosition, time);
         }
 
         if (m_QFS.GetShoot2()) {
-            G4DynamicParticle particle2(Light2Name, momentum_kine2_world, Energy2);
+            G4DynamicParticle particle2(Light2Name, momentum_kine2_world, TKE2);
             fastStep.CreateSecondaryTrack(particle2, localPosition, time);
         }
         if (m_QFS.GetShootB()) {
-            G4DynamicParticle particleB(HeavyName, momentum_kineB_world, EnergyB);
+            G4DynamicParticle particleB(HeavyName, momentum_kineB_world, TKEB);
             fastStep.CreateSecondaryTrack(particleB, localPosition, time);
         }
 
@@ -560,9 +560,9 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
         m_ReactionConditions->SetPhi(Phi2 / deg);
         m_ReactionConditions->SetPhi(PhiB / deg);
         // Energy 1,2 and B //
-        m_ReactionConditions->SetKineticEnergy(Energy1);
-        m_ReactionConditions->SetKineticEnergy(Energy2);
-        m_ReactionConditions->SetKineticEnergy(EnergyB);
+        m_ReactionConditions->SetKineticEnergy(TKE1);
+        m_ReactionConditions->SetKineticEnergy(TKE2);
+        m_ReactionConditions->SetKineticEnergy(TKEB);
         // ThetaCM and Ex//
         m_ReactionConditions->SetThetaCM(m_QFS.GetThetaCM() / deg);
         //m_ReactionConditions->SetExcitationEnergy3(m_QFS.GetExcitation3());
-- 
GitLab