From bfa557212e17a37f982ed9835cf615a0c2325d00 Mon Sep 17 00:00:00 2001
From: De Sereville Nicolas <deserevi@ipninter.in2p3.fr>
Date: Wed, 27 Jan 2016 14:30:12 +0100
Subject: [PATCH] + update input files

+ Add consistent use of target thickness in MUGAST analysis.cx

+ Add ShowResults.C macro in MUGAST project
---
 .../MUGAST_Manu.detector                      |  1 +
 Inputs/EventGenerator/30Pdp.reaction          |  2 +-
 Projects/MUGAST/Analysis.cxx                  | 74 ++++++++++++-------
 Projects/MUGAST/ShowResults.C                 | 22 ++++++
 4 files changed, 70 insertions(+), 29 deletions(-)
 create mode 100644 Projects/MUGAST/ShowResults.C

diff --git a/Inputs/DetectorConfiguration/MUGAST_Manu.detector b/Inputs/DetectorConfiguration/MUGAST_Manu.detector
index c7f70ebf7..9135644bd 100644
--- a/Inputs/DetectorConfiguration/MUGAST_Manu.detector
+++ b/Inputs/DetectorConfiguration/MUGAST_Manu.detector
@@ -9,6 +9,7 @@ Target
 	X= 0
 	Y= 0
 	Z= 0
+	NBLAYERS= 100
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/30Pdp.reaction b/Inputs/EventGenerator/30Pdp.reaction
index 07e4253bd..fcbada0df 100644
--- a/Inputs/EventGenerator/30Pdp.reaction
+++ b/Inputs/EventGenerator/30Pdp.reaction
@@ -3,7 +3,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Beam
    Particle= 30P
-   Energy= 180
+   Energy= 360
    SigmaEnergy= 0
    SigmaX= 0
    SigmaY= 0
diff --git a/Projects/MUGAST/Analysis.cxx b/Projects/MUGAST/Analysis.cxx
index 42349c3a5..88c65f990 100644
--- a/Projects/MUGAST/Analysis.cxx
+++ b/Projects/MUGAST/Analysis.cxx
@@ -21,10 +21,11 @@
  *****************************************************************************/
 #include<iostream>
 using namespace std;
-#include"Analysis.h"
-#include"NPAnalysisFactory.h"
-#include"NPDetectorManager.h"
-#include"NPOptionManager.h"
+
+#include "Analysis.h"
+#include "NPAnalysisFactory.h"
+#include "NPDetectorManager.h"
+#include "NPOptionManager.h"
 ////////////////////////////////////////////////////////////////////////////////
 Analysis::Analysis(){
 }
@@ -33,48 +34,65 @@ Analysis::~Analysis(){
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-void Analysis::Init(){
+void Analysis::Init() {
+  // initialize input and output branches
   InitOutputBranch();
   InitInputBranch();
   
-  M2  = (TMust2Physics*) m_DetectorManager -> GetDetector("MUST2Array");
-  GD = (GaspardTracker*)  m_DetectorManager -> GetDetector("GaspardTracker");
+  // get MUST2 and Gaspard objects
+  M2 = (TMust2Physics*)  m_DetectorManager -> GetDetector("MUST2Array");
+  GD = (GaspardTracker*) m_DetectorManager -> GetDetector("GaspardTracker");
+
+  // energy losses
   LightCD2 = EnergyLoss("proton_CD2.G4table","G4Table",100 );
   LightAl = EnergyLoss("proton_Al.G4table","G4Table",100);
   LightSi = EnergyLoss("proton_Si.G4table","G4Table",100);
-  BeamCD2 = EnergyLoss("Na24[0.0]_CD2.G4table","G4Table",100);
+//  BeamCD2 = EnergyLoss("Na24[0.0]_CD2.G4table","G4Table",100);
+  BeamCD2 = EnergyLoss("P30_CD2.G4table","G4Table",100);
+
+  // get reaction information
   myReaction = new NPL::Reaction();
   myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
-   TargetThickness = 18*micrometer;
   OriginalBeamEnergy = myReaction->GetBeamEnergy();
-   Rand = TRandom3();
-   DetectorNumber = 0 ;
-   ThetaNormalTarget = 0 ;
-   ThetaM2Surface = 0;
-   Si_E_M2 = 0 ;
-   CsI_E_M2 = 0 ;
-   Energy = 0;
-   E_M2 = 0;
-  
-   ThetaGDSurface = 0;
-   X_GD = 0 ;
-   Y_GD = 0 ;
-   Z_GD = 0 ;
-   Si_E_GD = 0 ;
-   E_GD = 0;
-   Si_X_GD = 0;
-   Si_Y_GD = 0;
+
+  // target thickness
+  TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
+
+  // initialize various parameters
+  Rand = TRandom3();
+  DetectorNumber = 0;
+  ThetaNormalTarget = 0;
+  ThetaM2Surface = 0;
+  Si_E_M2 = 0;
+  CsI_E_M2 = 0;
+  Energy = 0;
+  E_M2 = 0;
+
+  ThetaGDSurface = 0;
+  X_GD = 0;
+  Y_GD = 0;
+  Z_GD = 0;
+  Si_E_GD = 0;
+  E_GD = 0;
+  Si_X_GD = 0;
+  Si_Y_GD = 0;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-void Analysis::TreatEvent(){
+void Analysis::TreatEvent() {
   // Reinitiate calculated variable
   ReInitValue();
+
+  // assume beam is centered and parallel to z axis
   double XTarget = 0;
   double YTarget = 0;
   TVector3 BeamDirection = TVector3(0,0,1);
-  double BeamEnergy = BeamCD2.Slow(OriginalBeamEnergy,Rand.Uniform(0,TargetThickness),0);
+
+  // determine beam energy for a randomized interaction point in target
+//  double BeamEnergy = BeamCD2.Slow(OriginalBeamEnergy, Rand.Uniform(0,TargetThickness), 0);
+  double BeamEnergy = BeamCD2.Slow(OriginalBeamEnergy, TargetThickness/2., 0);
   myReaction->SetBeamEnergy(BeamEnergy);
+
   //////////////////////////// LOOP on MUST2 //////////////////
   for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
     /************************************************/
diff --git a/Projects/MUGAST/ShowResults.C b/Projects/MUGAST/ShowResults.C
new file mode 100644
index 000000000..aa9a4893a
--- /dev/null
+++ b/Projects/MUGAST/ShowResults.C
@@ -0,0 +1,22 @@
+void ShowResults()
+{
+   // get tree   
+   TFile *f = new TFile("../../Outputs/Analysis/PhysicsTree.root");
+   TTree *t = (TTree*) f->Get("PhysicsTree");
+
+   // draw 
+   TCanvas *c1 = new TCanvas("c1", "kinematic information", 600, 600);
+   c1->Draw();
+   TH2F *hk = new TH2F("hk", "hk", 90, 90, 180, 200, 0, 12);
+   hk->GetXaxis()->SetTitle("#Theta_{lab} (deg)");
+   hk->GetYaxis()->SetTitle("E_{p} (MeV)");
+   t->Draw("ELab:ThetaLab>>hk");
+   // inset
+   TPad *pad = new TPad("pad1", "excitation energy", 0.45, 0.45, 0.87, 0.87);
+   pad->Draw();
+   pad->cd();
+   TH1F *hx = new TH1F("hx", "hx", 150, 5, 8);
+   hx->SetXTitle("E_{X} (MeV)");
+   hx->SetYTitle("counts / (20 keV)");
+   t->Draw("Ex>>hx");
+}
-- 
GitLab