From 8b3dbdee3f0cef06db8bd95def685b9d49b51f32 Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Fri, 3 Jul 2020 11:40:35 +0200
Subject: [PATCH] * Fixing issue in Beam reaction         - all reaction
 happenned in the first slice

---
 NPSimulation/Process/BeamReaction.cc | 13 ++++++-------
 Projects/MUGAST/16Odp17O_gs.reaction | 29 ++++++++++++++++++++++++++++
 Projects/MUGAST/Analysis.cxx         |  3 ++-
 Projects/MUGAST/ShowResults.C        |  6 ++----
 Projects/MUGAST/e793.detector        |  9 +++++----
 Projects/MUGAST/run.mac              |  3 +++
 6 files changed, 47 insertions(+), 16 deletions(-)
 create mode 100644 Projects/MUGAST/16Odp17O_gs.reaction
 create mode 100644 Projects/MUGAST/run.mac

diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 66eea8347..14f8868f6 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -141,10 +141,6 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
   m_ReactionConditions->Clear();
   m_shoot=true;
   }
-  else if (((in-m_StepSize) <= 1E-20) && m_shoot) { // last step
-    return true;
-  }
-
 
   // If the condition is met, the event is generated
   if (ratio < rand) {
@@ -153,13 +149,17 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
     if(m_ReactionType=="QFSReaction"){
         if ( m_shoot && m_QFS.IsAllowed() ) {
             return true;
-        } 
+        }
+        else
+          m_shoot=false;
     }
    
     else if(m_ReactionType=="TwoBodyReaction"){
         if ( m_shoot && m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy()) ) {
             return true;
         } 
+        else
+          m_shoot=false;
     }
   }
 
@@ -170,8 +170,6 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
     m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
     m_PreviousDirection = PrimaryTrack->GetMomentumDirection();
   }
-  else
-    return true;
 
   return false;
 }
@@ -202,6 +200,7 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
     double rand   = G4RandFlat::shoot();
     double length = (1-rand) * (m_PreviousLength);
     double reac_energy = energy + (1-rand) * (m_PreviousEnergy - energy);
+
     G4ThreeVector ldir = m_PreviousDirection;
     ldir *= length;
     localPosition = localPosition - ldir;
diff --git a/Projects/MUGAST/16Odp17O_gs.reaction b/Projects/MUGAST/16Odp17O_gs.reaction
new file mode 100644
index 000000000..2fdc85b1b
--- /dev/null
+++ b/Projects/MUGAST/16Odp17O_gs.reaction
@@ -0,0 +1,29 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%% Reaction file: 56Ni(d,p)57Ni at VAMOS-AGATA %%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+  Particle= 16O
+  Energy= 96 MeV
+  SigmaEnergy= 0 MeV
+  SigmaX= 0 mm
+  SigmaY= 0 mm
+  SigmaThetaX= 0.0 deg
+  SigmaPhiY= 0.0 deg
+  MeanThetaX= 0 mm
+  MeanPhiY= 0 mm
+  MeanX= 0 mm
+  MeanY= 0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+TwoBodyReaction
+ Beam= 16O
+ Target= 2H
+ Light= 1H
+ Heavy= 17O
+ ExcitationEnergyLight= 0.0 MeV
+ ExcitationEnergyHeavy= 0.0 MeV
+ CrossSectionPath= 16Odp17O_s12_6AMeV.dat CS
+ ShootLight= 1
+ ShootHeavy= 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/MUGAST/Analysis.cxx b/Projects/MUGAST/Analysis.cxx
index 4c86b827c..d01b32e73 100644
--- a/Projects/MUGAST/Analysis.cxx
+++ b/Projects/MUGAST/Analysis.cxx
@@ -76,9 +76,10 @@ void Analysis::Init() {
   LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",100 );
   LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",100);
   LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",100);
-  BeamCD2 = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",100);
+  BeamCD2 = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",00);
 
   FinalBeamEnergy = BeamCD2.Slow(OriginalBeamEnergy,TargetThickness*0.5,0); 
+  //FinalBeamEnergy = OriginalBeamEnergy; 
   cout << "Original beam energy: " << OriginalBeamEnergy  << " MeV      Mid-target beam energy: " << FinalBeamEnergy << "MeV " << endl; 
   reaction.SetBeamEnergy(FinalBeamEnergy);
 
diff --git a/Projects/MUGAST/ShowResults.C b/Projects/MUGAST/ShowResults.C
index a25457483..619b3a238 100644
--- a/Projects/MUGAST/ShowResults.C
+++ b/Projects/MUGAST/ShowResults.C
@@ -72,10 +72,10 @@ void ShowResults(){
   reaction->GetKinematicLine3()->Draw("c");
 
   c1->cd(2);
-  TH1F* hEx = new TH1F("hEx", "hEx",1000, -1, 5);
+  TH1F* hEx = new TH1F("hEx", "hEx",240, -1, 5);
   t->Draw("Ex>>hEx","ThetaLab>100 && ThetaLab<156","col");
   hEx->GetXaxis()->SetTitle("Ex (MeV)");
-  hEx->GetYaxis()->SetTitle("counts/100 keV");
+  hEx->GetYaxis()->SetTitle("counts/25 keV");
  
   c1->cd(3);
   TH1F *hCM = new TH1F("hCM", "hCM", 36, 0, 180); 
@@ -87,7 +87,6 @@ void ShowResults(){
     }
   }
 
-/*  
   TCanvas *c2 = new TCanvas("c2", "Control", 1000, 1000);
   c2->Divide(2,2);
   c2->cd(1);
@@ -100,7 +99,6 @@ void ShowResults(){
   t->Draw("ELab-OriginalELab>>hcE","","col");
   TLine* lE = new TLine(0,0,60,60);
   //lE->Draw();
-*/
 
 }
 
diff --git a/Projects/MUGAST/e793.detector b/Projects/MUGAST/e793.detector
index 1bb1133dc..e25775f0f 100644
--- a/Projects/MUGAST/e793.detector
+++ b/Projects/MUGAST/e793.detector
@@ -1,13 +1,14 @@
 %%%%%%%%%%Detector%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
 Target
- THICKNESS= 0.01 micrometer
+ THICKNESS= 9.43 micrometer
  ANGLE= 0 deg
  RADIUS= 10 mm
  MATERIAL= CD2
- X= 0
- Y= 0
- Z= 0
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+ NbSlices= 100
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CATSDetector
diff --git a/Projects/MUGAST/run.mac b/Projects/MUGAST/run.mac
new file mode 100644
index 000000000..02b80e5ef
--- /dev/null
+++ b/Projects/MUGAST/run.mac
@@ -0,0 +1,3 @@
+/run/beamOn 10000
+/gen/open 16Odp17O_gs.reaction
+/run/beamOn 10000
-- 
GitLab