diff --git a/Inputs/DetectorConfiguration/NeutronWall.detector b/Inputs/DetectorConfiguration/NeutronWall.detector
index 871d2ed7997e17f7adfee49eb6440fbadc9e30bf..d930a84f0c0abb6ed848fb9d9ef34d2e3711064f 100644
--- a/Inputs/DetectorConfiguration/NeutronWall.detector
+++ b/Inputs/DetectorConfiguration/NeutronWall.detector
@@ -2,7 +2,7 @@
GeneralTarget
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Target
- THICKNESS= 10
+ THICKNESS= 10
RADIUS= 20
MATERIAL= CD2
ANGLE= 0
@@ -15,9 +15,26 @@ Target
% Y= -3.75
% Z= 566.54
% Rot= 16.61
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NeutronWall
- THETA= 20
+ THETA= 0
PHI= 0
- R= 4
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ R= 4000
+ BARS= 25
+ VETOWALL= 0
+ VWDISTANCE= 400
+ OVERLAP= 3
+ VWMATERIAL= BC400
+ NWMATERIAL= NE213
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+NeutronWall
+ THETA= 0
+ PHI= 0
+ R= 6000
+ BARS= 25
+ VETOWALL= 0
+ VWDISTANCE= 1000
+ OVERLAP= 3
+ VWMATERIAL= BC400
+ NWMATERIAL= NE213
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/3He.source b/Inputs/EventGenerator/3He.source
index b80e9a4a1fec2107cdf38b2507847f66ca9022fa..66efecc3e4be50185b581bdd2afda851991e014b 100644
--- a/Inputs/EventGenerator/3He.source
+++ b/Inputs/EventGenerator/3He.source
@@ -5,9 +5,9 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 710
+ EnergyHigh= 525
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 10
+ HalfOpenAngleMax= 90
x0= 0
y0= 0
z0= 0
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index 9c9bc63fd260253abb84404ccf4e4b31a487ef90..1c9997fde237f339d97cecb394ca6d5542eea6a5 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -5,9 +5,9 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 795
+ EnergyHigh= 590
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 10
+ HalfOpenAngleMax= 90
x0= 0
y0= 0
z0= 0
diff --git a/Inputs/EventGenerator/deuton.source b/Inputs/EventGenerator/deuton.source
index 7a864b0fb3018a654583b78a9f4ab1e2d1e00a19..4bd428cd0b4f89d52538700126a6914038b11e9a 100644
--- a/Inputs/EventGenerator/deuton.source
+++ b/Inputs/EventGenerator/deuton.source
@@ -5,9 +5,9 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 265
+ EnergyHigh= 205
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 10
+ HalfOpenAngleMax= 90
x0= 0
y0= 0
z0= 0
diff --git a/Inputs/EventGenerator/proton.source b/Inputs/EventGenerator/proton.source
index 5e33cde9ba90f3b2077da3ec30448bfe620cace0..75ce81c443c5a62b35a19132ac185721b8412dde 100644
--- a/Inputs/EventGenerator/proton.source
+++ b/Inputs/EventGenerator/proton.source
@@ -5,7 +5,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 195
+ EnergyHigh= 250
HalfOpenAngleMin= 0
HalfOpenAngleMax= 90
x0= 0
diff --git a/Inputs/EventGenerator/triton.source b/Inputs/EventGenerator/triton.source
index f67b1d5fb4d199a7907c33275ed56469da4ba73c..28ecc223289d59d56ff668874ffa5543fe304440 100644
--- a/Inputs/EventGenerator/triton.source
+++ b/Inputs/EventGenerator/triton.source
@@ -5,9 +5,9 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 315
+ EnergyHigh= 235
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 10
+ HalfOpenAngleMax= 90
x0= 0
y0= 0
z0= 0
diff --git a/NPLib/Detectors/NeutronWall/TNeutronWallData.h b/NPLib/Detectors/NeutronWall/TNeutronWallData.h
index caee4ad88eb9ddcf46f11fd553455fea86e27778..f36a9733a03fe23dbcf8ea88e46ded95ff15dded 100644
--- a/NPLib/Detectors/NeutronWall/TNeutronWallData.h
+++ b/NPLib/Detectors/NeutronWall/TNeutronWallData.h
@@ -36,12 +36,24 @@ class TNeutronWallData : public TObject {
private:
// Energy
vector<UShort_t> fNeutronWall_E_DetectorNbr;
+ vector<UShort_t> fNeutronWall_E_PadNbr;
vector<Double_t> fNeutronWall_Energy;
// Time
vector<UShort_t> fNeutronWall_T_DetectorNbr;
+ vector<UShort_t> fNeutronWall_T_PadNbr;
vector<Double_t> fNeutronWall_Time;
+ // Energy
+ vector<UShort_t> fVetoWall_E_DetectorNbr;
+ vector<UShort_t> fVetoWall_E_PadNbr;
+ vector<Double_t> fVetoWall_Energy;
+
+ // Time
+ vector<UShort_t> fVetoWall_T_DetectorNbr;
+ vector<UShort_t> fVetoWall_T_PadNbr;
+ vector<Double_t> fVetoWall_Time;
+
//////////////////////////////////////////////////////////////
// Constructor and destructor
@@ -64,36 +76,73 @@ class TNeutronWallData : public TObject {
// frequently used methods
// add //! to avoid ROOT creating dictionnary for the methods
public:
- ////////////////////// SETTERS ////////////////////////
+ ////////////////////// SETTERS NEUTRON WALL ////////////////////////
// Energy
inline void SetE_DetectorNbr(const UShort_t& DetNbr)
{fNeutronWall_E_DetectorNbr.push_back(DetNbr);} //!
+ inline void SetE_PadNbr(const UShort_t& PadNbr)
+ {fNeutronWall_E_PadNbr.push_back(PadNbr);} //!
inline void Set_Energy(const Double_t& Energy)
{fNeutronWall_Energy.push_back(Energy);}//!
// Prefer global setter so that all vectors have the same size
- inline void SetEnergy(const UShort_t& DetNbr,const Double_t& Energy) {
+ inline void SetEnergy(const UShort_t& DetNbr, const UShort_t& PadNbr, const Double_t& Energy) {
SetE_DetectorNbr(DetNbr);
+ SetE_PadNbr(PadNbr);
Set_Energy(Energy);
};//!
// Time
inline void SetT_DetectorNbr(const UShort_t& DetNbr)
{fNeutronWall_T_DetectorNbr.push_back(DetNbr);} //!
+ inline void SetT_PadNbr(const UShort_t& PadNbr)
+ {fNeutronWall_T_PadNbr.push_back(PadNbr);} //!
inline void Set_Time(const Double_t& Time)
{fNeutronWall_Time.push_back(Time);}//!
// Prefer global setter so that all vectors have the same size
- inline void SetTime(const UShort_t& DetNbr,const Double_t& Time) {
+ inline void SetTime(const UShort_t& DetNbr,const UShort_t& PadNbr,const Double_t& Time) {
SetT_DetectorNbr(DetNbr);
+ SetT_PadNbr(PadNbr);
Set_Time(Time);
};//!
+ ////////////////////// SETTERS VETO WALL ////////////////////////
+ // Energy
+ inline void SetE_VetoDetectorNbr(const UShort_t& DetNbr)
+ {fVetoWall_E_DetectorNbr.push_back(DetNbr);} //!
+ inline void SetE_VetoPadNbr(const UShort_t& PadNbr)
+ {fVetoWall_E_PadNbr.push_back(PadNbr);} //!
+ inline void Set_VetoEnergy(const Double_t& Energy)
+ {fVetoWall_Energy.push_back(Energy);}//!
+ // Prefer global setter so that all vectors have the same size
+ inline void SetVetoEnergy(const UShort_t& DetNbr,const UShort_t& PadNbr,const Double_t& Energy) {
+ SetE_VetoDetectorNbr(DetNbr);
+ SetE_VetoPadNbr(PadNbr);
+ Set_VetoEnergy(Energy);
+ };//!
+
+ // Time
+ inline void SetT_VetoDetectorNbr(const UShort_t& DetNbr)
+ {fVetoWall_T_DetectorNbr.push_back(DetNbr);} //!
+ inline void SetT_VetoPadNbr(const UShort_t& PadNbr)
+ {fVetoWall_T_PadNbr.push_back(PadNbr);} //!
+ inline void Set_VetoTime(const Double_t& Time)
+ {fVetoWall_Time.push_back(Time);}//!
+ // Prefer global setter so that all vectors have the same size
+ inline void SetVetoTime(const UShort_t& DetNbr,const UShort_t& PadNbr,const Double_t& Time) {
+ SetT_VetoDetectorNbr(DetNbr);
+ SetT_VetoPadNbr(PadNbr);
+ Set_VetoTime(Time);
+ };//!
+
- ////////////////////// GETTERS ////////////////////////
+ ////////////////////// GETTERS NEUTRON WALL ////////////////////////
// Energy
inline UShort_t GetMultEnergy() const
{return fNeutronWall_E_DetectorNbr.size();}
inline UShort_t GetE_DetectorNbr(const unsigned int &i) const
{return fNeutronWall_E_DetectorNbr[i];}//!
+ inline UShort_t GetE_PadNbr(const unsigned int &i) const
+ {return fNeutronWall_E_PadNbr[i];}//!
inline Double_t Get_Energy(const unsigned int &i) const
{return fNeutronWall_Energy[i];}//!
@@ -102,9 +151,30 @@ class TNeutronWallData : public TObject {
{return fNeutronWall_T_DetectorNbr.size();}
inline UShort_t GetT_DetectorNbr(const unsigned int &i) const
{return fNeutronWall_T_DetectorNbr[i];}//!
+ inline UShort_t GetT_PadNbr(const unsigned int &i) const
+ {return fNeutronWall_T_PadNbr[i];}//!
inline Double_t Get_Time(const unsigned int &i) const
{return fNeutronWall_Time[i];}//!
+ ////////////////////// GETTERS VETO WALL ////////////////////////
+ // Energy
+ inline UShort_t GetVetoMultEnergy() const
+ {return fVetoWall_E_DetectorNbr.size();}
+ inline UShort_t GetE_VetoDetectorNbr(const unsigned int &i) const
+ {return fVetoWall_E_DetectorNbr[i];}//!
+ inline UShort_t GetE_VetoPadNbr(const unsigned int &i) const
+ {return fVetoWall_E_PadNbr[i];}//!
+ inline Double_t Get_VetoEnergy(const unsigned int &i) const
+ {return fVetoWall_Energy[i];}//!
+
+ // Time
+ inline UShort_t GetVetoMultTime() const
+ {return fVetoWall_T_DetectorNbr.size();}
+ inline UShort_t GetT_VetoDetectorNbr(const unsigned int &i) const
+ {return fVetoWall_T_DetectorNbr[i];}//!
+ inline Double_t Get_VetoTime(const unsigned int &i) const
+ {return fVetoWall_Time[i];}//!
+
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
diff --git a/NPLib/Detectors/NeutronWall/TNeutronWallPhysics.cxx b/NPLib/Detectors/NeutronWall/TNeutronWallPhysics.cxx
index 5cb68a7cdb5c5a2d22b49e2cf61b1e2de723ffce..ddd29285e666c7ac50ce7848cdd0d6c0aa093dd5 100644
--- a/NPLib/Detectors/NeutronWall/TNeutronWallPhysics.cxx
+++ b/NPLib/Detectors/NeutronWall/TNeutronWallPhysics.cxx
@@ -71,6 +71,7 @@ void TNeutronWallPhysics::BuildPhysicalEvent() {
for (UShort_t t = 0; t < m_PreTreatedData->GetMultTime(); t++) {
if (m_PreTreatedData->GetE_DetectorNbr(e) == m_PreTreatedData->GetT_DetectorNbr(t)) {
DetectorNumber.push_back(m_PreTreatedData->GetE_DetectorNbr(e));
+ PadNumber.push_back(m_PreTreatedData->GetE_PadNbr(e));
Energy.push_back(m_PreTreatedData->Get_Energy(e));
Time.push_back(m_PreTreatedData->Get_Time(t));
}
@@ -94,7 +95,7 @@ void TNeutronWallPhysics::PreTreat() {
if (m_EventData->Get_Energy(i) > m_E_RAW_Threshold) {
Double_t Energy = Cal->ApplyCalibration("NeutronWall/ENERGY"+NPL::itoa(m_EventData->GetE_DetectorNbr(i)),m_EventData->Get_Energy(i));
if (Energy > m_E_Threshold) {
- m_PreTreatedData->SetEnergy(m_EventData->GetE_DetectorNbr(i), Energy);
+ m_PreTreatedData->SetEnergy(m_EventData->GetE_DetectorNbr(i),m_EventData->GetE_PadNbr(i), Energy);
}
}
}
@@ -102,7 +103,7 @@ void TNeutronWallPhysics::PreTreat() {
// Time
for (UShort_t i = 0; i < m_EventData->GetMultTime(); ++i) {
Double_t Time= Cal->ApplyCalibration("NeutronWall/TIME"+NPL::itoa(m_EventData->GetT_DetectorNbr(i)),m_EventData->Get_Time(i));
- m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i), Time);
+ m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i),m_EventData->GetT_PadNbr(i), Time);
}
}
@@ -175,6 +176,7 @@ void TNeutronWallPhysics::ReadAnalysisConfig() {
///////////////////////////////////////////////////////////////////////////
void TNeutronWallPhysics::Clear() {
DetectorNumber.clear();
+ PadNumber.clear();
Energy.clear();
Time.clear();
}
diff --git a/NPLib/Detectors/NeutronWall/TNeutronWallPhysics.h b/NPLib/Detectors/NeutronWall/TNeutronWallPhysics.h
index 340646e0c3fd7e9787b71162bf3a200a106a6ffd..c837d531c15a066c6b1b0f578fed318196d15489 100644
--- a/NPLib/Detectors/NeutronWall/TNeutronWallPhysics.h
+++ b/NPLib/Detectors/NeutronWall/TNeutronWallPhysics.h
@@ -64,6 +64,7 @@ class TNeutronWallPhysics : public TObject, public NPL::VDetector {
// output ROOT file
public:
vector<int> DetectorNumber;
+ vector<int> PadNumber;
vector<double> Energy;
vector<double> Time;
diff --git a/NPSimulation/Detectors/NeutronWall/NeutronWall.cc b/NPSimulation/Detectors/NeutronWall/NeutronWall.cc
index eef734ea40e5e2149e92feb9a75a39d26b58af51..20456b6ff7132d38ea561a02f145470f4e31f9fe 100644
--- a/NPSimulation/Detectors/NeutronWall/NeutronWall.cc
+++ b/NPSimulation/Detectors/NeutronWall/NeutronWall.cc
@@ -1,18 +1,18 @@
/*****************************************************************************
- * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
*****************************************************************************/
/*****************************************************************************
- * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu *
+ * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu *
* *
- * Creation Date : June 2016 *
+ * Creation Date : June 2016 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * This class describe NeutronWall simulation *
+ * This class describe NeutronWall simulation *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -23,10 +23,14 @@
#include <sstream>
#include <cmath>
#include <limits>
+#include <cstring>
+#include <string>
+
//G4 Geometry object
#include "G4Tubs.hh"
#include "G4Box.hh"
#include "G4Trd.hh"
+#include "G4SubtractionSolid.hh"
//G4 sensitive
#include "G4SDManager.hh"
@@ -58,20 +62,49 @@ namespace NeutronWall_NS{
const double EnergyThreshold = 0.1*MeV;
const double ResoTime = 4.5*ns ;
const double ResoEnergy = 5.0*MeV ;
- const double NS_X = 1.1*m;
- const double NS_Y = 1.1*m;
- const double NS_Z = 5*cm;
- const double Al_X = 1.1*m;
- const double Al_Y = 1.1*m;
- const double Al_Z = 0.2*m;
- const double NE213_X = 1.0*m;
- const double NE213_Y = 3.66*cm;
- const double NE213_Z = 3.025*cm;
- const double Py_X = 1.006*m;
- const double Py_Y = 3.81*cm;
- const double Py_Z = 3.175*cm;
-
- const string Scintillator = "NE213";
+ //The size of NS should depend on the distance between NeutronWall and plastic Bar right now
+ double NS_X = 2020.0*mm;
+ double NS_Y = 2020.0*mm;
+ //the front and back aluminum sheet are both 0.8 thick whereas 143.5 is user assumed heigh in z
+ double NS_Z = (143.5+0.8+0.8)*mm;
+ //using Alouter minus Alinner, one get an Al frame
+ const double Alinner_X = 2000.0*mm;
+ const double Alinner_Y = 2000.0*mm;
+ const double Alinner_Z = 143.5*mm;
+ const double Alouter_X = 2020.0*mm;
+ const double Alouter_Y = 2020.0*mm;
+ const double Alouter_Z = (143.5+0.8+0.8)*mm;
+
+ const double frame_thickness = 10*mm;
+
+ const double Scintillator_X = 1994.0*mm;
+ const double Scintillator_Y = 70.2*mm;
+ const double Scintillator_Z = 57.5*mm;
+ //do the same thing to Pyrex tube to create a volume to store Scintillator
+
+ const double Py_Xinner = 1994.0*mm;
+ const double Py_Yinner = 70.2*mm;
+ const double Py_Zinner = 57.5*mm;
+ const double Py_Xouter = 2000.0*mm;
+ const double Py_Youter = 76.2*mm;
+ const double Py_Zouter = 63.5*mm;
+
+ const double seperation_between_pyrex = 3.0*mm;
+ const double upper_gap = 10.0*mm;
+
+ //Add elements about the plastic bars
+ double PlasticBar_X = 0.0*mm;
+ double PlasticBar_Y = 0.0*mm;
+ double PlasticBar_Z = 10.0*mm;
+
+ //Add total height of neutronwall and vetowall for comparision
+ double TotalHeightOfNeutronWall = 0.0*mm;
+ double TotalHeightOfVetoWall = 0.0*mm;
+
+ //Scale down factor
+ double ScaleDownFactor = 0;
+
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -80,27 +113,30 @@ namespace NeutronWall_NS{
NeutronWall::NeutronWall(){
m_Event = new TNeutronWallData() ;
m_NeutronWallScorer = 0;
+ m_VetoWallScorer = 0;
m_NeutronWall_log = 0;
m_AlCase_log = 0;
m_Quartz_log = 0;
m_QuartzCap_log = 0;
m_PMTube_log = 0;
- m_NE213_log = 0;
+ m_Scintillator_log = 0;
m_ShadowBar_log = 0;
+ m_PlasticBar_log = 0;
// RGB Color + Transparency
- m_VisNE213 = new G4VisAttributes(G4Colour(0, 1, 0, 0.5));
- m_VisQuartz = new G4VisAttributes(G4Colour(0, 1, 1, 0.5));
- m_VisAl = new G4VisAttributes(G4Colour(0.5,0.5,0.5,0.5));
- m_VisNW = new G4VisAttributes(G4Colour(0.4,0.2,0.8,0.5));
+ m_VisScintillator = new G4VisAttributes(G4Colour(1, 0.843137, 0, 0.3)); //gold
+ m_VisQuartz = new G4VisAttributes(G4Colour(0,1, 0, 0.1)); //green
+ m_VisAl = new G4VisAttributes(G4Colour(173.0/255.0,178.0/255.0,189.0/255.0,0.3)); //Al
+ m_VisNW = new G4VisAttributes(G4Colour(0.972549,0.972549,1,0.1)); //ghostwhite
+ m_VisPlasticBar = new G4VisAttributes(G4Colour(0.9,0,0.9,1.0)); //pink
}
NeutronWall::~NeutronWall(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void NeutronWall::AddNeutronWall(double R, double Theta, double Phi, double X, double Y, double Z, double Rotation){
+void NeutronWall::AddNeutronWall(double R, double Theta, double Phi, double X, double Y, double Z, double Rotation, int Bars, string NWMaterial, double VWDistance, int VetoWall, string VWMaterial, double Overlap){
m_R.push_back(R);
m_Theta.push_back(Theta);
m_Phi.push_back(Phi);
@@ -108,12 +144,18 @@ void NeutronWall::AddNeutronWall(double R, double Theta, double Phi, double X
m_Y.push_back(Y);
m_Z.push_back(Z);
m_Rot.push_back(Rotation);
+ m_Bars.push_back(Bars);
+ m_NWMaterial.push_back(NWMaterial);
+ m_VWDistance.push_back(VWDistance);
+ m_VetoWall.push_back(VetoWall);
+ m_VWMaterial.push_back(VWMaterial);
+ m_Overlap.push_back(Overlap);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void NeutronWall::BuildDetector(){
-
-
+
+
}
@@ -134,6 +176,12 @@ void NeutronWall::ReadConfiguration(string Path){
double Theta = 0 , Phi = 0 , R = 0 ;
double X = 0 , Y = 0 , Z = 0 ;
double Rot =0;
+ int Bars = 0;
+ string NWMaterial = "NE213";
+ double VWDistance = 0.0;
+ int VetoWall = 0;
+ string VWMaterial = "BC400";
+ double Overlap = 3;
bool check_Theta = false ;
bool check_Phi = false ;
@@ -143,7 +191,11 @@ void NeutronWall::ReadConfiguration(string Path){
bool check_Y = false ;
bool check_Z = false ;
bool ReadingStatus = false ;
-
+ bool check_Bars = false ;
+ bool check_NWMaterial = false ;
+ //bool check_VWDistance = false ;
+ //bool check_VetoWall = false ;
+ //bool check_VWMaterial = false ;
while (!ConfigFile.eof()) {
getline(ConfigFile, LineBuffer);
@@ -197,8 +249,8 @@ void NeutronWall::ReadConfiguration(string Path){
check_R = true;
ConfigFile >> DataBuffer ;
R = atof(DataBuffer.c_str()) ;
- R = R * m;
- G4cout << "R: " << R/m << G4endl;
+ R = R * mm;
+ G4cout << "R: " << R/mm << G4endl;
}
//Position method
@@ -236,6 +288,56 @@ void NeutronWall::ReadConfiguration(string Path){
G4cout << "Rotation: " << Rot/deg << G4endl;
}
+ //Bar number
+ else if (DataBuffer.compare(0, 5, "BARS=") == 0){
+ check_Bars = true;
+ ConfigFile >> DataBuffer ;
+ Bars = atoi(DataBuffer.c_str()) ;
+ G4cout << "Bars: " << Bars << G4endl;
+ }
+
+
+ //Material type
+ else if (DataBuffer.compare(0, 11, "NWMATERIAL=") == 0){
+ check_NWMaterial = true;
+ ConfigFile >> DataBuffer ;
+ NWMaterial = DataBuffer;
+ G4cout << "NWMaterials: " << NWMaterial << G4endl;
+ }
+
+ //Distance
+ else if (DataBuffer.compare(0, 11, "VWDISTANCE=") == 0){
+ //check_VWDistance = true;
+ ConfigFile >> DataBuffer ;
+ VWDistance = atof(DataBuffer.c_str());
+ VWDistance = VWDistance * mm;
+ G4cout << "VWDistance: " << VWDistance << G4endl;
+ }
+
+ //Decide whether to add the vetowall or not, 1 means yes, 0 means no
+ else if (DataBuffer.compare(0, 9, "VETOWALL=") == 0){
+ //check_VetoWall = true;
+ ConfigFile >> DataBuffer ;
+ VetoWall = atoi(DataBuffer.c_str());
+ G4cout << "VetoWall: " << VetoWall << G4endl;
+ }
+
+ //VetoWall Material
+ else if (DataBuffer.compare(0, 11, "VWMATERIAL=") == 0){
+ //check_VWMaterial = true;
+ ConfigFile >> DataBuffer ;
+ VWMaterial = DataBuffer ;
+ G4cout << "VWMaterial: " << VWMaterial << G4endl;
+ }
+
+ //Overlap
+ else if (DataBuffer.compare(0, 8, "OVERLAP=") == 0){
+ ConfigFile >> DataBuffer ;
+ Overlap = atof(DataBuffer.c_str());
+ Overlap = Overlap*mm;
+ G4cout << "Overlap: " << Overlap << G4endl;
+ }
+
///////////////////////////////////////////////////
// If no Detector Token and no comment, toggle out
else{
@@ -246,9 +348,9 @@ void NeutronWall::ReadConfiguration(string Path){
/////////////////////////////////////////////////
// If All necessary information there, toggle out
- if (( check_Theta && check_Phi && check_R)
+ if (( check_Theta && check_Phi && check_R && check_Bars && check_NWMaterial)
||
- ( check_X && check_Y && check_Z && check_rotation)){
+ ( check_X && check_Y && check_Z && check_rotation && check_Bars && check_NWMaterial)){
// Convert Cartesian to Spherical (detector always face the target)
@@ -258,7 +360,7 @@ void NeutronWall::ReadConfiguration(string Path){
Phi = atan2(Y,X);
}
- AddNeutronWall(R,Theta,Phi,X,Y,Z,Rot);
+ AddNeutronWall(R,Theta,Phi,X,Y,Z,Rot, Bars, NWMaterial, VWDistance, VetoWall, VWMaterial, Overlap);
// Reinitialisation of Check Boolean
check_Theta = false ;
@@ -269,6 +371,11 @@ void NeutronWall::ReadConfiguration(string Path){
check_Y = false ;
check_Z = false ;
ReadingStatus = false ;
+ check_Bars = false ;
+ check_NWMaterial = false ;
+ //check_VWDistance = false ;
+ //check_VetoWall = false ;
+ //check_VWMaterial = false ;
G4cout << "///"<< G4endl ;
}
}
@@ -295,13 +402,27 @@ void NeutronWall::ConstructDetector(G4LogicalVolume* world){
G4ThreeVector Y(ii,jj,kk);
G4ThreeVector w = Det_pos.unit();
G4ThreeVector u = w.cross(Y);
- G4ThreeVector v = w.cross(u);
+ G4ThreeVector v = u.cross(w);
v = v.unit();
u = u.unit();
- G4RotationMatrix* Rot = new G4RotationMatrix(u,v,w);
- G4Material* ScintMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(NeutronWall_NS::Scintillator);
+
+ //Initialize scale down factor , measured from face to face but VWDistance is from center to center
+ NeutronWall_NS::ScaleDownFactor = ((m_R[i]+NeutronWall_NS::NS_Z*0.5-(m_VWDistance[i]-NeutronWall_NS::NS_Z*0.5+0.5*NeutronWall_NS::PlasticBar_Z))/(m_R[i]+NeutronWall_NS::NS_Z*0.5));
+
+ if (m_VetoWall[i] == 1){
+ //Initialize property about the plastic bar, if exists
+ NeutronWall_NS::PlasticBar_X = NeutronWall_NS::Py_Xouter;
+ NeutronWall_NS::PlasticBar_Y = NeutronWall_NS::Py_Youter;
+
+ // If VetoWall exists, then extend the space of NeutronWall_NS add 3mm to Z in order to house the 1mm seperation between front and back layer of the veto wall
+ NeutronWall_NS::NS_Z = 2.0*(m_VWDistance[i]+1.5*NeutronWall_NS::PlasticBar_Z+3.0*mm);
+ }
+
+ G4RotationMatrix* Rot = new G4RotationMatrix(v,u,w);
+
+ G4Material* ScintMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(m_NWMaterial[i]);
G4Material* vacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
G4Material* Aluminum = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
G4Material* Pyrex = MaterialManager::getInstance()->GetMaterialFromLibrary("Pyrex");
@@ -312,57 +433,90 @@ void NeutronWall::ConstructDetector(G4LogicalVolume* world){
m_NeutronWall_log = new G4LogicalVolume(NeutronWall_box,vacuum,"NeutronWall_Log",0,0,0);
m_NeutronWall_log->SetVisAttributes(m_VisNW);
- //Aluminum case
- G4Box* AlCase_box = new G4Box("AlCase_Box",NeutronWall_NS::Al_X*0.5,NeutronWall_NS::Al_Y*0.5,NeutronWall_NS::Al_Z*0.5);
- m_AlCase_log = new G4LogicalVolume(AlCase_box,Aluminum,"AlCase_Log",0,0,0);
+ //Aluminum inner box (subtractee)
+ G4Box* Alinner_box = new G4Box("Alinner_box",NeutronWall_NS::Alinner_X*0.5,NeutronWall_NS::Alinner_Y*0.5,NeutronWall_NS::Alinner_Z*0.5);
+ //Aluminum outer box (subtracter)
+ G4Box* Alouter_box = new G4Box("Alouter_box",NeutronWall_NS::Alouter_X*0.5,NeutronWall_NS::Alouter_Y*0.5,NeutronWall_NS::Alouter_Z*0.5);
+
+ G4SubtractionSolid* AlCase_frame = new G4SubtractionSolid("AlCase_frame", Alouter_box, Alinner_box);
+
+ m_AlCase_log = new G4LogicalVolume(AlCase_frame,Aluminum,"AlCase_Log",0,0,0);
m_AlCase_log->SetVisAttributes(m_VisAl);
- //Quartz tube volume
- G4Box* Quartz_box = new G4Box("Quartz_Box",NeutronWall_NS::Py_X*0.5,NeutronWall_NS::Py_Y*0.5,NeutronWall_NS::Py_Z*0.5);
+ //Quartz tube
+ G4Box* Quartz_boxinner = new G4Box("Quartz_Boxinner",NeutronWall_NS::Py_Xinner*0.5,NeutronWall_NS::Py_Yinner*0.5,NeutronWall_NS::Py_Zinner*0.5);
+ G4Box* Quartz_boxouter = new G4Box("Quartz_Box",NeutronWall_NS::Py_Xouter*0.5,NeutronWall_NS::Py_Youter*0.5,NeutronWall_NS::Py_Zouter*0.5);
+ G4SubtractionSolid* Quartz_box = new G4SubtractionSolid("Quartz_box", Quartz_boxouter, Quartz_boxinner);
+
m_Quartz_log = new G4LogicalVolume(Quartz_box,Pyrex,"Quartz_Log",0,0,0);
m_Quartz_log->SetVisAttributes(m_VisQuartz);
+ //??currently unused
G4Tubs* QuartzCap_tube = new G4Tubs("QuartsCap_Tube",0,7.99*cm,0.3175*cm,0.0*deg,360.0*deg);
m_QuartzCap_log = new G4LogicalVolume(QuartzCap_tube,Aluminum,"QuartzCap_Log",0,0,0);
- //NE213 Volume
- G4Box* NE213_box = new G4Box("NE213_Box",NeutronWall_NS::NE213_X*0.5,NeutronWall_NS::NE213_Y*0.5,NeutronWall_NS::NE213_Z*0.5);
- m_NE213_log = new G4LogicalVolume(NE213_box,ScintMaterial,"NE213_Log",0,0,0);
+ //Scintillator
+ G4Box* Scintillator_box = new G4Box("Scintillator_Box",NeutronWall_NS::Scintillator_X*0.5,NeutronWall_NS::Scintillator_Y*0.5,NeutronWall_NS::Scintillator_Z*0.5);
+ m_Scintillator_log = new G4LogicalVolume(Scintillator_box,ScintMaterial,"Scintillator_Log",0,0,0);
- m_NE213_log->SetVisAttributes(m_VisNE213);
- m_NE213_log->SetSensitiveDetector(m_NeutronWallScorer);
+ m_Scintillator_log->SetVisAttributes(m_VisScintillator);
+ m_Scintillator_log->SetSensitiveDetector(m_NeutronWallScorer);
- //Shadow bar construction
+ //Shadow bar construction??currently unused
G4Trd* ShadowBar_trd = new G4Trd("ShadowBar_Trd",6.51/2*cm, 7.23/2*cm, 6.79/2*cm, 7.66/2*cm, 30./2*cm);
m_ShadowBar_log = new G4LogicalVolume(ShadowBar_trd, Aluminum, "ShadowBar_Log");
- //******************* Placement *******************//
- m_NE213Tube_phys = new G4PVPlacement(0,G4ThreeVector(0,0,0),
- m_NE213_log,
- "NE213Tube_phys",
- m_Quartz_log,
- false,
- 0);
- for (int j = 0; j < 12; j++ ) {
- G4double detector_y = (j-6)*8.60*cm;
- m_Quartz_phys = new G4PVPlacement(0,G4ThreeVector(0,detector_y,0),m_Quartz_log,
- "Quartz_phys",m_NeutronWall_log,false,j);
+ if (m_VetoWall[i] == 1){
+ //Initialize total height
+ NeutronWall_NS::TotalHeightOfNeutronWall = m_Bars[i]*NeutronWall_NS::Py_Youter + (m_Bars[i]-1)*NeutronWall_NS::seperation_between_pyrex;
+ NeutronWall_NS::TotalHeightOfVetoWall = (m_Bars[i]+1)/2*NeutronWall_NS::PlasticBar_Y+(m_Bars[i]-1)/2*(NeutronWall_NS::PlasticBar_Y-2*m_Overlap[i]);
+ if (NeutronWall_NS::TotalHeightOfNeutronWall*NeutronWall_NS::ScaleDownFactor > NeutronWall_NS::TotalHeightOfVetoWall){
+ G4cout << "\t*************************************************************" <<endl;
+ G4cout << "\t* The shadow of VetoWall is not enough to cover NeutronWall *" <<endl;
+ G4cout << "\t* Re-input VWDistance and Overlap *" <<endl;
+ G4cout << "\t*************************************************************" <<endl;
+ exit(1);
+ }
+ //PlasticBar
+ G4Material* Plastic = MaterialManager::getInstance()->GetMaterialFromLibrary(m_VWMaterial[i]);
+ G4Box* PlasticBar_box = new G4Box("PlasticBar_Box", NeutronWall_NS::PlasticBar_X*0.5, NeutronWall_NS::PlasticBar_Y*0.5, NeutronWall_NS::PlasticBar_Z*0.5);
+ m_PlasticBar_log = new G4LogicalVolume(PlasticBar_box, Plastic, "PlasticBar_Log");
+ m_PlasticBar_log->SetSensitiveDetector(m_VetoWallScorer);
+ m_PlasticBar_log->SetVisAttributes(m_VisPlasticBar);
}
+
+ //******************* Placement *******************//
//----Aluminum Case----
- m_AlCase_phys = new G4PVPlacement(0,G4ThreeVector(0,0, 4.2*cm),m_AlCase_log,
+ m_AlCase_phys = new G4PVPlacement(0,G4ThreeVector(0,0,0),m_AlCase_log,
"AlCase_phys",m_NeutronWall_log,false,0);
+ for (int j = 0; j < 25; j++ ) {
+ m_ScintillatorTube_phys = new G4PVPlacement(0,G4ThreeVector(0,(NeutronWall_NS::NS_Y*0.5-NeutronWall_NS::frame_thickness-NeutronWall_NS::upper_gap-NeutronWall_NS::Py_Youter*0.5-j*(NeutronWall_NS::Py_Youter+NeutronWall_NS::seperation_between_pyrex)),0),m_Scintillator_log,"ScintillatorTube_phys",m_NeutronWall_log,false,j);
+ m_Quartz_phys = new G4PVPlacement(0,G4ThreeVector(0*mm,(NeutronWall_NS::NS_Y*0.5-NeutronWall_NS::frame_thickness-NeutronWall_NS::upper_gap-NeutronWall_NS::Py_Youter*0.5-j*(NeutronWall_NS::Py_Youter+NeutronWall_NS::seperation_between_pyrex)),0*mm),m_Quartz_log,
+ "Quartz_phys",m_NeutronWall_log,false,j);
+ }
+ for(int j = 0; j<25; j++){
+ if (m_VetoWall[i] == 1){
+ G4ThreeVector VetoTransB = G4ThreeVector(0,0.5*NeutronWall_NS::TotalHeightOfVetoWall-0.5*NeutronWall_NS::PlasticBar_Y-(j/2)*(NeutronWall_NS::PlasticBar_Y+NeutronWall_NS::PlasticBar_Y-2*m_Overlap[i]),-m_VWDistance[i]);//+Det_pos;
+ G4ThreeVector VetoTransF = G4ThreeVector(0,0.5*NeutronWall_NS::TotalHeightOfVetoWall-1.5*NeutronWall_NS::PlasticBar_Y+m_Overlap[i]-(j/2)*(NeutronWall_NS::PlasticBar_Y+NeutronWall_NS::PlasticBar_Y-2*m_Overlap[i]),-m_VWDistance[i]-1*mm-NeutronWall_NS::PlasticBar_Z);//+Det_pos;
+ if (j%2 == 0){
+ m_PlasticBar_phys = new G4PVPlacement(0,VetoTransB,m_PlasticBar_log,"PlasticBar_phys",m_NeutronWall_log,false,j);
+ }
+ else {
+ m_PlasticBar_phys = new G4PVPlacement(0, VetoTransF,m_PlasticBar_log,"PlasticBar_phys",m_NeutronWall_log,false,j);
+
+ }
+ }
+ }
/****************** Place the walls*************************/
m_NeutronWall_phys = new G4PVPlacement(G4Transform3D(*Rot, Det_pos),
- m_NeutronWall_log,
- "NeutronWall_phys",world,false,i);
-
-
+ m_NeutronWall_log,
+ "NeutronWall_phys",world,false,i);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -374,7 +528,7 @@ void NeutronWall::InitializeRootOutput(){
pTree->Branch("NeutronWall", "TNeutronWallData", &m_Event) ;
pTree->SetBranchAddress("NeutronWall", &m_Event) ;
}
-
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
@@ -398,18 +552,54 @@ void NeutronWall::ReadSensitive(const G4Event* event){
double Energy = RandGauss::shoot(Info[0],NeutronWall_NS::ResoEnergy);
if(Energy>NeutronWall_NS::EnergyThreshold){
double Time = RandGauss::shoot(Info[1],NeutronWall_NS::ResoTime);
- int DetectorNbr = (int) Info[2];
- m_Event->SetEnergy(DetectorNbr,Energy);
- m_Event->SetTime(DetectorNbr,Time);
+ int DetectorNbr = (int) Info[3];
+ int PadNbr = (int) Info[2];
+ //cout << Info[2] << " " << Info[3] << endl;
+ m_Event->SetEnergy(DetectorNbr,PadNbr,Energy);
+ m_Event->SetTime(DetectorNbr,PadNbr,Time);
}
}
// clear map for next event
CaloHitMap->clear();
+
+ ///////////
+ // Veto wall scorer
+ G4THitsMap<G4double*>* VetoHitMap;
+ std::map<G4int, G4double**>::iterator Veto_itr;
+
+ G4int VetoCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("VetoWallScorer/VetoCalorimeter");
+ VetoHitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(VetoCollectionID));
+
+
+ // Loop on the Calo map
+ for (Veto_itr = VetoHitMap->GetMap()->begin() ; Veto_itr != VetoHitMap->GetMap()->end() ; Veto_itr++){
+
+ G4double* Info = *(Veto_itr->second);
+ //(Info[0]/2.35)*((Info[0]*1.02)*pow((Info[0]*1.8),.5))
+ // double Energy = RandGauss::shoot(Info[0],((Info[0]*1000*1.02/2.35)*pow((Info[0]*1000*1.8),.5)) );
+ double Energy = RandGauss::shoot(Info[0],NeutronWall_NS::ResoEnergy);
+ if(Energy>NeutronWall_NS::EnergyThreshold){
+ double Time = RandGauss::shoot(Info[1],NeutronWall_NS::ResoTime);
+ int DetectorNbr = (int) Info[3];
+ int PadNbr = (int) Info[2];
+
+ m_Event->SetVetoEnergy(DetectorNbr,PadNbr,Energy);
+ m_Event->SetVetoTime(DetectorNbr,PadNbr,Time);
+ }
+ }
+ // clear map for next event
+ VetoHitMap->clear();
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////
void NeutronWall::InitializeScorers() {
+ // Otherwise the scorer is initialised
+ vector<int> level;
+ level.push_back(0);
+ level.push_back(1);
+
// This check is necessary in case the geometry is reloaded
bool already_exist = false;
m_NeutronWallScorer = CheckScorer("NeutronWallScorer",already_exist) ;
@@ -417,12 +607,23 @@ void NeutronWall::InitializeScorers() {
if(already_exist)
return ;
- // Otherwise the scorer is initialised
- vector<int> level; level.push_back(0);
- G4VPrimitiveScorer* Calorimeter= new CALORIMETERSCORERS::PS_Calorimeter("Calorimeter",level, 0) ;
- //and register it to the multifunctionnal detector
+ // Neutron Wall Scorer
+ G4VPrimitiveScorer* Calorimeter= new CALORIMETERSCORERS::PS_Calorimeter("Calorimeter",level,1) ;
+ //and register it to the multifunctional detector
m_NeutronWallScorer->RegisterPrimitive(Calorimeter);
G4SDManager::GetSDMpointer()->AddNewDetector(m_NeutronWallScorer) ;
+
+
+ // Veto Wall Scorer
+ already_exist = false;
+ m_VetoWallScorer = CheckScorer("VetoWallScorer",already_exist) ;
+ if(already_exist)
+ return;
+
+ G4VPrimitiveScorer* VetoCalorimeter= new CALORIMETERSCORERS::PS_Calorimeter("VetoCalorimeter",level,1) ;
+ //and register it to the multifunctional detector
+ m_VetoWallScorer->RegisterPrimitive(VetoCalorimeter);
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_VetoWallScorer) ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Detectors/NeutronWall/NeutronWall.hh b/NPSimulation/Detectors/NeutronWall/NeutronWall.hh
index 934ddec2d6a8532cf81a92883c20c07d32b634b9..82576bc2a60715471ff571848d1f94fe1c123609 100644
--- a/NPSimulation/Detectors/NeutronWall/NeutronWall.hh
+++ b/NPSimulation/Detectors/NeutronWall/NeutronWall.hh
@@ -55,7 +55,13 @@ public:
double X,
double Y,
double Z,
- double Rotation);
+ double Rotation,
+ int Bars,
+ string NWMaterial,
+ double VWDistance,
+ int VetoWall,
+ string VWMaterial,
+ double Overlap);
void BuildDetector();
@@ -65,13 +71,15 @@ private:
G4LogicalVolume* m_Quartz_log;
G4LogicalVolume* m_QuartzCap_log;
G4LogicalVolume* m_PMTube_log;
- G4LogicalVolume* m_NE213_log;
+ G4LogicalVolume* m_Scintillator_log;
G4LogicalVolume* m_ShadowBar_log;
+ G4LogicalVolume* m_PlasticBar_log;
- G4VPhysicalVolume* m_NE213Tube_phys;
+ G4VPhysicalVolume* m_ScintillatorTube_phys;
G4VPhysicalVolume* m_Quartz_phys;
G4VPhysicalVolume* m_AlCase_phys;
G4VPhysicalVolume* m_NeutronWall_phys;
+ G4VPhysicalVolume* m_PlasticBar_phys;
////////////////////////////////////////////////////
////// Inherite from NPS::VDetector class /////////
@@ -99,6 +107,7 @@ public: // Scorer
// Associated Scorer
G4MultiFunctionalDetector* m_NeutronWallScorer ;
+ G4MultiFunctionalDetector* m_VetoWallScorer ;
////////////////////////////////////////////////////
///////////Event class to store Data////////////////
////////////////////////////////////////////////////
@@ -117,13 +126,20 @@ private: // Geometry
vector<double> m_Y;
vector<double> m_Z;
vector<double> m_Rot;
+ vector<int> m_Bars;
+ vector<string> m_NWMaterial;
+ vector<double> m_VWDistance;
+ vector<int> m_VetoWall;
+ vector<string> m_VWMaterial;
+ vector<double> m_Overlap;
// Visualisation Attribute
- G4VisAttributes* m_VisNE213;
+ G4VisAttributes* m_VisScintillator;
G4VisAttributes* m_VisQuartz;
G4VisAttributes* m_VisAl;
G4VisAttributes* m_VisNW;
+ G4VisAttributes* m_VisPlasticBar;
// Needed for dynamic loading of the library
public:
diff --git a/Projects/Lassa/Analysis.cxx b/Projects/Lassa/Analysis.cxx
index 7b40f5a52b8986cefb9fa650fdf3f429ebf5caf3..3e560153cc647bbaa08f67b73bb2711174b7f2bf 100644
--- a/Projects/Lassa/Analysis.cxx
+++ b/Projects/Lassa/Analysis.cxx
@@ -40,9 +40,6 @@ void Analysis::Init(){
InitialConditions = new TInitialConditions();
InitOutputBranch();
InitInputBranch();
- totalEvents = 0;
- detectedEvents = 0;
- peakEvents = 0;
f_proton = new TF1("f_proton","1 - TMath::Exp(-(x-0.0918309)/27.7746)",0,10);
f_deuton = new TF1("f_deuton","1 - TMath::Exp(-(x-0.0434552)/21.134)",0,10);
@@ -52,8 +49,14 @@ void Analysis::Init(){
Proton_CsI = EnergyLoss("proton_CsI.G4table","G4Table",100 );
Deuton_CsI = EnergyLoss("deuteron_CsI.G4table","G4Table",100 );
Triton_CsI = EnergyLoss("triton_CsI.G4table","G4Table",100 );
+ He3_CsI = EnergyLoss("He3_CsI.G4table","G4Table",100 );
+
proton = new NPL::Nucleus("1H");
+ deuton = new NPL::Nucleus("2H");
+ triton = new NPL::Nucleus("3H");
+ helium3 = new NPL::Nucleus("3He");
+ alpha = new NPL::Nucleus("4He");
beam = new NPL::Nucleus("112Sn");
target = new NPL::Nucleus("112Sn");
}
@@ -79,21 +82,22 @@ void Analysis::TreatEvent(){
double BetaCM = fEnergyImpulsionLab_total.Beta();
- double ParticleEnergy = InitialConditions->GetKineticEnergy(0);
+ InitialEnergy = InitialConditions->GetKineticEnergy(0);
double EDelta = 2.0;
- InitialEnergy = ParticleEnergy;
- //if(Lassa->ThickSi_E.size()>0) InitialEnergy = ParticleEnergy;
+ if(Lassa->ThickSi_E.size()>0) InitialEnergy_Lassa = InitialEnergy;
double phi_in = acos(InitialConditions->GetMomentumDirectionX(0)/sin(InitialConditions->GetThetaCM(0)*deg));
- if(InitialEnergy>0){
- ECM_initial = proton->GetEnergyCM(InitialEnergy, InitialConditions->GetThetaCM(0)*deg, phi_in, BetaCM);
+
+ ECM_initial = alpha->GetEnergyCM(InitialEnergy, InitialConditions->GetThetaCM(0)*deg, phi_in, BetaCM);
+
+ if(Lassa->ThickSi_E.size()>0){
+ ECM_initial_Lassa = alpha->GetEnergyCM(InitialEnergy_Lassa, InitialConditions->GetThetaCM(0)*deg, phi_in, BetaCM);
}
- else ECM_initial = -100;
+ else ECM_initial_Lassa = -100;
///////////////////////////LOOP on Lassa Hit//////////////////////////////////
if(Lassa->ThickSi_E.size() == 1){
detectedEvents++;
-
//for(unsigned int countLassa = 0 ; countLassa < Lassa->ThickSi_E.size(); countLassa++){
TelescopeNumber = Lassa->TelescopeNumber[0];
@@ -133,6 +137,9 @@ void Analysis::TreatEvent(){
if(Lassa->CsI_E.size()==1){
E_CsI = Lassa->CsI_E[0];
ELab += E_CsI;
+
+ PID = pow(E_ThickSi+E_CsI,1.78)-pow(E_CsI,1.78);
+
//Try to simulate the nuclear reaction loss
//ThicknessCsI = Proton_CsI.EvaluateMaterialThickness(0*MeV, Lassa->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
//ThicknessCsI = Deuton_CsI.EvaluateMaterialThickness(0*MeV, Lassa->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
@@ -148,14 +155,14 @@ void Analysis::TreatEvent(){
}
- if(fabs(ParticleEnergy-ELab)>EDelta){
+ if(fabs(InitialEnergy-ELab)>EDelta){
ELab = -100;
}
- if(fabs(ParticleEnergy-ELab_nucl)>EDelta) ELab_nucl = -100;
+ if(fabs(InitialEnergy-ELab_nucl)>EDelta) ELab_nucl = -100;
if(ELab>0){
- ECM = proton->GetEnergyCM(ELab, ThetaLab, PhiLab, BetaCM);
+ ECM = alpha->GetEnergyCM(ELab, ThetaLab, PhiLab, BetaCM);
}
else ECM = -100;
@@ -189,10 +196,13 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab,"PhiLab/D");
RootOutput::getInstance()->GetTree()->Branch("InitialEnergy",&InitialEnergy,"InitialEnergy/D");
+ RootOutput::getInstance()->GetTree()->Branch("InitialEnergy_Lassa",&InitialEnergy_Lassa,"InitialEnergy_Lassa/D");
RootOutput::getInstance()->GetTree()->Branch("ECM_initial",&ECM_initial,"ECM_initial/D");
+ RootOutput::getInstance()->GetTree()->Branch("ECM_initial_Lassa",&ECM_initial_Lassa,"ECM_initial_Lassa/D");
RootOutput::getInstance()->GetTree()->Branch("E_ThickSi",&E_ThickSi,"E_ThickSi/D");
RootOutput::getInstance()->GetTree()->Branch("E_CsI",&E_CsI,"E_CsI/D");
RootOutput::getInstance()->GetTree()->Branch("R_alpha",&R_alpha,"R_alpha/D");
+ RootOutput::getInstance()->GetTree()->Branch("PID",&PID,"PID/D");
// RootOutput::getInstance()->GetTree()->Branch("peakEvents",&peakEvents,"peakEvents/I");
}
@@ -217,9 +227,12 @@ void Analysis::ReInitValue(){
Z = -100;
TelescopeNumber = -1;
InitialEnergy = -100;
+ InitialEnergy_Lassa = -100;
+ ECM_initial_Lassa = -100;
ECM_initial = -100;
ThicknessCsI = -1;
R_alpha = -100;
+ PID = -1;
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/Lassa/Analysis.h b/Projects/Lassa/Analysis.h
index 56796c57dc7d0e0e6ed656dd585dd1655b827019..14ccdd9f335b1c89db89468fca13a917794d255d 100644
--- a/Projects/Lassa/Analysis.h
+++ b/Projects/Lassa/Analysis.h
@@ -54,10 +54,13 @@ class Analysis: public NPL::VAnalysis{
double TelescopeNumber;
double thresholdEnergy;
double InitialEnergy;
+ double InitialEnergy_Lassa;
double ECM_initial;
+ double ECM_initial_Lassa;
double ThicknessCsI;
double R_alpha;
-
+ double PID;
+
int totalEvents;
int detectedEvents;
int peakEvents;
@@ -76,8 +79,13 @@ class Analysis: public NPL::VAnalysis{
NPL::EnergyLoss Proton_CsI;
NPL::EnergyLoss Deuton_CsI;
NPL::EnergyLoss Triton_CsI;
+ NPL::EnergyLoss He3_CsI;
NPL::Nucleus *proton;
+ NPL::Nucleus *deuton;
+ NPL::Nucleus *triton;
+ NPL::Nucleus *helium3;
+ NPL::Nucleus *alpha;
NPL::Nucleus *beam;
NPL::Nucleus *target;
diff --git a/Projects/Lassa/RunToTreat.txt b/Projects/Lassa/RunToTreat.txt
index ad4003db5e8eb9348ba767b403e2a391f15319fd..e13a324518cca0d05be9944222d179618b3206c2 100755
--- a/Projects/Lassa/RunToTreat.txt
+++ b/Projects/Lassa/RunToTreat.txt
@@ -1,7 +1,7 @@
TTreeName
SimulatedTree
RootFileName
- ../../Outputs/Simulation/test2.root
+ ../../Outputs/Simulation/e09042_alphasource_nucl.root
%../../Outputs/Simulation/lassa_pid_d.root
%../../Outputs/Simulation/lassa_pid_t.root
diff --git a/Projects/macros/GeometricalEfficiency.C b/Projects/macros/GeometricalEfficiency.C
index b82c1ba513d5a86501adfb9f587d713ad7b820a4..67ee4183bd70ab6246e9da0a2e96c93b9af8b594 100644
--- a/Projects/macros/GeometricalEfficiency.C
+++ b/Projects/macros/GeometricalEfficiency.C
@@ -48,7 +48,7 @@
using namespace std;
-void GeometricalEfficiency(const char * fname = "hiraU_flat_nucl"){
+void GeometricalEfficiency(const char * fname = "hiraU_flat_p_nucl"){
// Open output ROOT file from NPTool simulation run
TString path = gSystem->Getenv("NPTOOL");
path += "/Outputs/Simulation/";