diff --git a/.gitignore b/.gitignore
index 33a4376b8db4f13a17178a5f2bc9c74a4dfd7c9d..908556d1c514536da6b9f731b228ff3fc5199e4b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -48,4 +48,5 @@ Projects/T40/19Fdp.reaction
Projects/T40/configs/*.dat
Outputs
Outputs/*
-Outputs/*/*.*
\ No newline at end of file
+Outputs/*/*.*
+NPLib/Utility/npspectra_test
\ No newline at end of file
diff --git a/NPLib/Detectors/LightPipe/TLightPipeData.cxx b/NPLib/Detectors/LightPipe/TLightPipeData.cxx
index 4203e20bc1eb8f234d6f2c1646c7cbcd89d56658..b31f980e4dbaed99e2b4d3fa5e40f92aac2304ae 100644
--- a/NPLib/Detectors/LightPipe/TLightPipeData.cxx
+++ b/NPLib/Detectors/LightPipe/TLightPipeData.cxx
@@ -40,40 +40,4 @@ TLightPipeData::TLightPipeData() {
TLightPipeData::~TLightPipeData() {
}
-
-
//////////////////////////////////////////////////////////////////////
-void TLightPipeData::Clear() {
- // Energy
- fLightPipe_E_DetectorNbr.clear();
- fLightPipe_Energy.clear();
- // Time
- fLightPipe_T_DetectorNbr.clear();
- fLightPipe_Time.clear();
-}
-
-
-
-//////////////////////////////////////////////////////////////////////
-void TLightPipeData::Dump() const {
- // This method is very useful for debuging and worth the dev.
- cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event [TLightPipeData::Dump()] XXXXXXXXXXXXXXXXX" << endl;
-
- // Energy
- size_t mysize = fLightPipe_E_DetectorNbr.size();
- cout << "LightPipe_E_Mult: " << mysize << endl;
-
- for (size_t i = 0 ; i < mysize ; i++){
- cout << "DetNbr: " << fLightPipe_E_DetectorNbr[i]
- << " Energy: " << fLightPipe_Energy[i];
- }
-
- // Time
- mysize = fLightPipe_T_DetectorNbr.size();
- cout << "LightPipe_T_Mult: " << mysize << endl;
-
- for (size_t i = 0 ; i < mysize ; i++){
- cout << "DetNbr: " << fLightPipe_T_DetectorNbr[i]
- << " Time: " << fLightPipe_Time[i];
- }
-}
diff --git a/NPLib/Detectors/LightPipe/TLightPipeData.h b/NPLib/Detectors/LightPipe/TLightPipeData.h
index 659fcbc1dd7194f60fd21a4c525f689100f6fb52..7c1b8e2d7088dbd71a27f0d9f11d68239a48e13a 100644
--- a/NPLib/Detectors/LightPipe/TLightPipeData.h
+++ b/NPLib/Detectors/LightPipe/TLightPipeData.h
@@ -24,7 +24,7 @@
// STL
#include <vector>
-using namespace std;
+#include <iostream>
// ROOT
#include "TObject.h"
@@ -33,29 +33,75 @@ class TLightPipeData : public TObject {
//////////////////////////////////////////////////////////////
// data members are hold into vectors in order
// to allow multiplicity treatment
- private:
- // Energy
- vector<UShort_t> fLightPipe_E_DetectorNbr;
- vector<Double_t> fLightPipe_Energy;
-
- // Time
- vector<UShort_t> fLightPipe_T_DetectorNbr;
- vector<Double_t> fLightPipe_Time;
-
-
+public:
+ // Energy (# of photons)
+ // Left side
+ std::vector<UShort_t> fLightPipe_Left_LayerNbr;
+ std::vector<UShort_t> fLightPipe_Left_RowNbr;
+ std::vector<Double_t> fLightPipe_Left_Energy;
+ // Right side
+ std::vector<UShort_t> fLightPipe_Right_LayerNbr;
+ std::vector<UShort_t> fLightPipe_Right_RowNbr;
+ std::vector<Double_t> fLightPipe_Right_Energy;
+ // Top
+ std::vector<UShort_t> fLightPipe_Top_LayerNbr;
+ std::vector<UShort_t> fLightPipe_Top_RowNbr;
+ std::vector<Double_t> fLightPipe_Top_Energy;
+ // Bottom
+ std::vector<UShort_t> fLightPipe_Bottom_LayerNbr;
+ std::vector<UShort_t> fLightPipe_Bottom_RowNbr;
+ std::vector<Double_t> fLightPipe_Bottom_Energy;
+
+public:
+ enum { Left_t, Right_t, Top_t, Bottom_t };
+
//////////////////////////////////////////////////////////////
// Constructor and destructor
- public:
- TLightPipeData();
- ~TLightPipeData();
+public:
+ TLightPipeData();
+ ~TLightPipeData();
//////////////////////////////////////////////////////////////
// Inherited from TObject and overriden to avoid warnings
- public:
- void Clear();
- void Clear(const Option_t*) {};
- void Dump() const;
+public:
+ void Clear()
+ {
+ fLightPipe_Left_LayerNbr.clear();
+ fLightPipe_Left_RowNbr.clear();
+ fLightPipe_Left_Energy.clear();
+
+ fLightPipe_Right_LayerNbr.clear();
+ fLightPipe_Right_RowNbr.clear();
+ fLightPipe_Right_Energy.clear();
+
+ fLightPipe_Top_LayerNbr.clear();
+ fLightPipe_Top_RowNbr.clear();
+ fLightPipe_Top_Energy.clear();
+
+ fLightPipe_Bottom_LayerNbr.clear();
+ fLightPipe_Bottom_RowNbr.clear();
+ fLightPipe_Bottom_Energy.clear();
+ }
+ void Clear(const Option_t*) {};
+ void Dump() const
+ {
+#if 0
+ // This method is very useful for debuging and worth the dev.
+ std::cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event [TLightPipeData::Dump()] XXXXXXXXXXXXXXXXX\n";
+ // Energy
+ size_t mysize = fLightPipe_E_DetectorNbr.size();
+ std::cout << "LightPipe_E_Mult: " << mysize << "\n";
+ for (size_t i = 0 ; i < mysize ; i++){
+ std::cout << "DetNbr: " << fLightPipe_E_DetectorNbr[i]
+ << "RowNbr (x): " << fLightPipe_E_RowNbr[i]
+ << "ColumnNbr (y): " << fLightPipe_E_ColumnNbr[i]
+ << "DetNbr (z): " << fLightPipe_E_DetectorNbr[i]
+ << " Energy: " << fLightPipe_Energy[i];
+ }
+ std::cout << "\n";
+#endif
+ }
//////////////////////////////////////////////////////////////
@@ -63,42 +109,40 @@ class TLightPipeData : public TObject {
// Prefer inline declaration to avoid unnecessary called of
// frequently used methods
// add //! to avoid ROOT creating dictionnary for the methods
- public:
- ////////////////////// SETTERS ////////////////////////
- // Energy
- inline void SetEnergy(const UShort_t& DetNbr,const Double_t& Energy){
- fLightPipe_E_DetectorNbr.push_back(DetNbr);
- fLightPipe_Energy.push_back(Energy);
- };//!
-
- // Time
- inline void SetTime(const UShort_t& DetNbr,const Double_t& Time) {
- fLightPipe_T_DetectorNbr.push_back(DetNbr);
- fLightPipe_Time.push_back(Time);
- };//!
-
-
- ////////////////////// GETTERS ////////////////////////
- // Energy
- inline UShort_t GetMultEnergy() const
- {return fLightPipe_E_DetectorNbr.size();}
- inline UShort_t GetE_DetectorNbr(const unsigned int &i) const
- {return fLightPipe_E_DetectorNbr[i];}//!
- inline Double_t Get_Energy(const unsigned int &i) const
- {return fLightPipe_Energy[i];}//!
-
- // Time
- inline UShort_t GetMultTime() const
- {return fLightPipe_T_DetectorNbr.size();}
- inline UShort_t GetT_DetectorNbr(const unsigned int &i) const
- {return fLightPipe_T_DetectorNbr[i];}//!
- inline Double_t Get_Time(const unsigned int &i) const
- {return fLightPipe_Time[i];}//!
-
+public:
+ ////////////////////// SETTERS ////////////////////////
+ // Energy
+ void SetEnergy(Int_t side, UShort_t layer, UShort_t row, Double_t energy) {
+ switch(side) {
+ case Left_t:
+ fLightPipe_Left_LayerNbr.push_back(layer);
+ fLightPipe_Left_RowNbr.push_back(row);
+ fLightPipe_Left_Energy.push_back(energy);
+ break;
+ case Right_t:
+ fLightPipe_Right_LayerNbr.push_back(layer);
+ fLightPipe_Right_RowNbr.push_back(row);
+ fLightPipe_Right_Energy.push_back(energy);
+ break;
+ case Top_t:
+ fLightPipe_Top_LayerNbr.push_back(layer);
+ fLightPipe_Top_RowNbr.push_back(row);
+ fLightPipe_Top_Energy.push_back(energy);
+ break;
+ case Bottom_t:
+ fLightPipe_Bottom_LayerNbr.push_back(layer);
+ fLightPipe_Bottom_RowNbr.push_back(row);
+ fLightPipe_Bottom_Energy.push_back(energy);
+ break;
+ default:
+ std::cerr << "WARNING: Invalid side: " << side << " specified to TLightPipeData::SetEnergy(); nothing will be set!\n";
+ break;
+ }
+ };//!
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
- ClassDef(TLightPipeData,1) // LightPipeData structure
+ ClassDef(TLightPipeData,1); // LightPipeData structure
};
#endif
diff --git a/NPLib/Detectors/LightPipe/TLightPipePhysics.cxx b/NPLib/Detectors/LightPipe/TLightPipePhysics.cxx
index dc81b1ab50cee0b636f881ecfa5f561fe38a27f2..34553e07c8e648abea55b8e291eb5cbc4918a641 100644
--- a/NPLib/Detectors/LightPipe/TLightPipePhysics.cxx
+++ b/NPLib/Detectors/LightPipe/TLightPipePhysics.cxx
@@ -28,6 +28,8 @@
#include <cmath>
#include <stdlib.h>
#include <limits>
+#include <algorithm>
+#include <iterator>
using namespace std;
// NPL
@@ -81,19 +83,13 @@ void TLightPipePhysics::BuildSimplePhysicalEvent() {
void TLightPipePhysics::BuildPhysicalEvent() {
// apply thresholds and calibration
PreTreat();
-
- // match energy and time together
- unsigned int mysizeE = m_PreTreatedData->GetMultEnergy();
- unsigned int mysizeT = m_PreTreatedData->GetMultTime();
- for (UShort_t e = 0; e < mysizeE ; e++) {
- for (UShort_t t = 0; t < mysizeT ; t++) {
- if (m_PreTreatedData->GetE_DetectorNbr(e) == m_PreTreatedData->GetT_DetectorNbr(t)) {
- DetectorNumber.push_back(m_PreTreatedData->GetE_DetectorNbr(e));
- Energy.push_back(m_PreTreatedData->Get_Energy(e));
- Time.push_back(m_PreTreatedData->Get_Time(t));
- }
- }
- }
+#if 0
+ // copy calibrated data
+ copy(m_PreTreatedData->GetE_RowNbr().begin(), m_PreTreatedData->GetE_RowNbr().end(), back_inserter(RowNumber));
+ copy(m_PreTreatedData->GetE_ColumnNbr().begin(), m_PreTreatedData->GetE_ColumnNbr().end(), back_inserter(ColumnNumber));
+ copy(m_PreTreatedData->GetE_DetectorNbr().begin(), m_PreTreatedData->GetE_DetectorNbr().end(), back_inserter(DetectorNumber));
+ copy(m_PreTreatedData->Get_Energy().begin(), m_PreTreatedData->Get_Energy().end(), back_inserter(Energy));
+#endif
}
///////////////////////////////////////////////////////////////////////////
@@ -107,23 +103,18 @@ void TLightPipePhysics::PreTreat() {
// instantiate CalibrationManager
static CalibrationManager* Cal = CalibrationManager::getInstance();
+#if 0
// Energy
unsigned int mysize = m_EventData->GetMultEnergy();
for (UShort_t i = 0; i < mysize ; ++i) {
if (m_EventData->Get_Energy(i) > m_E_RAW_Threshold) {
- Double_t Energy = Cal->ApplyCalibration("LightPipe/ENERGY"+NPL::itoa(m_EventData->GetE_DetectorNbr(i)),m_EventData->Get_Energy(i));
+ Double_t Energy = m_EventData->Get_Energy(i); // ORIGINAL:: Cal->ApplyCalibration("LightPipe/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_RowNbr(i), m_EventData->GetE_ColumnNbr(i), m_EventData->GetE_DetectorNbr(i), Energy);
}
}
}
-
- // Time
- mysize = m_EventData->GetMultTime();
- for (UShort_t i = 0; i < mysize; ++i) {
- Double_t Time= Cal->ApplyCalibration("LightPipe/TIME"+NPL::itoa(m_EventData->GetT_DetectorNbr(i)),m_EventData->Get_Time(i));
- m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i), Time);
- }
+#endif
}
@@ -194,9 +185,10 @@ void TLightPipePhysics::ReadAnalysisConfig() {
///////////////////////////////////////////////////////////////////////////
void TLightPipePhysics::Clear() {
+ RowNumber.clear();
+ ColumnNumber.clear();
DetectorNumber.clear();
Energy.clear();
- Time.clear();
}
diff --git a/NPLib/Detectors/LightPipe/TLightPipePhysics.h b/NPLib/Detectors/LightPipe/TLightPipePhysics.h
index 01932fea9df41325bf11a9713dc25ceef96c9c64..6f241ffad80b8dba9f867366c309a81d30efd44e 100644
--- a/NPLib/Detectors/LightPipe/TLightPipePhysics.h
+++ b/NPLib/Detectors/LightPipe/TLightPipePhysics.h
@@ -46,25 +46,26 @@ class TLightPipeSpectra;
class TLightPipePhysics : public TObject, public NPL::VDetector {
//////////////////////////////////////////////////////////////
// constructor and destructor
- public:
- TLightPipePhysics();
- ~TLightPipePhysics() {};
+public:
+ TLightPipePhysics();
+ ~TLightPipePhysics() {};
//////////////////////////////////////////////////////////////
// Inherited from TObject and overriden to avoid warnings
- public:
- void Clear();
- void Clear(const Option_t*) {};
+public:
+ void Clear();
+ void Clear(const Option_t*) {};
//////////////////////////////////////////////////////////////
// data obtained after BuildPhysicalEvent() and stored in
// output ROOT file
- public:
- vector<int> DetectorNumber;
- vector<double> Energy;
- vector<double> Time;
+public:
+ vector<int> RowNumber; // x
+ vector<int> ColumnNumber; // y
+ vector<int> DetectorNumber; // z
+ vector<double> Energy;
/// A usefull method to bundle all operation to add a detector
void AddDetector(TVector3 POS, string shape);
@@ -72,109 +73,109 @@ class TLightPipePhysics : public TObject, public NPL::VDetector {
//////////////////////////////////////////////////////////////
// methods inherited from the VDetector ABC class
- public:
- // read stream from ConfigFile to pick-up detector parameters
- void ReadConfiguration(NPL::InputParser);
+public:
+ // read stream from ConfigFile to pick-up detector parameters
+ void ReadConfiguration(NPL::InputParser);
- // add parameters to the CalibrationManger
- void AddParameterToCalibrationManager();
+ // add parameters to the CalibrationManger
+ void AddParameterToCalibrationManager();
- // method called event by event, aiming at extracting the
- // physical information from detector
- void BuildPhysicalEvent();
+ // method called event by event, aiming at extracting the
+ // physical information from detector
+ void BuildPhysicalEvent();
- // same as BuildPhysicalEvent() method but with a simpler
- // treatment
- void BuildSimplePhysicalEvent();
+ // same as BuildPhysicalEvent() method but with a simpler
+ // treatment
+ void BuildSimplePhysicalEvent();
- // same as above but for online analysis
- void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
+ // same as above but for online analysis
+ void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
- // activate raw data object and branches from input TChain
- // in this method mother branches (Detector) AND daughter leaves
- // (fDetector_parameter) have to be activated
- void InitializeRootInputRaw();
+ // activate raw data object and branches from input TChain
+ // in this method mother branches (Detector) AND daughter leaves
+ // (fDetector_parameter) have to be activated
+ void InitializeRootInputRaw();
- // activate physics data object and branches from input TChain
- // in this method mother branches (Detector) AND daughter leaves
- // (fDetector_parameter) have to be activated
- void InitializeRootInputPhysics();
+ // activate physics data object and branches from input TChain
+ // in this method mother branches (Detector) AND daughter leaves
+ // (fDetector_parameter) have to be activated
+ void InitializeRootInputPhysics();
- // create branches of output ROOT file
- void InitializeRootOutput();
+ // create branches of output ROOT file
+ void InitializeRootOutput();
- // clear the raw and physical data objects event by event
- void ClearEventPhysics() {Clear();}
- void ClearEventData() {m_EventData->Clear();}
+ // clear the raw and physical data objects event by event
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() {m_EventData->Clear();}
- // methods related to the TLightPipeSpectra class
- // instantiate the TLightPipeSpectra class and
- // declare list of histograms
- void InitSpectra();
+ // methods related to the TLightPipeSpectra class
+ // instantiate the TLightPipeSpectra class and
+ // declare list of histograms
+ void InitSpectra();
- // fill the spectra
- void FillSpectra();
+ // fill the spectra
+ void FillSpectra();
- // used for Online mainly, sanity check for histograms and
- // change their color if issues are found, for example
- void CheckSpectra();
+ // used for Online mainly, sanity check for histograms and
+ // change their color if issues are found, for example
+ void CheckSpectra();
- // used for Online only, clear all the spectra
- void ClearSpectra();
+ // used for Online only, clear all the spectra
+ void ClearSpectra();
- // write spectra to ROOT output file
- void WriteSpectra();
+ // write spectra to ROOT output file
+ void WriteSpectra();
//////////////////////////////////////////////////////////////
// specific methods to LightPipe array
- public:
- // remove bad channels, calibrate the data and apply thresholds
- void PreTreat();
+public:
+ // remove bad channels, calibrate the data and apply thresholds
+ void PreTreat();
- // clear the pre-treated object
- void ClearPreTreatedData() {m_PreTreatedData->Clear();}
+ // clear the pre-treated object
+ void ClearPreTreatedData() {m_PreTreatedData->Clear();}
- // read the user configuration file. If no file is found, load standard one
- void ReadAnalysisConfig();
+ // read the user configuration file. If no file is found, load standard one
+ void ReadAnalysisConfig();
- // give and external TLightPipeData object to TLightPipePhysics.
- // needed for online analysis for example
- void SetRawDataPointer(TLightPipeData* rawDataPointer) {m_EventData = rawDataPointer;}
+ // give and external TLightPipeData object to TLightPipePhysics.
+ // needed for online analysis for example
+ void SetRawDataPointer(TLightPipeData* rawDataPointer) {m_EventData = rawDataPointer;}
// objects are not written in the TTree
- private:
- TLightPipeData* m_EventData; //!
- TLightPipeData* m_PreTreatedData; //!
- TLightPipePhysics* m_EventPhysics; //!
+private:
+ TLightPipeData* m_EventData; //!
+ TLightPipeData* m_PreTreatedData; //!
+ TLightPipePhysics* m_EventPhysics; //!
// getters for raw and pre-treated data object
- public:
- TLightPipeData* GetRawData() const {return m_EventData;}
- TLightPipeData* GetPreTreatedData() const {return m_PreTreatedData;}
+public:
+ TLightPipeData* GetRawData() const {return m_EventData;}
+ TLightPipeData* GetPreTreatedData() const {return m_PreTreatedData;}
// parameters used in the analysis
- private:
- // thresholds
- double m_E_RAW_Threshold; //!
- double m_E_Threshold; //!
+private:
+ // thresholds
+ double m_E_RAW_Threshold; //!
+ double m_E_Threshold; //!
// number of detectors
- private:
- int m_NumberOfDetectors; //!
+private:
+ int m_NumberOfDetectors; //!
// spectra class
- private:
- TLightPipeSpectra* m_Spectra; // !
+private:
+ TLightPipeSpectra* m_Spectra; // !
// spectra getter
- public:
- map<string, TH1*> GetSpectra();
+public:
+ map<string, TH1*> GetSpectra();
// Static constructor to be passed to the Detector Factory
- public:
- static NPL::VDetector* Construct();
+public:
+ static NPL::VDetector* Construct();
- ClassDef(TLightPipePhysics,1) // LightPipePhysics structure
+ ClassDef(TLightPipePhysics,1) // LightPipePhysics structure
};
#endif
diff --git a/NPLib/Detectors/LightPipe/TLightPipeSpectra.cxx b/NPLib/Detectors/LightPipe/TLightPipeSpectra.cxx
index dcd074d9a3e6b4ba753b8f8da92fd0b6289edbc5..26e904c5e0a1f35a7ac51e66eee4fab1b511f185 100644
--- a/NPLib/Detectors/LightPipe/TLightPipeSpectra.cxx
+++ b/NPLib/Detectors/LightPipe/TLightPipeSpectra.cxx
@@ -108,67 +108,67 @@ void TLightPipeSpectra::InitPhysicsSpectra() {
////////////////////////////////////////////////////////////////////////////////
void TLightPipeSpectra::FillRawSpectra(TLightPipeData* RawData) {
- static string name;
- static string family;
-
- // Energy
- unsigned int sizeE = RawData->GetMultEnergy();
- for (unsigned int i = 0; i < sizeE; i++) {
- name = "LightPipe"+NPL::itoa(RawData->GetE_DetectorNbr(i))+"_ENERGY_RAW";
- family = "LightPipe/RAW";
-
- FillSpectra(family,name,RawData->Get_Energy(i));
- }
-
- // Time
- unsigned int sizeT = RawData->GetMultTime();
- for (unsigned int i = 0; i < sizeT; i++) {
- name = "LightPipe"+NPL::itoa(RawData->GetT_DetectorNbr(i))+"_TIME_RAW";
- family = "LightPipe/RAW";
-
- FillSpectra(family,name,RawData->Get_Time(i));
- }
+ // static string name;
+ // static string family;
+
+ // // Energy
+ // unsigned int sizeE = RawData->GetMultEnergy();
+ // for (unsigned int i = 0; i < sizeE; i++) {
+ // name = "LightPipe"+NPL::itoa(RawData->GetE_DetectorNbr(i))+"_ENERGY_RAW";
+ // family = "LightPipe/RAW";
+
+ // FillSpectra(family,name,RawData->Get_Energy(i));
+ // }
+
+ // // Time
+ // unsigned int sizeT = RawData->GetMultTime();
+ // for (unsigned int i = 0; i < sizeT; i++) {
+ // name = "LightPipe"+NPL::itoa(RawData->GetT_DetectorNbr(i))+"_TIME_RAW";
+ // family = "LightPipe/RAW";
+
+ // FillSpectra(family,name,RawData->Get_Time(i));
+ // }
}
////////////////////////////////////////////////////////////////////////////////
void TLightPipeSpectra::FillPreTreatedSpectra(TLightPipeData* PreTreatedData) {
- static string name;
- static string family;
+ // static string name;
+ // static string family;
- // Energy
- unsigned int sizeE = PreTreatedData->GetMultEnergy();
- for (unsigned int i = 0; i < sizeE; i++) {
- name = "LightPipe"+NPL::itoa(PreTreatedData->GetE_DetectorNbr(i))+"_ENERGY_CAL";
- family = "LightPipe/CAL";
-
- FillSpectra(family,name,PreTreatedData->Get_Energy(i));
- }
-
- // Time
- unsigned int sizeT = PreTreatedData->GetMultTime();
- for (unsigned int i = 0; i < sizeT; i++) {
- name = "LightPipe"+NPL::itoa(PreTreatedData->GetT_DetectorNbr(i))+"_TIME_CAL";
- family = "LightPipe/CAL";
-
- FillSpectra(family,name,PreTreatedData->Get_Time(i));
- }
+ // // Energy
+ // unsigned int sizeE = PreTreatedData->GetMultEnergy();
+ // for (unsigned int i = 0; i < sizeE; i++) {
+ // name = "LightPipe"+NPL::itoa(PreTreatedData->GetE_DetectorNbr(i))+"_ENERGY_CAL";
+ // family = "LightPipe/CAL";
+
+ // FillSpectra(family,name,PreTreatedData->Get_Energy(i));
+ // }
+
+ // // Time
+ // unsigned int sizeT = PreTreatedData->GetMultTime();
+ // for (unsigned int i = 0; i < sizeT; i++) {
+ // name = "LightPipe"+NPL::itoa(PreTreatedData->GetT_DetectorNbr(i))+"_TIME_CAL";
+ // family = "LightPipe/CAL";
+
+ // FillSpectra(family,name,PreTreatedData->Get_Time(i));
+ // }
}
////////////////////////////////////////////////////////////////////////////////
void TLightPipeSpectra::FillPhysicsSpectra(TLightPipePhysics* Physics) {
- static string name;
- static string family;
- family= "LightPipe/PHY";
-
- // Energy vs time
- unsigned int sizeE = Physics->Energy.size();
- for(unsigned int i = 0 ; i < sizeE ; i++){
- name = "LightPipe_ENERGY_TIME";
- FillSpectra(family,name,Physics->Energy[i],Physics->Time[i]);
- }
+ // static string name;
+ // static string family;
+ // family= "LightPipe/PHY";
+
+ // // Energy vs time
+ // unsigned int sizeE = Physics->Energy.size();
+ // for(unsigned int i = 0 ; i < sizeE ; i++){
+ // name = "LightPipe_ENERGY_TIME";
+ // FillSpectra(family,name,Physics->Energy[i],Physics->Time[i]);
+ // }
}
diff --git a/NPSimulation/Detectors/LightPipe/LightPipe.cc b/NPSimulation/Detectors/LightPipe/LightPipe.cc
index b058407638661d993828feee4cfce10e38fc4aac..f29b9ac21eac8f226e375a8fcdc3e0ebed98f947 100644
--- a/NPSimulation/Detectors/LightPipe/LightPipe.cc
+++ b/NPSimulation/Detectors/LightPipe/LightPipe.cc
@@ -1,18 +1,18 @@
/*****************************************************************************
- * Copyright (C) 2009-2018 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2018 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: Greg Christian contact address: gchristian@tamu.edu *
+ * Original Author: Greg Christian contact address: gchristian@tamu.edu *
* *
- * Creation Date : July 2018 *
+ * Creation Date : July 2018 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * This class describe LightPipe simulation *
+ * This class describe LightPipe simulation *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -23,6 +23,8 @@
#include <sstream>
#include <cmath>
#include <limits>
+#include <fstream>
+#include <algorithm>
//G4 Geometry object
#include "G4Tubs.hh"
#include "G4Box.hh"
@@ -37,10 +39,18 @@
#include "G4PVPlacement.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
+#include "G4OpticalSurface.hh"
+#include "G4LogicalBorderSurface.hh"
+#include "G4LogicalSkinSurface.hh"
+
+// ROOT
+#include "TSystem.h"
+#include "TMath.h"
// NPTool header
#include "LightPipe.hh"
#include "CalorimeterScorers.hh"
+#include "PhotoDiodeScorers.hh"
#include "RootOutput.h"
#include "MaterialManager.hh"
#include "NPSDetectorFactory.hh"
@@ -53,120 +63,74 @@ using namespace std;
using namespace CLHEP;
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-namespace LightPipe_NS{
- // Energy and time Resolution
- const double EnergyThreshold = 0.1*MeV;
- const double ResoTime = 4.5*ns ;
- const double ResoEnergy = 1.0*MeV ;
- const double Radius = 50*mm ;
- const double Width = 100*mm ;
- const double Thickness = 300*mm ;
- const string Material = "BC400";
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// LightPipe Specific Method
LightPipe::LightPipe(){
m_Event = new TLightPipeData() ;
m_LightPipeScorer = 0;
- m_SquareDetector = 0;
- m_CylindricalDetector = 0;
-
// RGB Color + Transparency
- m_VisSquare = new G4VisAttributes(G4Colour(0, 1, 0, 0.5));
- m_VisCylinder = new G4VisAttributes(G4Colour(0, 0, 1, 0.5));
-
+ m_VisSquare = new G4VisAttributes(G4Colour(0, 0, 1, 0.5));
+ m_VisPipe = new G4VisAttributes(G4Colour(0, 1, 0, 0.5));
+ m_VisPD = new G4VisAttributes(G4Colour(0.1, 0.2, 0.3));
+ m_ScintillatorMaterial = CreateScintillatorMaterial();
+ m_PipeMaterial = CreatePipeMaterial();
+ //m_Wrapping = CreateWrappingMaterial();
+ m_ReflectiveSurface = CreateReflectiveSurface();
+
+ m_VisSquare->SetForceWireframe(true);
+ m_VisPipe->SetForceWireframe(true);
}
LightPipe::~LightPipe(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void LightPipe::AddDetector(G4ThreeVector POS, string Shape){
- // Convert the POS value to R theta Phi as Spherical coordinate is easier in G4
- m_R.push_back(POS.mag());
- m_Theta.push_back(POS.theta());
- m_Phi.push_back(POS.phi());
- m_Shape.push_back(Shape);
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void LightPipe::AddDetector(double R, double Theta, double Phi, string Shape){
- m_R.push_back(R);
- m_Theta.push_back(Theta);
- m_Phi.push_back(Phi);
- m_Shape.push_back(Shape);
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* LightPipe::BuildSquareDetector(){
- if(!m_SquareDetector){
- G4Box* box = new G4Box("LightPipe_Box",LightPipe_NS::Width*0.5,
- LightPipe_NS::Width*0.5,LightPipe_NS::Thickness*0.5);
-
- G4Material* DetectorMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(LightPipe_NS::Material);
- m_SquareDetector = new G4LogicalVolume(box,DetectorMaterial,"logic_LightPipe_Box",0,0,0);
- m_SquareDetector->SetVisAttributes(m_VisSquare);
- m_SquareDetector->SetSensitiveDetector(m_LightPipeScorer);
- }
- return m_SquareDetector;
+void LightPipe::AddDetector(G4int nrow, G4int ncol, G4int nlayer, G4double width, G4double thickness, G4double pipe_width, G4double pipe_thickness){
+ m_Detector.emplace_back(make_tuple(nrow,ncol,nlayer,width,thickness,pipe_width,pipe_thickness));
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* LightPipe::BuildCylindricalDetector(){
- if(!m_CylindricalDetector){
- G4Tubs* tub = new G4Tubs("LightPipe_Cyl",0,LightPipe_NS::Radius,LightPipe_NS::Thickness*0.5,0,360*deg);
- G4Material* DetectorMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(LightPipe_NS::Material);
- m_CylindricalDetector = new G4LogicalVolume(tub,DetectorMaterial,"logic_LightPipe_tub",0,0,0);
- m_CylindricalDetector->SetVisAttributes(m_VisSquare);
- m_CylindricalDetector->SetSensitiveDetector(m_LightPipeScorer);
- }
- return m_CylindricalDetector;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Virtual Method of NPS::VDetector class
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
void LightPipe::ReadConfiguration(NPL::InputParser parser){
+
+ G4double width = 10*mm;
+ G4double thickness = 10*mm;
+ G4double pipe_width = 3*mm;
+ G4double pipe_thickness = 1*mm;
+ G4int nrow = 1;
+ G4int ncol = 1;
+ G4int nlayer = 1;
+
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("LightPipe");
- if(NPOptionManager::getInstance()->GetVerboseLevel())
+ if(true || NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
-
- vector<string> cart = {"POS","Shape"};
- vector<string> sphe = {"R","Theta","Phi","Shape"};
+ vector<string> lp = {"WIDTH","THICKNESS","PIPE_WIDTH","PIPE_THICKNESS","NROW","NCOL","NLAYER"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(cart)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// LightPipe " << i+1 << endl;
-
- G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
- string Shape = blocks[i]->GetString("Shape");
- AddDetector(Pos,Shape);
- }
- else if(blocks[i]->HasTokenList(sphe)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
+ if(blocks[i]->HasTokenList(lp)){
+ if(true || NPOptionManager::getInstance()->GetVerboseLevel()){
cout << endl << "//// LightPipe " << i+1 << endl;
- double R = blocks[i]->GetDouble("R","mm");
- double Theta = blocks[i]->GetDouble("Theta","deg");
- double Phi = blocks[i]->GetDouble("Phi","deg");
- string Shape = blocks[i]->GetString("Shape");
- AddDetector(R,Theta,Phi,Shape);
+ }
+ width = blocks[i]->GetDouble("WIDTH", "mm");
+ thickness = blocks[i]->GetDouble("THICKNESS", "mm");
+ pipe_width = blocks[i]->GetDouble("PIPE_WIDTH", "mm");
+ pipe_thickness = blocks[i]->GetDouble("PIPE_THICKNESS", "mm");
+ nrow = blocks[i]->GetInt("NROW");
+ ncol = blocks[i]->GetInt("NCOL");
+ nlayer = blocks[i]->GetInt("NLAYER");
}
else{
cout << "ERROR: check your input file formatting " << endl;
exit(1);
}
}
+ //
+ AddDetector(nrow, ncol, nlayer, width, thickness, pipe_width, pipe_thickness);
}
@@ -175,39 +139,230 @@ void LightPipe::ReadConfiguration(NPL::InputParser parser){
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
void LightPipe::ConstructDetector(G4LogicalVolume* world){
- for (unsigned short i = 0 ; i < m_R.size() ; i++) {
-
- G4double wX = m_R[i] * sin(m_Theta[i] ) * cos(m_Phi[i] ) ;
- G4double wY = m_R[i] * sin(m_Theta[i] ) * sin(m_Phi[i] ) ;
- G4double wZ = m_R[i] * cos(m_Theta[i] ) ;
- G4ThreeVector Det_pos = G4ThreeVector(wX, wY, wZ) ;
- // So the face of the detector is at R instead of the middle
- Det_pos+=Det_pos.unit()*LightPipe_NS::Thickness*0.5;
- // Building Detector reference frame
- G4double ii = cos(m_Theta[i]) * cos(m_Phi[i]);
- G4double jj = cos(m_Theta[i]) * sin(m_Phi[i]);
- G4double kk = -sin(m_Theta[i]);
- G4ThreeVector Y(ii,jj,kk);
- G4ThreeVector w = Det_pos.unit();
- G4ThreeVector u = w.cross(Y);
- G4ThreeVector v = w.cross(u);
- v = v.unit();
- u = u.unit();
-
- G4RotationMatrix* Rot = new G4RotationMatrix(u,v,w);
-
- if(m_Shape[i] == "Cylindrical"){
- new G4PVPlacement(G4Transform3D(*Rot,Det_pos),
- BuildCylindricalDetector(),
- "LightPipe",world,false,i+1);
- }
+ bool warnOverlap = false;
- else if(m_Shape[i] == "Square"){
- new G4PVPlacement(G4Transform3D(*Rot,Det_pos),
- BuildSquareDetector(),
- "LightPipe",world,false,i+1);
- }
- }
+
+ //Create experimental hall
+ G4Material* Vacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum") ;
+ G4double expHall_x = 10.*m;
+ G4double expHall_y = 10.*m;
+ G4double expHall_z = 10.*m;
+
+ G4Box* fExperimentalHall_box = new G4Box("expHall_box",expHall_x,expHall_y,expHall_z);
+ G4LogicalVolume* fExperimentalHall_log = new G4LogicalVolume(fExperimentalHall_box,
+ Vacuum,"expHall_log",0,0,0);
+ G4VPhysicalVolume* fExperimentalHall_phys = new G4PVPlacement(0,G4ThreeVector(),
+ fExperimentalHall_log,"expHall",0,false,0);
+
+ fExperimentalHall_log->SetVisAttributes(G4VisAttributes::Invisible);
+
+ auto BuildRectangle = [this](G4double width, G4double length, G4double thickness, G4Material* material){
+ G4Box* box = new G4Box("LightPipe_Box",width*0.5,
+ length*0.5,thickness*0.5);
+ G4LogicalVolume* Detector = new G4LogicalVolume(box,material,"logic_LightPipe_Box",0,0,0);
+ Detector->SetVisAttributes(this->m_VisSquare);
+// Detector->SetSensitiveDetector(this->m_LightPipeScorer);
+ return Detector;
+ };
+ auto getCenter = [&](int i, int imax, double width){
+ return (i - (imax/2.))*width + width/2.;
+ };
+
+ int i=0, j=0, k=0;
+
+ int iPipeX=1, iPipeY=1, iDet=1;
+ for(const auto& det : m_Detector) {
+ const G4int& nrow = get<0>(det);
+ const G4int& ncol = get<1>(det);
+ const G4int& nlayer = get<2>(det);
+ const G4double& width = get<3>(det);
+ const G4double& thickness = get<4>(det);
+ const G4double& pipe_width = get<5>(det);
+ const G4double& pipe_thickness = get<6>(det);
+ const G4double pd_thickness = 1*mm;
+
+ vector<vector<G4PVPlacement*> > physVol(nrow);
+ for(auto& v : physVol) { v.resize(ncol); }
+
+
+ auto buildRow = [&](G4int irow, G4double z){
+ G4double rowWidthX = nrow*width;
+ G4double pipe_length = width*ncol + 1*cm;
+ //
+ // Build light pipe above detectors
+ //
+ // Create geometric object
+ G4ThreeVector pipePos(
+ 0, getCenter(irow, nrow, width) + width/2. + pipe_thickness/2., z);
+ auto pipe = BuildRectangle(
+ pipe_length, pipe_width, pipe_thickness, m_PipeMaterial);
+ pipe->SetVisAttributes(m_VisPipe);
+ // Rotate it
+ G4RotationMatrix* myRotation = new G4RotationMatrix();
+ myRotation->rotateX(90.*deg);
+ // Create PV Placement
+ G4PVPlacement* pv = new G4PVPlacement(
+ myRotation, pipePos, pipe, "LightPipe_PipeX", world, false, iPipeX++, warnOverlap);
+
+ std::vector<G4PVPlacement*> pvRow;
+ for(int icol=0; icol< ncol; ++icol){
+ //
+ // Build row of detectors
+ //
+ //
+ // Create geometric object
+ G4ThreeVector Det_pos
+ (getCenter(icol,ncol,width), getCenter(irow,nrow,width), z);
+ auto Scintillator = BuildRectangle(width, width, thickness, m_ScintillatorMaterial);
+// Scintillator->SetSensitiveDetector(this->m_LightPipeScorer);
+ // Create PV placement
+ pvRow.push_back(
+ new G4PVPlacement(
+ 0, Det_pos, Scintillator, "LightPipe_Detector", world, false, iDet++, warnOverlap) );
+ }
+ // Create reflective surfaces between detectors
+ for(int icol=0; icol< ncol; ++icol){
+ // to the left
+ if(icol != 0) {
+ new G4LogicalBorderSurface("CrystalSurface", pvRow.at(icol-1), pvRow.at(icol), m_ReflectiveSurface);
+ } else {
+ new G4LogicalBorderSurface("CrystalSurface", fExperimentalHall_phys, pvRow.at(icol), m_ReflectiveSurface);
+ }
+ // to the right
+ if(icol != ncol-1) {
+ new G4LogicalBorderSurface("CrystalSurface", pvRow.at(icol), pvRow.at(icol+1), m_ReflectiveSurface);
+ } else {
+ new G4LogicalBorderSurface("CrystalSurface", pvRow.at(icol), fExperimentalHall_phys, m_ReflectiveSurface);
+ }
+ }
+ return pvRow;
+ };
+
+ buildRow(5,0);
+#if 0
+
+ int pdNum = 1;
+ for(int ilayer = 0; ilayer< nlayer; ++ilayer){
+ G4double detZ = getCenter(ilayer,nlayer,thickness) + pipe_thickness*ilayer;
+ for(int irow = 0; irow< nrow; ++irow){
+ for(int icol = 0; icol< ncol; ++icol){
+ //
+ // Build detectors
+ G4ThreeVector Det_pos
+ (getCenter(irow,nrow,width),
+ getCenter(icol,ncol,width),
+ detZ);
+ auto Rect = BuildRectangle(width, width, thickness, m_ScintillatorMaterial);
+ auto pvScint = new G4PVPlacement(0, Det_pos, Rect,"LightPipe_Detector",world,false,i+1,warnOverlap);
+ physVol.at(irow).at(icol) = pvScint;
+ ++i;
+ } // for(icol...)
+ } // for(irow...)
+ //
+ // Reflection between detectors
+ for(int irow = -1; irow< nrow; ++irow){
+ for(int icol = -1; icol< ncol; ++icol){
+ if(icol > -1) {
+ auto pvLeft = irow != -1 ? physVol.at(irow).at(icol) : fExperimentalHall_phys;
+ auto pvRight = irow != nrow-1 ? physVol.at(irow+1).at(icol) : fExperimentalHall_phys;
+ new G4LogicalBorderSurface("CrystalSurface", pvLeft, pvRight, m_ReflectiveSurface);
+ new G4LogicalBorderSurface("CrystalSurface", pvRight, pvLeft, m_ReflectiveSurface);
+ }
+ if(irow > -1) {
+ auto pvUp = icol != -1 ? physVol.at(irow).at(icol) : fExperimentalHall_phys;
+ auto pvDn = icol != ncol - 1 ? physVol.at(irow).at(icol+1) : fExperimentalHall_phys;
+ new G4LogicalBorderSurface("CrystalSurface", pvUp, pvDn, m_ReflectiveSurface);
+ new G4LogicalBorderSurface("CrystalSurface", pvDn, pvUp, m_ReflectiveSurface);
+ }
+ //new G4LogicalSkinSurface("m_ReflectiveSurface",logicCsI,m_ReflectiveSurface);
+ }
+ }
+ //
+ // Build Pipes
+ auto buildPipe = [&](bool horizontal, G4double layer_center) {
+ G4PVPlacement* pvLast = 0;
+ G4double pipe_center_z = layer_center + thickness*0.5 + pipe_thickness*0.5;
+
+ int ircmax = horizontal ? ncol : nrow;
+ for(int irc=0; irc< ircmax; ++irc){
+ G4double pipe_center_xy = getCenter(irc,ircmax,width);
+ G4ThreeVector LG_pos = horizontal ?
+ G4ThreeVector(0,pipe_center_xy-thickness/2.,pipe_center_z-thickness/2.) :
+ G4ThreeVector(pipe_center_xy,0,pipe_center_z);
+ G4Material* mat = m_PipeMaterial;
+
+ G4double pipe_length = horizontal ? width*nrow : width*ncol;
+ pipe_length += 1*cm;
+ auto pipe = horizontal ?
+ BuildRectangle(pipe_length, pipe_width, pipe_thickness, mat) :
+ BuildRectangle(pipe_width, pipe_length, pipe_thickness, mat) ;
+ pipe->SetVisAttributes(m_VisPipe);
+ G4PVPlacement* pv = 0;
+ if(horizontal){
+ G4RotationMatrix* myRotation = new G4RotationMatrix();
+ myRotation->rotateX(90.*deg);
+ pv = new G4PVPlacement(myRotation, LG_pos, pipe, "LightPipe_PipeX",world,false, j+1, warnOverlap);
+ ++j;
+ } else {
+ pv = new G4PVPlacement(0, LG_pos, pipe, "LightPipe_PipeY",world,false, k+1, warnOverlap);
+ ++k;
+ }
+ new G4LogicalBorderSurface("CrystalSurface", pv, fExperimentalHall_phys, m_ReflectiveSurface);
+ if(pvLast) {
+ new G4LogicalBorderSurface("CrystalSurface", pvLast, pv, m_ReflectiveSurface);
+ new G4LogicalBorderSurface("CrystalSurface", pv, pvLast, m_ReflectiveSurface);
+ }
+ pvLast = pv;
+#if 0
+ int icrmax = horizontal ? nrow : ncol;
+ for(int icr = 0; icr< icrmax; ++icr) {
+ auto pvDet = horizontal ? physVol.at(icr).at(irc) : physVol.at(irc).at(icr);
+ new G4LogicalBorderSurface("CrystalSurface", pv, pvDet, m_ReflectiveSurface);
+ }
+#endif
+#if 0
+ //
+ // Photodiode
+ G4String NamePD = "PhotoDiode";
+ G4Material* PDMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
+ G4Box* solidPhotoDiode = horizontal ?
+ new G4Box(NamePD, 0.5*pipe_thickness, 0.5*pipe_width, 0.5*pd_thickness) :
+ new G4Box(NamePD, 0.5*pipe_width, 0.5*pipe_thickness, 0.5*pd_thickness) ;
+
+ G4LogicalVolume* logicPD = new G4LogicalVolume(solidPhotoDiode, PDMaterial, NamePD,0,0,0);
+ logicPD->SetSensitiveDetector(m_PDScorer);
+ logicPD->SetVisAttributes(m_VisPD);
+
+ G4RotationMatrix* myRotation = new G4RotationMatrix();
+ if(horizontal) { myRotation->rotateY(90.*deg); }
+ else { myRotation->rotateX(90.*deg); }
+ //
+ // RIGHT / TOP
+ G4int side = horizontal ? TLightPipeData::Right_t : TLightPipeData::Top_t;
+ G4ThreeVector PD_Pos = horizontal ?
+ G4ThreeVector(pipe_length*0.5 + 0.5*pd_thickness, pipe_center_xy, pipe_center_z) :
+ G4ThreeVector(pipe_center_xy, pipe_length*0.5 + 0.5*pd_thickness, pipe_center_z) ;
+ new G4PVPlacement(myRotation, PD_Pos, logicPD, NamePD, world, false, pdNum, warnOverlap);
+ m_DetectorMap[pdNum++] = make_tuple(side, ilayer, irc);
+ //
+ // LEFT / BOTTOM
+ side = horizontal ? TLightPipeData::Left_t : TLightPipeData::Bottom_t;
+ if(horizontal) { PD_Pos.set(-PD_Pos.x(),+PD_Pos.y(),PD_Pos.z()); }
+ else { PD_Pos.set(+PD_Pos.x(),-PD_Pos.y(),PD_Pos.z()); }
+ new G4PVPlacement(myRotation, PD_Pos, logicPD, NamePD, world, false, pdNum, warnOverlap);
+ m_DetectorMap[pdNum++] = make_tuple(side, ilayer, irc);
+#endif
+ }
+ };
+ //
+ if(!(ilayer%2)) { // even
+ // even layers get vertical light pipes downstream
+ buildPipe(false, detZ);
+ }
+ buildPipe(true, detZ);
+ } // for(ilayer)
+#endif
+ }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
@@ -224,41 +379,249 @@ void LightPipe::InitializeRootOutput(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
-void LightPipe::ReadSensitive(const G4Event* ){
+void LightPipe::ReadSensitive(const G4Event* event){
m_Event->Clear();
-
+ return;
///////////
// Calorimeter scorer
- CalorimeterScorers::PS_Calorimeter* Scorer= (CalorimeterScorers::PS_Calorimeter*) m_LightPipeScorer->GetPrimitive(0);
+ NPS::HitsMap<G4double*>* CaloHitMap;
+ std::map<G4int, G4double**>::iterator Calo_itr;
+
+ G4int CaloCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("LightPipeScorer/Calorimeter");
+ CaloHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CaloCollectionID));
+
+ // Loop on the Calo map
+#if 0
+ for (Calo_itr = CaloHitMap->GetMap()->begin() ; Calo_itr != CaloHitMap->GetMap()->end() ; Calo_itr++){
- unsigned int size = Scorer->GetMult();
- for(unsigned int i = 0 ; i < size ; i++){
- double Energy = RandGauss::shoot(Scorer->GetEnergy(i),LightPipe_NS::ResoEnergy);
+ G4double* Info = *(Calo_itr->second);
+ double Energy = RandGauss::shoot(Info[0],LightPipe_NS::ResoEnergy);
if(Energy>LightPipe_NS::EnergyThreshold){
- double Time = RandGauss::shoot(Scorer->GetTime(i),LightPipe_NS::ResoTime);
- int DetectorNbr = Scorer->GetLevel(i)[0];
- m_Event->SetEnergy(DetectorNbr,Energy);
- m_Event->SetTime(DetectorNbr,Time);
+ //double Time = RandGauss::shoot(Info[1],LightPipe_NS::ResoTime);
+ int DetectorNbr = (int) Info[2];
+ m_Event->SetEnergy(0,0,DetectorNbr,Energy);
}
}
+#endif
+ // clear map for next event
+ CaloHitMap->clear();
+
+ // PhotoDiode //
+ NPS::HitsMap<G4double*>* PhotoDiodeHitMap;
+ std::map<G4int, G4double**>::iterator PhotoDiode_itr;
+
+ G4int PhotoDiodeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("PDScorer/PhotoDiode");
+ PhotoDiodeHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(PhotoDiodeCollectionID));
+
+ // Loop on the PhotoDiode map
+ map<int, int> NumberOfOpticalPhoton; // <det no, # photons>
+ for(const auto& hit : *PhotoDiodeHitMap->GetMap()) {
+ G4double* Info = *(hit.second);
+ G4int detectorNumber = Info[7];
+ if(NumberOfOpticalPhoton.find(detectorNumber) == NumberOfOpticalPhoton.end()){
+ NumberOfOpticalPhoton[detectorNumber] = 0;
+ }
+ NumberOfOpticalPhoton[detectorNumber] += Info[8];
+ }
+
+ for(const auto& pd : NumberOfOpticalPhoton) {
+ int det = pd.first; // detector number
+ int numPhoton = pd.second;
+ auto detMap = m_DetectorMap.find(det);
+ if(detMap != m_DetectorMap.end()){
+ const auto& side = get<0>(detMap->second);
+ const auto& layer = get<1>(detMap->second);
+ const auto& row = get<2>(detMap->second);
+ m_Event->SetEnergy(side, layer, row, numPhoton);
+ } else {
+ std::cerr << "WARNING:: Detector number encountered without map! The number is: " << det << "\nSkipping event...\n";
+ }
+ }
+
+ PhotoDiodeHitMap->clear();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////
-void LightPipe::InitializeScorers() {
+void LightPipe::InitializeScorers() {
// This check is necessary in case the geometry is reloaded
- bool already_exist = false;
+ bool already_exist = false;
+ vector<G4int> NestingLevel;
+ NestingLevel.push_back(0);
+ NestingLevel.push_back(1);
m_LightPipeScorer = CheckScorer("LightPipeScorer",already_exist) ;
+ m_PDScorer = CheckScorer("PDScorer",already_exist) ;
- if(already_exist)
+ 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) ;
+ G4VPrimitiveScorer* Calorimeter= new CalorimeterScorers::PS_Calorimeter("Calorimeter",NestingLevel, 0) ;
//and register it to the multifunctionnal detector
m_LightPipeScorer->RegisterPrimitive(Calorimeter);
+
+ G4VPrimitiveScorer* PDScorer =
+ new PHOTODIODESCORERS::PS_PhotoDiode_Rectangle
+ ("PhotoDiode",0,GetPipeWidth(0),GetPipeWidth(0),1,1);
+ m_PDScorer->RegisterPrimitive(PDScorer);
+
G4SDManager::GetSDMpointer()->AddNewDetector(m_LightPipeScorer) ;
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_PDScorer) ;
+}
+
+G4Material* LightPipe::CreateScintillatorMaterial() const {
+ // p-Terphenyl
+ // Taken from Proteus, Inc. specs
+ // See here: http://people.physics.tamu.edu/christian/files/p-terphenyl.png
+ //
+ G4double specificGravity = 1.23; // from proteus
+ G4double densityReference = 0.999972*g/cm3; // water @4-deg C [ the standard ]
+ G4double density = specificGravity*densityReference;
+
+ G4Material* material = new G4Material("p_Terphenyl_Scint", density, 2);
+ material->AddElement(MaterialManager::getInstance()->GetElementFromLibrary("H"), 10);
+ material->AddElement(MaterialManager::getInstance()->GetElementFromLibrary("C"), 14);
+
+ // Adding Scintillation property:
+ vector<double> energy, scint, fast, slow, rindx, atten;
+ {
+ // Read emission spectrum from datfile (+ add constant parameters)
+ ifstream ifs(gSystem->ExpandPathName("$NPTOOL/NPSimulation/Detectors/LightPipe/p-terphenyl_emission.dat"));
+ if(!ifs.good()){
+ std::cerr << "ERROR: Couldn't open file: \"$NPTOOL/NPSimulation/Detectors/LightPipe/p-terphenyl_emission.dat\"\n";
+ exit(1);
+ }
+ double wl, pr;
+ while(ifs >> wl >> pr) {
+ energy.emplace_back( h_Planck*c_light / (wl*nm) ); // convert to energy
+ scint.emplace_back(pr); // scintillation probability
+ rindx.emplace_back(1.65); // refractive index
+ fast.emplace_back(3*ns); // FAST component ???
+ slow.emplace_back(100*ns); // SLOW component ???
+ atten.emplace_back(4.73*mm); // Attenuation length (from https://arxiv.org/pdf/1305.0442.pdf)
+ }
+ }
+
+ // Set Material Properties
+ G4int numPoints = energy.size();
+ G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
+ MPT->AddConstProperty("SCINTILLATIONYIELD", 27000/MeV); // from proteus
+ MPT->AddProperty("SCINTILLATION", &energy[0], &scint[0], numPoints);
+ MPT->AddProperty("RINDEX", &energy[0], &rindx[0], numPoints) ;
+ MPT->AddProperty("ABSLENGTH", &energy[0], &atten[0], numPoints);
+ MPT->AddProperty("FASTCOMPONENT", &energy[0], &fast[0], numPoints);
+ MPT->AddProperty("SLOWCOMPONENT", &energy[0], &slow[0], numPoints);
+ MPT->AddConstProperty("RESOLUTIONSCALE", 1.0);
+ MPT->AddConstProperty("FASTTIMECONSTANT", 20*ns); // ?????
+ MPT->AddConstProperty("SLOWTIMECONSTANT", 100*ns); // ?????
+ MPT->AddConstProperty("YIELDRATIO",1.0); // ?????
+ material->SetMaterialPropertiesTable(MPT);
+ return material;
+}
+
+G4Material* LightPipe::CreatePipeMaterial() const {
+ // Bicron BC-482A
+ // https://www.crystals.saint-gobain.com/sites/imdf.crystals.com/files/documents/bc482a-bc484-data-sheet.pdf
+ // https://www.crystals.saint-gobain.com/sites/imdf.crystals.com/files/documents/organics-plastic-scintillators.pdf
+ //
+ G4double density = 1.03*g/cm3;
+
+ G4Material* material = new G4Material("BC482A_WLS", density, 2);
+ material->AddElement(MaterialManager::getInstance()->GetElementFromLibrary("H"), 10); // H:C ratio - 1.110
+ material->AddElement(MaterialManager::getInstance()->GetElementFromLibrary("C"), 9);
+
+ // Adding WLS property
+ auto readDatfile = [](const char* fname, vector<double>& energy, vector<double>& abs, vector<double>& emit){
+ ifstream ifs(gSystem->ExpandPathName(fname));
+ if(!ifs.good()) {
+ std::cerr << "ERROR: no file: \"" << fname << "\"\n";
+ exit(1);
+ }
+ // skip header
+ std::string dummy;
+ std::getline(ifs, dummy);
+
+ double wl, pr_a, pr_e;
+ while(ifs >> wl >> pr_a >> pr_e) {
+ energy.emplace_back( h_Planck*c_light / (wl*nm) ); // convert to energy
+ abs.emplace_back(pr_a); // absorption probability
+ emit.emplace_back(pr_e); // emission probability
+ }
+ // SORT IN ORDER OF INCREASING ENERGY
+ vector<int> isort(energy.size());
+ vector<double> e0 = energy, pa0 = abs, pe0 = emit;
+ TMath::Sort((int)energy.size(), &energy[0], &isort[0]);
+ for(size_t i=0; i< energy.size(); ++i){
+ energy.at(i) = e0.at(isort.at(i));
+ abs.at(i) = pa0.at(isort.at(i));
+ emit.at(i) = pe0.at(isort.at(i));
+ }
+ };
+ // Absorption & Emission
+ vector<double> energy, p_abs, p_emit;
+ readDatfile("$NPTOOL/NPSimulation/Detectors/LightPipe/BC482A_properties.dat", energy, p_abs, p_emit);
+ //
+ // Absorption is given as a probability, but GEANT4 needs a length
+ // Invert and set the minimum to 0.4336 mm, which is the (measured)
+ // attenuation length for EJ-280 @peak absorption.
+ // This is not exact, but it's close.
+ // For absorption of 0, set attenuation length very long (5 m)
+ for(auto&& p : p_abs) {
+ p = p > 0 ? (0.4336*mm) / p : 5*m;
+ }
+
+ //
+ const size_t numPoints = energy.size();
+ vector<double> rindx(numPoints, 1.59);
+ vector<double> abslength(numPoints, 400*cm); // BULK attenuation length
+
+ // Set Material Properties
+ G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
+ MPT->AddProperty("RINDEX", &energy[0], &rindx[0], numPoints);
+ MPT->AddProperty("ABSLENGTH", &energy[0], &abslength[0], numPoints);
+ MPT->AddProperty("WLSABSLENGTH", &energy[0], &p_abs[0], numPoints);
+ MPT->AddProperty("WLSCOMPONENT", &energy[0], &p_emit[0], numPoints);
+ MPT->AddConstProperty("WLSTIMECONSTANT", 12.*ns);
+ MPT->AddConstProperty("WLSMEANNUMBEROFPHOTONS", 0.86);
+
+ material->SetMaterialPropertiesTable(MPT);
+ return material;
+}
+
+G4Material* LightPipe::CreateWrappingMaterial() const {
+ // Teflon (C2F4) -- for now
+ //
+ G4double density = 2.2*g/cm3;
+
+ G4Material* material = new G4Material("TEFLON", density, 2);
+ material->AddElement(MaterialManager::getInstance()->GetElementFromLibrary("F"), 4); // H:C ratio - 1.110
+ material->AddElement(MaterialManager::getInstance()->GetElementFromLibrary("C"), 2);
+ return material;
+}
+
+G4OpticalSurface* LightPipe::CreateReflectiveSurface() const {
+ G4OpticalSurface* OpticalCrystalSurface = new G4OpticalSurface("CrystalSurface");
+ OpticalCrystalSurface->SetType(dielectric_metal);
+ //polished: smooth perfectly polished surcface
+ //ground: rough surface
+ OpticalCrystalSurface->SetFinish(polished);
+ //unified
+ //glisur
+ OpticalCrystalSurface->SetModel(glisur);
+
+ G4double pp[] = {0.01*eV, 10*eV};
+ const G4int num = sizeof(pp)/sizeof(G4double);
+ G4double reflectivity[] = {1., 1.};
+ G4double efficiency[] = {1., 1.};
+
+ G4MaterialPropertiesTable* OpticalCrystalSurfaceProperty = new G4MaterialPropertiesTable();
+
+ OpticalCrystalSurfaceProperty->AddProperty("REFLECTIVITY",pp,reflectivity,num);
+ OpticalCrystalSurfaceProperty->AddProperty("EFFICIENCY",pp,efficiency,num);
+ OpticalCrystalSurface->SetMaterialPropertiesTable(OpticalCrystalSurfaceProperty);
+ return OpticalCrystalSurface;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -276,13 +639,13 @@ NPS::VDetector* LightPipe::Construct(){
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern"C" {
- class proxy_nps_LightPipe{
- public:
- proxy_nps_LightPipe(){
- NPS::DetectorFactory::getInstance()->AddToken("LightPipe","LightPipe");
- NPS::DetectorFactory::getInstance()->AddDetector("LightPipe",LightPipe::Construct);
- }
- };
-
- proxy_nps_LightPipe p_nps_LightPipe;
+class proxy_nps_LightPipe{
+public:
+ proxy_nps_LightPipe(){
+ NPS::DetectorFactory::getInstance()->AddToken("LightPipe","LightPipe");
+ NPS::DetectorFactory::getInstance()->AddDetector("LightPipe",LightPipe::Construct);
+ }
+};
+
+proxy_nps_LightPipe p_nps_LightPipe;
}
diff --git a/NPSimulation/Detectors/LightPipe/LightPipe.hh b/NPSimulation/Detectors/LightPipe/LightPipe.hh
index 72d7eb87c156c1e6df57527b80d8654f15a7f22f..9980c8d220a36915901e9bfd29a84dd4259f0ae5 100644
--- a/NPSimulation/Detectors/LightPipe/LightPipe.hh
+++ b/NPSimulation/Detectors/LightPipe/LightPipe.hh
@@ -24,7 +24,8 @@
// C++ header
#include <string>
#include <vector>
-using namespace std;
+#include <tuple>
+#include <map>
// G4 headers
#include "G4ThreeVector.hh"
@@ -37,81 +38,91 @@ using namespace std;
#include "TLightPipeData.h"
#include "NPInputParser.h"
+class G4OpticalSurface;
+
class LightPipe : public NPS::VDetector{
////////////////////////////////////////////////////
/////// Default Constructor and Destructor /////////
////////////////////////////////////////////////////
- public:
- LightPipe() ;
- virtual ~LightPipe() ;
-
- ////////////////////////////////////////////////////
- /////// Specific Function of this Class ///////////
- ////////////////////////////////////////////////////
- public:
- // Cartesian
- void AddDetector(G4ThreeVector POS, string Shape);
- // Spherical
- void AddDetector(double R,double Theta,double Phi,string Shape);
-
+public:
+ LightPipe() ;
+ virtual ~LightPipe() ;
- G4LogicalVolume* BuildSquareDetector();
- G4LogicalVolume* BuildCylindricalDetector();
-
- private:
- G4LogicalVolume* m_SquareDetector;
- G4LogicalVolume* m_CylindricalDetector;
+ ////////////////////////////////////////////////////
+ /////// Specific Function of this Class ///////////
+ ////////////////////////////////////////////////////
+public:
+ // Cartesian
+ void AddDetector(G4int nrow, G4int ncol, G4int nlayer, G4double width, G4double thickness, G4double pipe_width, G4double pipe_thickness);
- ////////////////////////////////////////////////////
- ////// Inherite from NPS::VDetector class /////////
- ////////////////////////////////////////////////////
- public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- // Called in DetecorConstruction::ReadDetextorConfiguration Method
- void ReadConfiguration(NPL::InputParser) ;
+ ////////////////////////////////////////////////////
+ ////// Inherite from NPS::VDetector class /////////
+ ////////////////////////////////////////////////////
+public:
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(NPL::InputParser) ;
- // Construct detector and inialise sensitive part.
- // Called After DetecorConstruction::AddDetector Method
- void ConstructDetector(G4LogicalVolume* world) ;
+ // Construct detector and inialise sensitive part.
+ // Called After DetecorConstruction::AddDetector Method
+ void ConstructDetector(G4LogicalVolume* world) ;
- // Add Detector branch to the EventTree.
- // Called After DetecorConstruction::AddDetector Method
- void InitializeRootOutput() ;
+ // Add Detector branch to the EventTree.
+ // Called After DetecorConstruction::AddDetector Method
+ void InitializeRootOutput() ;
- // Read sensitive part and fill the Root tree.
- // Called at in the EventAction::EndOfEventAvtion
- void ReadSensitive(const G4Event* event) ;
+ // Read sensitive part and fill the Root tree.
+ // Called at in the EventAction::EndOfEventAvtion
+ void ReadSensitive(const G4Event* event) ;
- public: // Scorer
- // Initialize all Scorer used by the MUST2Array
- void InitializeScorers() ;
+private:
+ G4Material* m_ScintillatorMaterial;
+ G4Material* m_PipeMaterial;
+ G4Material* m_WrappingMaterial;
+ G4OpticalSurface* m_ReflectiveSurface;
- // Associated Scorer
- G4MultiFunctionalDetector* m_LightPipeScorer ;
- ////////////////////////////////////////////////////
- ///////////Event class to store Data////////////////
- ////////////////////////////////////////////////////
- private:
- TLightPipeData* m_Event ;
+ G4Material* CreateScintillatorMaterial() const;
+ G4Material* CreatePipeMaterial() const;
+ G4Material* CreateWrappingMaterial() const;
+ G4OpticalSurface* CreateReflectiveSurface() const;
+
+public: // Scorer
+ // Initialize all Scorer used by the MUST2Array
+ void InitializeScorers() ;
- ////////////////////////////////////////////////////
- ///////////////Private intern Data//////////////////
- ////////////////////////////////////////////////////
- private: // Geometry
- // Detector Coordinate
- vector<double> m_R;
- vector<double> m_Theta;
- vector<double> m_Phi;
-
- // Shape type
- vector<string> m_Shape ;
-
- // Visualisation Attribute
- G4VisAttributes* m_VisSquare;
- G4VisAttributes* m_VisCylinder;
+ // Associated Scorer
+ G4MultiFunctionalDetector* m_LightPipeScorer ;
+ G4MultiFunctionalDetector* m_PDScorer;
+ ////////////////////////////////////////////////////
+ ///////////Event class to store Data////////////////
+ ////////////////////////////////////////////////////
+private:
+ TLightPipeData* m_Event ;
+ std::map<int, std::tuple<int, int, int> > m_DetectorMap; // <detNo, <side, layer, row/col> >
+
+public:
+ G4int GetNrow(size_t i) const { return std::get<0>(m_Detector.at(i)); }
+ G4int GetNcol(size_t i) const { return std::get<1>(m_Detector.at(i)); }
+ G4int GetNlayer(size_t i) const { return std::get<2>(m_Detector.at(i)); }
+ G4double GetWidth(size_t i) const { return std::get<3>(m_Detector.at(i)); }
+ G4double GetThickness(size_t i) const { return std::get<4>(m_Detector.at(i)); }
+ G4double GetPipeWidth(size_t i) const { return std::get<5>(m_Detector.at(i)); }
+ G4double GetPipeThickness(size_t i) const { return std::get<6>(m_Detector.at(i)); }
+ ////////////////////////////////////////////////////
+ ///////////////Private intern Data//////////////////
+ ////////////////////////////////////////////////////
+private: // Geometry
+ // Detector Coordinate
+ std::vector<std::tuple<G4int, G4int, G4int, G4double, G4double, G4double, G4double> > m_Detector;
+
+private:
+ // Visualisation Attribute
+ G4VisAttributes* m_VisSquare;
+ G4VisAttributes* m_VisPipe;
+ G4VisAttributes* m_VisPD;
// Needed for dynamic loading of the library
- public:
- static NPS::VDetector* Construct();
+public:
+ static NPS::VDetector* Construct();
};
#endif
diff --git a/NPSimulation/Detectors/LightPipe/WLS_absorption.dat b/NPSimulation/Detectors/LightPipe/WLS_absorption.dat
new file mode 100644
index 0000000000000000000000000000000000000000..ea2dc5dbd125b6154fd71f5bbc8dbdaea8977a63
--- /dev/null
+++ b/NPSimulation/Detectors/LightPipe/WLS_absorption.dat
@@ -0,0 +1,46 @@
+324 0.23543424030519855
+330 0.2797744973266716
+336 0.3568982562212688
+342 0.45794399557184606
+348 0.5720861560960396
+354 0.6423607325337347
+360 0.7237046979297348
+366 0.8169704975047025
+372 0.8954133588231912
+378 0.9457181027326789
+384 0.9901303540388469
+390 0.9955050870093485
+396 0.9538469297065142
+402 0.8898202898772857
+408 0.8584075622446017
+414 0.8364139962860144
+420 0.7851614701046612
+426 0.7082888913779626
+432 0.6392903396510657
+438 0.5515063263380027
+444 0.4256967788043867
+450 0.2415679335503349
+456 0.13310591472926148
+462 0.030562155168096572
+468 0.008148621487096941
+474 0.002747241342357043
+480 0.
+486 0.
+492 0.
+498 0.
+504 0.
+510 0.
+516 0.
+522 0.
+528 0.
+534 0.
+540 0.
+546 0.
+552 0.
+558 0.
+564 0.
+570 0.
+576 0.
+582 0.
+588 0.
+594 0.
diff --git a/NPSimulation/Detectors/LightPipe/WLS_emission.dat b/NPSimulation/Detectors/LightPipe/WLS_emission.dat
new file mode 100644
index 0000000000000000000000000000000000000000..5ead863e4071735e0520bac20d13c7104d18f96f
--- /dev/null
+++ b/NPSimulation/Detectors/LightPipe/WLS_emission.dat
@@ -0,0 +1,46 @@
+324 0.
+330 0.
+336 0.
+342 0.
+348 0.
+354 0.
+360 0.
+366 0.
+372 0.
+378 0.
+384 0.
+390 0.
+396 0.
+402 0.
+408 0.
+414 0.
+420 0.
+426 0.
+432 0.
+438 0.
+444 0.08657197078612922
+450 0.24244440143865642
+456 0.463427444632436
+462 0.8093744108256966
+468 0.9283384225343975
+474 0.9789958445998295
+480 0.9952580190378026
+486 0.9683221897945588
+492 0.8580880346078893
+498 0.695394404265841
+504 0.5420850105165227
+510 0.44045761435208164
+516 0.37538391688413775
+522 0.32514287498636296
+528 0.2727348311625266
+534 0.21490510479855662
+540 0.16099793158878195
+546 0.11888268966239912
+552 0.08963478399179192
+558 0.06809158627158829
+564 0.0527319299208866
+570 0.03171523603459048
+576 0.020248021089946944
+582 0.014252509594427476
+588 0.010999675726623082
+594 0.008252310074156277
diff --git a/NPSimulation/Detectors/LightPipe/WLS_properties.dat b/NPSimulation/Detectors/LightPipe/WLS_properties.dat
new file mode 100644
index 0000000000000000000000000000000000000000..d1993d1bbf5160ca0fde3bd40b9c318efc34832d
--- /dev/null
+++ b/NPSimulation/Detectors/LightPipe/WLS_properties.dat
@@ -0,0 +1,47 @@
+Wavelength:Absorption:Emission/D
+324 0.23543424030519855 0.
+330 0.2797744973266716 0.
+336 0.3568982562212688 0.
+342 0.45794399557184606 0.
+348 0.5720861560960396 0.
+354 0.6423607325337347 0.
+360 0.7237046979297348 0.
+366 0.8169704975047025 0.
+372 0.8954133588231912 0.
+378 0.9457181027326789 0.
+384 0.9901303540388469 0.
+390 0.9955050870093485 0.
+396 0.9538469297065142 0.
+402 0.8898202898772857 0.
+408 0.8584075622446017 0.
+414 0.8364139962860144 0.
+420 0.7851614701046612 0.
+426 0.7082888913779626 0.
+432 0.6392903396510657 0.
+438 0.5515063263380027 0.
+444 0.4256967788043867 0.08657197078612922
+450 0.2415679335503349 0.24244440143865642
+456 0.13310591472926148 0.463427444632436
+462 0.030562155168096572 0.8093744108256966
+468 0.008148621487096941 0.9283384225343975
+474 0.002747241342357043 0.9789958445998295
+480 0. 0.9952580190378026
+486 0. 0.9683221897945588
+492 0. 0.8580880346078893
+498 0. 0.695394404265841
+504 0. 0.5420850105165227
+510 0. 0.44045761435208164
+516 0. 0.37538391688413775
+522 0. 0.32514287498636296
+528 0. 0.2727348311625266
+534 0. 0.21490510479855662
+540 0. 0.16099793158878195
+546 0. 0.11888268966239912
+552 0. 0.08963478399179192
+558 0. 0.06809158627158829
+564 0. 0.0527319299208866
+570 0. 0.03171523603459048
+576 0. 0.020248021089946944
+582 0. 0.014252509594427476
+588 0. 0.010999675726623082
+594 0. 0.008252310074156277
diff --git a/NPSimulation/Detectors/LightPipe/p-terphenyl_emission.dat b/NPSimulation/Detectors/LightPipe/p-terphenyl_emission.dat
new file mode 100644
index 0000000000000000000000000000000000000000..911fd86e3b592819efa3fa5f6e107bd40df7c4b4
--- /dev/null
+++ b/NPSimulation/Detectors/LightPipe/p-terphenyl_emission.dat
@@ -0,0 +1,75 @@
+361.445321354552 1.444922829859e-02
+363.533013158026 1.992315208414e-02
+364.607099230735 2.869952020107e-02
+365.934586736136 4.028891270207e-02
+366.827511784590 8.985356445997e-02
+366.261348583594 6.378451194015e-02
+365.727105563070 5.495775945116e-02
+367.906637885795 1.267952419141e-01
+368.982123966420 1.732616789934e-01
+370.072450130948 2.544713850373e-01
+371.160256281229 3.635928111621e-01
+372.695112959061 6.233185502033e-01
+371.718019434732 4.647609552326e-01
+376.097990865051 8.630744205872e-01
+376.469086296499 1.000000000000e+00
+380.495509061244 7.496087058979e-01
+382.834082283160 5.997447928985e-01
+384.493651666098 4.822893537483e-01
+388.027551646199 3.892674207881e-01
+391.648532118639 4.936935058515e-01
+392.715898153355 6.229178224304e-01
+394.292307066117 7.753631725814e-01
+397.402564651012 9.681631000138e-01
+401.972470488523 8.510528339580e-01
+404.372028055235 6.591217885477e-01
+405.966693071214 5.398317811562e-01
+408.375826692067 4.341851531419e-01
+412.195608288497 3.754834092222e-01
+415.889669174122 4.479017417278e-01
+420.794456905001 5.864537571663e-01
+425.369682772591 5.261713873970e-01
+428.389779847734 4.187470504375e-01
+430.391511165200 3.370199593801e-01
+432.729244382367 2.634128119145e-01
+435.922942439021 2.100928551832e-01
+441.537870184918 2.063859494648e-01
+447.182926101154 2.396809152272e-01
+452.266074840458 2.095415871399e-01
+455.798014809477 1.641183146772e-01
+457.388423801394 1.272126261154e-01
+460.140279359942 1.002024376139e-01
+463.331457402348 7.687915721502e-02
+468.408166105241 6.350386965489e-02
+474.543560793774 6.616764087099e-02
+480.143816456717 5.912486759856e-02
+484.190175334178 4.587216052248e-02
+487.721835301614 3.443216612122e-02
+491.759738131266 2.484442744244e-02
+496.309539855112 1.787785615624e-02
+501.894115429403 1.402851166410e-02
+508.004869978624 1.227757393170e-02
+514.108344486685 1.015310425803e-02
+519.689224040079 7.631673360851e-03
+525.266407572576 5.522321716898e-03
+530.837767072145 3.765871563834e-03
+536.425030661633 2.901162607052e-03
+542.543625255181 2.656083660074e-03
+548.662779851895 2.441689797166e-03
+554.776894420114 2.145581012474e-03
+560.883728947172 1.763139708609e-03
+567.004003550218 1.631609711317e-03
+573.125398159597 1.524506269596e-03
+579.255192816468 1.537699988162e-03
+585.381627454342 1.504113933399e-03
+591.514782130211 1.564396165356e-03
+597.640656764919 1.522476355408e-03
+603.771571428122 1.551429813673e-03
+609.908086122988 1.663345549497e-03
+616.039560789357 1.703296236145e-03
+622.174955477890 1.806729639767e-03
+628.309230160091 1.896706428831e-03
+634.442944839125 1.980801015043e-03
+640.574419505495 2.027543975111e-03
+646.709254190862 2.137886140935e-03
+650.970318282232 2.232582668682e-03
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index b7cf24e0dad55050812cfa12de1bccc001f2f7d5..27158ff4ecf1c7916123d31012d9cab47c6b33b2 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -184,10 +184,17 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep)
G4ParticleDefinition* LightName;
- if(m_Reaction.GetUseExInGeant4())
- LightName = IonTable->GetIon(LightZ, LightA, m_Reaction.GetExcitation3()*MeV);
- else
- LightName = IonTable->GetIon(LightZ, LightA);
+ if(LightZ == 0 && LightA == 1) // neutron is special case
+ {
+ LightName = G4Neutron::Definition();
+ }
+ else
+ {
+ if(m_Reaction.GetUseExInGeant4())
+ LightName = IonTable->GetIon(LightZ, LightA, m_Reaction.GetExcitation3()*MeV);
+ else
+ LightName = IonTable->GetIon(LightZ, LightA);
+ }
// Nucleus 4
G4int HeavyZ = m_Reaction.GetNucleus4()->GetZ() ;