diff --git a/Examples/Example1/energy_loss/.gitignore b/Examples/Example1/energy_loss/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Examples/Example1/project.config b/Examples/Example1/project.config
new file mode 100644
index 0000000000000000000000000000000000000000..151c3cc51df56fac51e913bfa50e2f415d4563fc
--- /dev/null
+++ b/Examples/Example1/project.config
@@ -0,0 +1,5 @@
+Project Example1
+ AnalysisOutput= ./root/analysis
+ SimulationOutput= ./root/simulation
+ EnergyLoss= ./energy_loss
+
diff --git a/Examples/Example1/root/analysis/.gitignore b/Examples/Example1/root/analysis/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Examples/Example1/root/simulation/.gitignore b/Examples/Example1/root/simulation/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Examples/Example4/run.mac b/Examples/Example4/run.mac
index ced78f72769fe42cd7cc053d0387e39f4e254d76..19fd56e6724ca337abf96c30b34f604643ab526b 100644
--- a/Examples/Example4/run.mac
+++ b/Examples/Example4/run.mac
@@ -1 +1 @@
-/run/beamOn 10
+/run/beamOn 1
diff --git a/NPLib/Core/NPDetectorManager.cxx b/NPLib/Core/NPDetectorManager.cxx
index 6805986b6aa34c48fa9b93a3a380f605f1732cf4..cfb35ff5d874d55ba0394b64078690683dd6067a 100644
--- a/NPLib/Core/NPDetectorManager.cxx
+++ b/NPLib/Core/NPDetectorManager.cxx
@@ -634,7 +634,6 @@ void NPL::DetectorManager::CheckSpectraServer(){
////////////////////////////////////////////////////////////////////////////////
void NPL::DetectorManager::FillOutputTree(){
- //std::cout << "Root Output " << m_RootOutput << std::endl;
m_RootOutput->Fill();
}
diff --git a/NPLib/Core/NPOptionManager.cxx b/NPLib/Core/NPOptionManager.cxx
index deec47209103fce307bf7af0d8dabbad412a8f90..9b9f6c3c408df7e8aad222fbc283cec8ba5cc41d 100644
--- a/NPLib/Core/NPOptionManager.cxx
+++ b/NPLib/Core/NPOptionManager.cxx
@@ -71,6 +71,9 @@ void NPOptionManager::ReadProjectConfigFile(){
if(blocks[i]->HasToken("EnergyLoss"))
m_EnergyLossPath = blocks[i]->GetString("EnergyLoss");
+
+ if(blocks[i]->HasToken("Cuts"))
+ m_CutsPath = blocks[i]->GetString("Cuts");
}
}
diff --git a/NPLib/Core/NPOptionManager.h b/NPLib/Core/NPOptionManager.h
index 30c6283a8ecd274ba9b19bab22b836b6a2a9638e..d37abeca53aaa88d782c01aed469c846e9886278 100644
--- a/NPLib/Core/NPOptionManager.h
+++ b/NPLib/Core/NPOptionManager.h
@@ -125,6 +125,7 @@ class NPOptionManager{
std::string GetSimulationOutputPath(){return m_SimulationOutputPath;};
std::string GetCalibrationOutputPath(){return m_CalibrationOutputPath;};
std::string GetEnergyLossPath(){return m_EnergyLossPath;};
+ std::string GetCutsPath(){return m_CutsPath;};
// Setters
void SetReactionFile(const std::string& name) {fReactionFileName = name;CheckEventGenerator();}
void SetDetectorFile(const std::string& name) {fDetectorFileName = name;CheckDetectorConfiguration();}
@@ -187,6 +188,7 @@ class NPOptionManager{
std::string m_SimulationOutputPath;// output path of simulated tree
std::string m_CalibrationOutputPath;// output path of calibration histograms
std::string m_EnergyLossPath;// input/output path of energy loss table
+ std::string m_CutsPath;// CutsPath
};
diff --git a/NPLib/Core/RootHistogramsCalib.h b/NPLib/Core/RootHistogramsCalib.h
index 1dab362cd3c91ee073fe5830cfe7e57fee18356f..11de2783aa161ab86d678d0978ddf7ebb3e860f5 100644
--- a/NPLib/Core/RootHistogramsCalib.h
+++ b/NPLib/Core/RootHistogramsCalib.h
@@ -33,6 +33,7 @@
#include "TGraphErrors.h"
#include "TH1F.h"
#include "TH2F.h"
+#include "TCutG.h"
#include "TFile.h"
#include "TTree.h"
#include "TList.h"
@@ -75,6 +76,8 @@ public:
std::map<TString,std::map<TString,TH1F*>>* GetTH1Map(){return TH1Map;};
std::map<TString,std::map<TString,TH2F*>>* GetTH2Map(){return TH2Map;};
std::map<TString,std::map<TString,TGraphErrors*>>* GetTGraphMap(){return TGraphMap;};
+ std::map<TString,std::map<TString,TCutG*>>* GetTCutGMap(){return TCutGMap;};
+ std::map<TString,std::map<TString,TFile*>>* GetTFileMap(){return TFileMap;};
void Fill();
private:
@@ -85,6 +88,8 @@ private:
std::map<TString,std::map<TString,TH1F*>>* TH1Map = new std::map<TString,std::map<TString,TH1F*>>;
std::map<TString,std::map<TString,TH2F*>>* TH2Map = new std::map<TString,std::map<TString,TH2F*>>;
std::map<TString,std::map<TString,TGraphErrors*>>* TGraphMap = new std::map<TString,std::map<TString,TGraphErrors*>>;
+ std::map<TString,std::map<TString,TCutG*>>* TCutGMap = new std::map<TString,std::map<TString,TCutG*>>;
+ std::map<TString,std::map<TString,TFile*>>* TFileMap = new std::map<TString,std::map<TString,TFile*>>;
};
#endif // ROOTHISTOCALIB_HH
diff --git a/NPLib/Detectors/Catana/TCatanaPhysics.h b/NPLib/Detectors/Catana/TCatanaPhysics.h
index f9fc98d14c71bc3b39a9c9502031cf4d4c9742fd..a238c294fba4493861614b134460aef56450b022 100644
--- a/NPLib/Detectors/Catana/TCatanaPhysics.h
+++ b/NPLib/Detectors/Catana/TCatanaPhysics.h
@@ -18,174 +18,177 @@
* *
*---------------------------------------------------------------------------*
* Comment: *
- * *
+ * *
* *
*****************************************************************************/
-// C++ headers
-#include <vector>
+// C++ headers
#include <map>
#include <string>
+#include <vector>
using namespace std;
// ROOT headers
-#include "TObject.h"
#include "TH1.h"
+#include "TObject.h"
#include "TVector3.h"
// NPTool headers
-#include "TCatanaData.h"
-#include "TCatanaSpectra.h"
#include "NPCalibrationManager.h"
-#include "NPVDetector.h"
#include "NPInputParser.h"
+#include "NPVDetector.h"
+#include "TCatanaData.h"
+#include "TCatanaSpectra.h"
// forward declaration
class TCatanaSpectra;
-
-
class TCatanaPhysics : public TObject, public NPL::VDetector {
//////////////////////////////////////////////////////////////
// constructor and destructor
- public:
- TCatanaPhysics();
- ~TCatanaPhysics() {};
-
+ public:
+ TCatanaPhysics();
+ ~TCatanaPhysics(){};
//////////////////////////////////////////////////////////////
// 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> DetectorNumber;
+ vector<double> Energy;
+ vector<double> Time;
/// A usefull method to bundle all operation to add a detector
- void AddDetector(double X, double Y, double Z, double Theta, double Phi, int ID, int Type);
- void ReadCSV(string path);
-
+ void AddDetector(double X, double Y, double Z, double Theta, double Phi, int ID, int Type);
+ void ReadCSV(string path);
+
// Position method and variable
- public:
- map<int,TVector3> m_Position;//!
- map<int,double> m_Theta;//!
- map<int,double> m_Phi;//!
- map<int,int> m_Type;//!
- TVector3 m_Ref;//!
- TVector3 GetPositionOfInteraction(int& i);//!
- // Return index of the closest hit to line defined by v1 and v2
- unsigned int FindClosestHitToLine(const TVector3& v1, const TVector3& v2, double& d);
+ public:
+ map<int, TVector3> m_Position; //!
+ map<int, double> m_Theta; //!
+ map<int, double> m_Phi; //!
+ map<int, int> m_Type; //!
+ TVector3 m_Ref; //!
+ TVector3 GetPositionOfInteraction(int& i); //!
+ // Return index of the closest hit to line defined by v1 and v2
+ unsigned int FindClosestHitToLine(const TVector3& v1, const TVector3& v2, double& d);
//////////////////////////////////////////////////////////////
// methods inherited from the VDetector ABC class
- public:
- // read stream from ConfigFile to pick-up detector parameters
- void ReadConfiguration(NPL::InputParser);
-
- // add parameters to the CalibrationManger
- void AddParameterToCalibrationManager();
+ public:
+ // read stream from ConfigFile to pick-up detector parameters
+ void ReadConfiguration(NPL::InputParser);
- // method called event by event, aiming at extracting the
- // physical information from detector
- void BuildPhysicalEvent();
+ // add parameters to the CalibrationManger
+ void AddParameterToCalibrationManager();
- // same as BuildPhysicalEvent() method but with a simpler
- // treatment
- void BuildSimplePhysicalEvent();
+ // method called event by event, aiming at extracting the
+ // physical information from detector
+ void BuildPhysicalEvent();
- // same as above but for online analysis
- void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
+ // same as BuildPhysicalEvent() method but with a simpler
+ // treatment
+ void BuildSimplePhysicalEvent();
- // 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();
+ // same as above but for online analysis
+ void BuildOnlinePhysicalEvent() { BuildPhysicalEvent(); };
- // 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 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();
- // create branches of output ROOT file
- void InitializeRootOutput();
+ // 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();
- // clear the raw and physical data objects event by event
- void ClearEventPhysics() {Clear();}
- void ClearEventData() {m_EventData->Clear();}
+ // create branches of output ROOT file
+ void InitializeRootOutput();
- // methods related to the TCatanaSpectra class
- // instantiate the TCatanaSpectra class and
- // declare list of histograms
- void InitSpectra();
+ // clear the raw and physical data objects event by event
+ void ClearEventPhysics() { Clear(); }
+ void ClearEventData() { m_EventData->Clear(); }
- // fill the spectra
- void FillSpectra();
+ // methods related to the TCatanaSpectra class
+ // instantiate the TCatanaSpectra class and
+ // declare list of histograms
+ void InitSpectra();
- // used for Online mainly, sanity check for histograms and
- // change their color if issues are found, for example
- void CheckSpectra();
+ // fill the spectra
+ void FillSpectra();
- // used for Online only, clear all the spectra
- void ClearSpectra();
+ // used for Online mainly, sanity check for histograms and
+ // change their color if issues are found, for example
+ void CheckSpectra();
- // write spectra to ROOT output file
- void WriteSpectra();
+ // used for Online only, clear all the spectra
+ void ClearSpectra();
+ // write spectra to ROOT output file
+ void WriteSpectra();
//////////////////////////////////////////////////////////////
// specific methods to Catana 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(); } //!
+
+ // read the user configuration file. If no file is found, load standard one
+ void ReadAnalysisConfig(); //!
- // clear the pre-treated object
- void ClearPreTreatedData() {m_PreTreatedData->Clear();}
+ // give and external TCatanaData object to TCatanaPhysics.
+ // needed for online analysis for example
+ void SetRawDataPointer(TCatanaData* rawDataPointer) { m_EventData = rawDataPointer; } //!
- // read the user configuration file. If no file is found, load standard one
- void ReadAnalysisConfig();
+ double GetTotalEnergy() {
+ double res = 0;
+ for (auto e : Energy) {
+ res += e;
+ }
+ return res;
+ };
- // give and external TCatanaData object to TCatanaPhysics.
- // needed for online analysis for example
- void SetRawDataPointer(TCatanaData* rawDataPointer) {m_EventData = rawDataPointer;}
-
// objects are not written in the TTree
- private:
- TCatanaData* m_EventData; //!
- TCatanaData* m_PreTreatedData; //!
- TCatanaPhysics* m_EventPhysics; //!
+ private:
+ TCatanaData* m_EventData; //!
+ TCatanaData* m_PreTreatedData; //!
+ TCatanaPhysics* m_EventPhysics; //!
// getters for raw and pre-treated data object
- public:
- TCatanaData* GetRawData() const {return m_EventData;}
- TCatanaData* GetPreTreatedData() const {return m_PreTreatedData;}
+ public:
+ TCatanaData* GetRawData() const { return m_EventData; }
+ TCatanaData* 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:
- TCatanaSpectra* m_Spectra; // !
+ private:
+ TCatanaSpectra* 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(TCatanaPhysics,1) // CatanaPhysics structure
+ ClassDef(TCatanaPhysics, 1) // CatanaPhysics structure
};
#endif
diff --git a/NPLib/Detectors/MUST2/CMakeLists.txt b/NPLib/Detectors/MUST2/CMakeLists.txt
index c674ad17c3f488d42371ddb9c6e4aefac98c773d..54e15a86f0bf9d7703f127f58d650edce4c99879 100644
--- a/NPLib/Detectors/MUST2/CMakeLists.txt
+++ b/NPLib/Detectors/MUST2/CMakeLists.txt
@@ -3,6 +3,6 @@ add_custom_command(OUTPUT TMust2DataDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts
add_custom_command(OUTPUT TMust2PhysicsReaderDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TMust2PhysicsReader.h TMust2PhysicsReaderDict.cxx TMust2PhysicsReader.rootmap libNPMUST2.dylib DEPENDS TMust2PhysicsReader.h)
add_library(NPMUST2 SHARED TMust2Data.cxx TMust2Physics.cxx TMust2PhysicsReader.cxx TMust2DataDict.cxx TMust2PhysicsDict.cxx TMust2PhysicsReaderDict.cxx TMust2Spectra.cxx)
#add_library(NPMUST2 SHARED TMust2Data.cxx TMust2Physics.cxx TMust2DataDict.cxx TMust2PhysicsDict.cxx TMust2Spectra.cxx)
-target_link_libraries(NPMUST2 ${ROOT_LIBRARIES} -lTreePlayer -lSpectrum NPCore)
+target_link_libraries(NPMUST2 ${ROOT_LIBRARIES} -lTreePlayer -lSpectrum NPCore NPPhysics)
#install(FILES TMust2Data.h TMust2Physics.h TMust2Spectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
install(FILES TMust2Data.h TMust2Physics.h TMust2Spectra.h TMust2PhysicsReader.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Detectors/MUST2/TMust2Physics.cxx b/NPLib/Detectors/MUST2/TMust2Physics.cxx
index 51a9a3754b3755ff77d040dbd7069d7054c0b582..30253a4f32cd5b2975fb74f70a4cb4aa87f3380f 100644
--- a/NPLib/Detectors/MUST2/TMust2Physics.cxx
+++ b/NPLib/Detectors/MUST2/TMust2Physics.cxx
@@ -417,6 +417,7 @@ void TMust2Physics::BuildPhysicalEvent() {
CsI_N.push_back(m_PreTreatedData->GetMMCsIECristalNbr(j));
CsI_E.push_back(m_PreTreatedData->GetMMCsIEEnergy(j));
CsI_T.push_back(-1000);
+ // std::cout << "BuildPhysical CSI " << CsI_N.at(CsI_N.size()-1) << " " << CsI_E.at(CsI_E.size()-1) << " " << CsI_E.size() << " " << "\n";
// Look for associate Time
for (unsigned int k = 0; k < m_CsITMult; ++k) {
// Same Cristal, Same Detector
@@ -469,6 +470,8 @@ void TMust2Physics::BuildPhysicalEvent() {
}
} // loop on event multiplicity
EventMultiplicity = TelescopeNumber.size();
+ //std::cout << "BuildPhysical SI " << Si_X.at(Si_X.size()-1) << " " << Si_E.at(Si_E.size()-1) << " " << "\n";
+ // std::cout << couple_size << std::endl;
return;
}
@@ -528,9 +531,13 @@ void TMust2Physics::PreTreat() {
for (unsigned int i = 0; i < m_CsIEMult; ++i) {
if (m_EventData->GetMMCsIEEnergy(i) > m_CsI_E_RAW_Threshold &&
IsValidChannel(3, m_EventData->GetMMCsIEDetectorNbr(i), m_EventData->GetMMCsIECristalNbr(i))) {
- double ECsI = fCsI_E(m_EventData, i);
- if (ECsI > m_CsI_E_Threshold)
+ // double ECsI = fCsI_E(m_EventData, i);
+ double ECsI = m_EventData->GetMMCsIEEnergy(i);
+ if (ECsI > 8192)
+ {
m_PreTreatedData->SetCsIE(m_EventData->GetMMCsIEDetectorNbr(i), m_EventData->GetMMCsIECristalNbr(i), ECsI);
+ // std::cout << "Pretreat " << m_PreTreatedData->GetMMCsIEDetectorNbr(i) << " " << m_PreTreatedData->GetMMCsIECristalNbr(i) << " " << ECsI << "\n";
+ }
}
}
@@ -1563,9 +1570,20 @@ void TMust2Physics::InitializeRootHistogramsCSIF(Int_t DetectorNumber){
TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
GATCONFMASTER_ = new TTreeReaderValue<unsigned short>(*inputTreeReader,"GATCONFMASTER");
}
+ string EnergyLossPath = NPOptionManager::getInstance()->GetEnergyLossPath();
+ string CutsPath = NPOptionManager::getInstance()->GetCutsPath();
+ DoCSIFit = true;
+
+ if(DoCSIFit){
+ for(unsigned int i = 0; i < ParticleType.size(); i++){
+ ParticleSi[ParticleType[i].c_str()] = new NPL::EnergyLoss(EnergyLossPath+ ParticleType[i].c_str()+"_Si.G4table", "G4Table", 100);
+ }
+ }
unsigned int NbCSI = 16;
auto TH2Map = RootHistogramsCalib::getInstance()->GetTH2Map();
auto TGraphMap = RootHistogramsCalib::getInstance()->GetTGraphMap();
+ auto TCutGMap = RootHistogramsCalib::getInstance()->GetTCutGMap();
+ auto TFileMap = RootHistogramsCalib::getInstance()->GetTFileMap();
for (Int_t j = 0; j < NbCSI; j++) {
TString hnameCSIE = Form("hMM%d_CSI_E%d", DetectorNumber, j+1);
@@ -1577,6 +1595,31 @@ void TMust2Physics::InitializeRootHistogramsCSIF(Int_t DetectorNumber){
(*TGraphMap)["MUST2"][hnameFITCSIE] = new TGraphErrors(3);
(*TGraphMap)["MUST2"][hnameFITCSIE]->SetTitle(htitleFITCSIE);
(*TGraphMap)["MUST2"][hnameFITCSIE]->SetName(hnameFITCSIE);
+
+ if(DoCSIFit){
+ string EnergyLossPath = NPOptionManager::getInstance()->GetEnergyLossPath();
+ string CutsPath = NPOptionManager::getInstance()->GetCutsPath();
+
+ for(unsigned int i = 0; i < ParticleType.size(); i++){
+ std::cout << ParticleType[i] << "\n";
+ TString CutName = Form("%s_hMM%u_CSI%u",ParticleType[i].c_str(), DetectorNumber, j+1);
+ TString cFileName = CutName+".root";
+ std::cout << CutName << " " << cFileName << " " << CutsPath+cFileName << "\n";
+
+ htitleCSIE = Form("%s_MM%u_CSI%u",ParticleType[i].c_str(), DetectorNumber, j+1);
+ (*TH2Map)["MUST2"][CutName] = new TH2F(CutName, htitleCSIE, 8192, 8192, 16384, 2000, 0, 200);
+
+ if((*TFileMap)["MUST2"][CutName] = new TFile(CutsPath+cFileName))
+ {
+ (*TCutGMap)["MUST2"][CutName]= (TCutG*)(*TFileMap)["MUST2"][CutName]->FindObjectAny(CutName);
+ }
+
+
+ }
+
+
+
+ }
}
}
@@ -1648,10 +1691,15 @@ void TMust2Physics::FillHistogramsCalib(){
if(NPOptionManager::getInstance()->IsReader())
m_EventData = &(**r_ReaderEventData);
FillHistogramsCalibEnergyF();
- std::cout << "Test pas bien" << std::endl;
}
+ //if(**GATCONFMASTER_ > 0)
+ //{
+ //test0++;
+ //std::cout << "test0 " << test0 << std::endl;
+ //}
if(IsCalibCSI){
+ // if(!IsCalibEnergy && NPOptionManager::getInstance()->IsReader())
if(!IsCalibEnergy && NPOptionManager::getInstance()->IsReader() && **(GATCONFMASTER_) > 0)
{
m_EventData = &(**r_ReaderEventData);
@@ -1705,6 +1753,7 @@ void TMust2Physics::FillHistogramsCalibEnergyF(){
void TMust2Physics::FillHistogramsCalibCSIF(){
DoCalibrationCSIPreTreat();
auto TH2Map = RootHistogramsCalib::getInstance()->GetTH2Map();
+ auto TCutGMap = RootHistogramsCalib::getInstance()->GetTCutGMap();
double matchSigma = m_StripEnergyMatchingSigma;
double NmatchSigma = m_StripEnergyMatchingNumberOfSigma;
@@ -1712,6 +1761,8 @@ void TMust2Physics::FillHistogramsCalibCSIF(){
unsigned int StripXMult =m_PreTreatedData->GetMMStripXEMult();
unsigned int StripYMult =m_PreTreatedData->GetMMStripYEMult();
+ // test1++;
+ // std::cout << "test1 " << test1 << std::endl;
// for(unsigned int ix = 0; ix < StripXMult; ix++){
// for(unsigned int iy = 0; iy < StripYMult; iy++){
if(StripXMult == 1 && StripYMult == 1){
@@ -1739,6 +1790,8 @@ void TMust2Physics::FillHistogramsCalibCSIF(){
unsigned int CSIDetNbr =m_PreTreatedData->GetMMCsIEDetectorNbr(icsi);
if(StripXDetNbr == CSIDetNbr){
+ // test2++;
+ // std::cout << "test2 " << test2 << std::endl;
unsigned int Cristal =m_PreTreatedData->GetMMCsIECristalNbr(icsi);
@@ -1746,6 +1799,26 @@ void TMust2Physics::FillHistogramsCalibCSIF(){
unsigned int CSIE =m_PreTreatedData->GetMMCsIEEnergy(icsi);
TString hnameCSIE = Form("hMM%d_CSI_E%d", CSIDetNbr, Cristal);
(*TH2Map)["MUST2"][hnameCSIE]->Fill(CSIE,StripXEnergy);
+ if(DoCSIFit){
+ Si_X.clear();
+ Si_Y.clear();
+ TelescopeNumber.clear();
+ Si_X.push_back(StripXNbr);
+ Si_Y.push_back(StripYNbr);
+ TelescopeNumber.push_back(DetNbr);
+ TVector3 BeamImpact = TVector3(0,0,0);
+ TVector3 HitDirection = GetPositionOfInteraction(0) - BeamImpact ;
+ double ThetaM2Surface = HitDirection.Angle(-GetTelescopeNormal(0) );
+ for(unsigned int i = 0; i < ParticleType.size(); i++){
+
+ TString CutName = Form("%s_hMM%u_CSI%u",ParticleType[i].c_str(), CSIDetNbr, Cristal);
+ if((*TCutGMap)["MUST2"][CutName] != 0 && (*TCutGMap)["MUST2"][CutName]->IsInside(CSIE,StripXEnergy)){
+ // test3++;
+ // std::cout << "test3 " << test3 << std::endl;
+ (*TH2Map)["MUST2"][CutName]->Fill(CSIE,(ParticleSi[ParticleType[i].c_str()]->EvaluateEnergyFromDeltaE(StripXEnergy, 300*um, ThetaM2Surface, 6.0 * MeV, 300.0 * MeV,0.001 * MeV, 10000)) - StripXEnergy);
+ }
+ }
+ }
}
}
}
@@ -1828,33 +1901,39 @@ void TMust2Physics::WriteHistogramsCalib(){
void TMust2Physics::WriteHistogramsCSIF(){
auto File = RootHistogramsCalib::getInstance()->GetFile();
auto TH2Map = RootHistogramsCalib::getInstance()->GetTH2Map();
+ auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();
auto TGraphMap = RootHistogramsCalib::getInstance()->GetTGraphMap();
unsigned int NbCSI = 16;
if(!File->GetDirectory("MUST2"))
File->mkdir("MUST2");
File->cd("MUST2");
- std::cout << "////////////////// TEST1\n";
map<int, bool>::iterator it;
for (it = DoCalibrationCsI.begin(); it != DoCalibrationCsI.end(); it++)
{
- std::cout << "////////////////// TEST2\n";
- std::cout << it->second <<"\n";
+ std::cout << "Writing Calibs for MUST2 : " << it->first << "\n";
if(it->second)
{
if(!gDirectory->GetDirectory(Form("M2_Telescope%d",it->first)))
{
- std::cout << "////////////////// TEST3\n";
gDirectory->mkdir(Form("M2_Telescope%d",it->first));
}
- std::cout << "////////////////// TEST4\n";
gDirectory->cd(Form("M2_Telescope%d",it->first));
gDirectory->mkdir("CSI");
//gDirectory->mkdir("Time");
gDirectory->cd("CSI");
for (Int_t j = 1; j <= NbCSI; j++) {
TString hnameCSIE = Form("hMM%d_CSI_E%d", it->first, j);
- (*TH2Map)["MUST2"][hnameCSIE]->Write();
+ (*TH2Map)["MUST2"][hnameCSIE]->Write();
+ if(DoCSIFit){
+ for(unsigned int i = 0; i < ParticleType.size(); i++){
+ TString CutName = Form("%s_hMM%u_CSI%u",ParticleType[i].c_str(), it->first, j);
+ (*TH2Map)["MUST2"][CutName]->Write();
+ (*TH1Map)["MUST2"][CutName+"_0"]->Write();
+ (*TH1Map)["MUST2"][CutName+"_1"]->Write();
+ (*TH1Map)["MUST2"][CutName+"_2"]->Write();
+ }
+ }
//(*TGraphMap)["MUST2"][hnameFITXE]->Write();
//(*TGraphMap)["MUST2"][hnameFITYE]->Write();
}
@@ -1917,7 +1996,7 @@ void TMust2Physics::DoCalibration(){
{
if(it->second)
{
- MakeCalibFolders();
+ MakeEnergyCalibFolders();
DoCalibrationEnergyF(it->first);
}
}
@@ -1931,7 +2010,8 @@ void TMust2Physics::DoCalibration(){
for (it = DoCalibrationCsI.begin(); it != DoCalibrationCsI.end(); it++)
{
if(it->second)
- {
+ {
+ MakeCSICalibFolders();
DoCalibrationCsIF(it->first);
}
}
@@ -2224,7 +2304,7 @@ Double_t TMust2Physics::source_Cm(Double_t *x, Double_t *par)
return fitval;
}
-void TMust2Physics::MakeCalibFolders(){
+void TMust2Physics::MakeEnergyCalibFolders(){
std::string Path = NPOptionManager::getInstance()->GetCalibrationOutputPath();
std::string OutputName = NPOptionManager::getInstance()->GetOutputFile();
if(OutputName.size() > 5){
@@ -2241,6 +2321,21 @@ void TMust2Physics::MakeCalibFolders(){
system(make_folder+"/latex");
system(make_folder+"/latex/pictures");
}
+
+void TMust2Physics::MakeCSICalibFolders(){
+ std::string Path = NPOptionManager::getInstance()->GetCalibrationOutputPath();
+ std::string OutputName = NPOptionManager::getInstance()->GetOutputFile();
+ if(OutputName.size() > 5){
+ if(OutputName.substr(OutputName.size()-5,OutputName.size()) == ".root"){
+ OutputName = OutputName.substr(0,OutputName.size()-5);
+ }
+ }
+ TString test_folder = "test -f "+Path+OutputName;
+ TString make_folder = "mkdir "+Path+OutputName;
+
+ system(make_folder);
+}
+
void TMust2Physics::CreateCalibrationEnergyFiles(unsigned int DetectorNumber, TString side, ofstream *calib_file, ofstream *dispersion_file){
std::string Path = NPOptionManager::getInstance()->GetCalibrationOutputPath();
std::string OutputName = NPOptionManager::getInstance()->GetOutputFile();
@@ -2256,7 +2351,7 @@ void TMust2Physics::CreateCalibrationEnergyFiles(unsigned int DetectorNumber, TS
(*dispersion_file).open( ( (string)(Path+OutputName+"/dispersion/"+Filename+".dispersion") ).c_str() );
}
-void TMust2Physics::CreateCalibrationCSIFiles(unsigned int DetectorNumber, ofstream *calib_file){
+void TMust2Physics::CreateCalibrationCSIFiles(unsigned int DetectorNumber, ofstream *calib_file, TString ParticleType){
std::string Path = NPOptionManager::getInstance()->GetCalibrationOutputPath();
std::string OutputName = NPOptionManager::getInstance()->GetOutputFile();
if(OutputName.size() > 5){
@@ -2264,8 +2359,14 @@ void TMust2Physics::CreateCalibrationCSIFiles(unsigned int DetectorNumber, ofstr
OutputName = OutputName.substr(0,OutputName.size()-5);
}
}
- TString Filename = "CsI_MM"+std::to_string(DetectorNumber);
+ TString Filename = "CsI_MM"+std::to_string(DetectorNumber)+"_"+ParticleType;
+ std::cout << Filename << std::endl;
(*calib_file).open( ( (string)(Path+OutputName+"/"+Filename+".txt") ).c_str() );
+ std::cout << (string)(Path+OutputName+"/"+Filename+".txt") << std::endl;
+}
+
+void TMust2Physics::CloseCalibrationCSIFiles(ofstream *calib_file){
+ calib_file->close();
}
void TMust2Physics::CloseCalibrationEnergyFiles(ofstream *calib_file, ofstream *dispersion_file){
@@ -2281,12 +2382,49 @@ void TMust2Physics::DoCalibrationTimeF(Int_t DetectorNumber){
}
void TMust2Physics::DoCalibrationCsIF(Int_t DetectorNumber){
+ TF1 *Gaus = new TF1("Gaus","gaus",0,200);
+ TF1 *f1 = new TF1("f1","pol1",8192,16384);
+ f1->SetParLimits(2,0,1e5);
+ f1->SetParameter(2,0.01);
+ auto File = new TFile("./FitSlices.root","RECREATE");
+ auto TH2Map = RootHistogramsCalib::getInstance()->GetTH2Map();
+ auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();
+ auto TGraphMap = RootHistogramsCalib::getInstance()->GetTGraphMap();
+ unsigned int NbCSI = 16;
+
+ ofstream *calib_file = new ofstream;
+
+ for(unsigned int j = 0; j < ParticleType.size(); j++){
+ CreateCalibrationCSIFiles(DetectorNumber, calib_file, ParticleType[j]);
+ for(unsigned int i = 1; i <= NbCSI; i++){
+ TString CutName = Form("%s_hMM%u_CSI%u",ParticleType[j].c_str(), DetectorNumber, i);
+ TString htitleCSIE = Form("%s_hMM%u_CSI%u",ParticleType[j].c_str(), DetectorNumber, i);
+ double a = 0;
+ double b = 0;
+ if((*TH2Map)["MUST2"][CutName] != 0){
+ std::cout << "test1" << std::endl;
+ (*TH2Map)["MUST2"][CutName]->FitSlicesY(Gaus,0,2000,0,"QNG5",0);
+ std::cout << "test2" << std::endl;
+ (*TH1Map)["MUST2"][CutName+"_0"]= (TH1F*)File->Get(htitleCSIE+"_0");
+ File->Write();
+ (*TH1Map)["MUST2"][CutName+"_0"]->Draw();
+ (*TH1Map)["MUST2"][CutName+"_1"]= (TH1F*)File->Get(htitleCSIE+"_1");
+ (*TH1Map)["MUST2"][CutName+"_2"]= (TH1F*)File->Get(htitleCSIE+"_2");
+ std::cout << "test3" << std::endl;
+ (*TH1Map)["MUST2"][CutName+"_1"]->Fit(f1,"","",8192,16384);
+ a = f1->GetParameter(0);
+ b = f1->GetParameter(1);
+ }
+ *calib_file << "MUST2_T" << DetectorNumber << "_CsI" << i << "_E " << a << " " << b << endl ;
+
+ }
+ CloseCalibrationCSIFiles(calib_file);
+ }
}
/*
auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();
auto TGraphMap = RootHistogramsCalib::getInstance()->GetTGraphMap();
unsigned int NbCsI = 16;
- ofstream *calib_file = new ofstream;
// ExtractCutsAndFillTree();
diff --git a/NPLib/Detectors/MUST2/TMust2Physics.h b/NPLib/Detectors/MUST2/TMust2Physics.h
index 57575a21163efed5e7c6884d551c5c16989bc99f..26a5ce69b63de09720c644fb2ae9ac6e4c2e6420 100644
--- a/NPLib/Detectors/MUST2/TMust2Physics.h
+++ b/NPLib/Detectors/MUST2/TMust2Physics.h
@@ -28,6 +28,7 @@
#include <stdlib.h>
// NPL
#include "NPCalibrationManager.h"
+#include "NPEnergyLoss.h"
#include "NPInputParser.h"
#include "NPVDetector.h"
#include "NPVTreeReader.h"
@@ -62,7 +63,7 @@ public:
bool Match_Si_CsI(int X, int Y, int CristalNbr, int DetectorNbr);
bool Match_Si_SiLi(int X, int Y, int PadNbr);
bool ResolvePseudoEvent();
- bool m_reader = true;
+ bool m_reader = true;//!
public:
// Provide Physical Multiplicity
@@ -79,7 +80,9 @@ public:
vector<double> Si_E;
vector<double> Si_T;
vector<int> Si_X;
+ vector<int> Si_StripNumberX;
vector<int> Si_Y;
+ vector<int> Si_StripNumberY;
// Use for checking purpose
vector<double> Si_EX;
@@ -154,13 +157,15 @@ public: // Innherited from VDetector Class
void DoCalibrationCsIF(Int_t DetectorNumber);//!
- void MakeCalibFolders();//!
+ void MakeEnergyCalibFolders();//!
+
+ void MakeCSICalibFolders();//!
void CreateCalibrationEnergyFiles(unsigned int DetectorNumber, TString side, ofstream *calib_file, ofstream *dispersion_file);//!
- void CreateCalibrationCSIFiles(unsigned int DetectorNumber, ofstream *calib_file);//!
+ void CreateCalibrationCSIFiles(unsigned int DetectorNumber, ofstream *calib_file, TString ParticleType);//!
- void CloseCalibrationCSIFiles(ofstream *calib_file){};//!
+ void CloseCalibrationCSIFiles(ofstream *calib_file);//!
void CloseCalibrationEnergyFiles(ofstream *calib_file, ofstream *dispersion_file);//!
@@ -400,9 +405,14 @@ private:
map<int,double> CSIEnergyYThreshold;//!
map<int,double> CSIEThreshold;//!
TTreeReaderValue<unsigned short>* GATCONFMASTER_;//!
+ bool DoCSIFit;//!
+ std::map<TString,NPL::EnergyLoss*> ParticleSi;//!
+ // std::vector<string> ParticleType{"proton","deuteron","triton","3He","alpha"};
+ std::vector<string> ParticleType{"proton","deuteron","triton","alpha"};//!
// map<int,std::string> CalibFile;//!
+
private: // Spectra Class
TMust2Spectra* m_Spectra; //!
diff --git a/NPLib/Detectors/SEASON/TSEASONPhysics.cxx b/NPLib/Detectors/SEASON/TSEASONPhysics.cxx
index 7c432dc1e42dc5bd834f3662ad805545ae94dd68..c0c6d19e7e06505d9cbae2e23107e1830f384b7e 100644
--- a/NPLib/Detectors/SEASON/TSEASONPhysics.cxx
+++ b/NPLib/Detectors/SEASON/TSEASONPhysics.cxx
@@ -17,494 +17,463 @@
* *
*---------------------------------------------------------------------------*
* Comment: *
- * *
+ * *
* *
*****************************************************************************/
-
- #include "TSEASONPhysics.h"
-
- // STL
- #include <sstream>
- #include <iostream>
- #include <cmath>
- #include <stdlib.h>
- #include <limits>
- using namespace std;
-
- // NPL
- #include "RootInput.h"
- #include "RootOutput.h"
- #include "NPDetectorFactory.h"
- #include "NPOptionManager.h"
-
- // ROOT
- #include "TChain.h"
-
- ClassImp(TSEASONPhysics)
-
-
- ///////////////////////////////////////////////////////////////////////////
- TSEASONPhysics::TSEASONPhysics()
- : m_EventData(new TSEASONData),
- m_PreTreatedData(new TSEASONData),
- m_EventPhysics(this),
- m_Spectra(0),
- m_E_RAW_Threshold(0), // adc channels
- m_E_Threshold(0), // MeV
- m_NumberOfDetectors(0) {
- }
-
- ///////////////////////////////////////////////////////////////////////////
- /// A usefull method to bundle all operation to add a detector
- void TSEASONPhysics::AddDetector(TVector3 X1_Y1,
- TVector3 X1_YMax,
- TVector3 XMax_Y1,
- TVector3 XMax_YMax,
- int NumberOfStripsX,
- int NumberOfStripsY){
- // In That simple case nothing is done
- // Typically for more complex detector one would calculate the relevant
- // positions (stripped silicon) or angles (gamma array)
- m_NumberOfDetectors++;
- m_NumberOfStripsX = NumberOfStripsX;
- m_NumberOfStripsY = NumberOfStripsY;
- // Vector U on Telescope Face (parallele to 0X axis) (parallele to X Strip) (NB: remember that Y strip are allong X axis)
- TVector3 U = XMax_Y1 - X1_Y1 ;
- double FaceX = U.Mag();
- double Ushift = (U.Mag()-FaceX)/2.;
- U=U.Unit();
-
- // Vector V on Telescope Face (parallele to 0Y axis) (parallele to Y Strip)
- TVector3 V = X1_YMax - X1_Y1 ;
- double FaceY = V.Mag();
- double Vshift = (V.Mag()-FaceY)/2.;
- V=V.Unit();
-
- // Position Vector of Strip Center
- TVector3 StripCenter = TVector3(0,0,0);
- // Position Vector of X=1 Y=1 Strip
- TVector3 Strip_1_1;
-
- // Geometry Parameter
- double StripPitchX = FaceX/NumberOfStripsX ; //mm
- double StripPitchY = FaceY/NumberOfStripsY ; //mm
-
- // buffer object to fill Position Array
- vector<double> lineX ; vector<double> lineY ; vector<double> lineZ ;
-
- vector< vector < double > > OneTelescopeStripPositionX;
- vector< vector < double > > OneTelescopeStripPositionY;
- vector< vector < double > > OneTelescopeStripPositionZ;
-
- // Moving StripCEnter to 1.1 corner
-
- Strip_1_1 = X1_Y1 + U*(StripPitchX/2.) + V*(StripPitchY/2.);
- Strip_1_1 += U*Ushift+V*Vshift;
-
- for(int i = 0 ; i < NumberOfStripsX ; ++i){
- lineX.clear();
- lineY.clear();
- lineZ.clear();
-
- for( int j = 0 ; j < NumberOfStripsY ; ++j){
- StripCenter = Strip_1_1 + StripPitchX*i*U + StripPitchY*i*V;
- lineX.push_back( StripCenter.X());
- lineY.push_back( StripCenter.Y());
- lineZ.push_back( StripCenter.Z());
- }
-
- OneTelescopeStripPositionX.push_back(lineX);
- OneTelescopeStripPositionY.push_back(lineY);
- OneTelescopeStripPositionZ.push_back(lineZ);
-
- }
-
- m_StripPositionX.push_back(OneTelescopeStripPositionX);
- m_StripPositionY.push_back(OneTelescopeStripPositionY);
- m_StripPositionZ.push_back(OneTelescopeStripPositionZ);
-
- }
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::BuildSimplePhysicalEvent() {
- BuildPhysicalEvent();
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::BuildPhysicalEvent() {
- // apply thresholds and calibration
- PreTreat();
-
- // match energy and time together
- unsigned int mysizeXE = m_PreTreatedData->GetXMultEnergy();
- unsigned int mysizeYE = m_PreTreatedData->GetYMultEnergy();
- unsigned int mysizeXT = m_PreTreatedData->GetXMultTime();
- unsigned int mysizeYT = m_PreTreatedData->GetYMultTime();
- unsigned int mysizeE = m_PreTreatedData->GetMultEnergy();
- unsigned int mysizeT = m_PreTreatedData->GetMultTime();
-
- for (UShort_t e = 0; e < mysizeXE ; e++) {
- for (UShort_t t = 0; t < mysizeXT ; t++) {
- if (m_PreTreatedData->GetXE_DetectorNbr(e) == m_PreTreatedData->GetXT_DetectorNbr(t)
- && m_PreTreatedData->GetXE_StripNbr(e) == m_PreTreatedData->GetXT_StripNbr(t)) {
- DetectorNumberX.push_back(m_PreTreatedData->GetXE_DetectorNbr(e));
- StripNumberX.push_back(m_PreTreatedData->GetXE_StripNbr(e));
- EnergyX.push_back(m_PreTreatedData->GetX_Energy(e));
- TimeX.push_back(m_PreTreatedData->GetX_Time(t));
- }
- }
- }
- for (UShort_t e = 0; e < mysizeYE ; e++) {
- for (UShort_t t = 0; t < mysizeYT ; t++) {
- if (m_PreTreatedData->GetYE_DetectorNbr(e) == m_PreTreatedData->GetYT_DetectorNbr(t)
- && m_PreTreatedData->GetYE_StripNbr(e) == m_PreTreatedData->GetYT_StripNbr(t)) {
- DetectorNumberY.push_back(m_PreTreatedData->GetYE_DetectorNbr(e));
- StripNumberY.push_back(m_PreTreatedData->GetYE_StripNbr(e));
- EnergyY.push_back(m_PreTreatedData->GetY_Energy(e));
- TimeY.push_back(m_PreTreatedData->GetY_Time(t));
- }
- }
- }
-
- 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)
- && m_PreTreatedData->GetE_StripNbrX(e) == m_PreTreatedData->GetT_StripNbrX(t)
- && m_PreTreatedData->GetE_StripNbrY(e) == m_PreTreatedData->GetT_StripNbrY(t)) {
- DetectorNumber.push_back(m_PreTreatedData->GetE_DetectorNbr(e));
- StripX.push_back(m_PreTreatedData->GetE_StripNbrX(e));
- StripY.push_back(m_PreTreatedData->GetE_StripNbrY(e));
- Energy.push_back(m_PreTreatedData->Get_Energy(e));
- Time.push_back(m_PreTreatedData->Get_Time(t));
- }
- }
- }
- //DetectorNumber = DetectorNumberX;
- }
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::PreTreat() {
- // This method typically applies thresholds and calibrations
- // Might test for disabled channels for more complex detector
-
- // clear pre-treated object
- ClearPreTreatedData();
- unsigned int mysizeX = m_EventData->GetXMultEnergy();
- unsigned int mysizeY = m_EventData->GetYMultEnergy();
-
- // instantiate CalibrationManager
- static CalibrationManager* Cal = CalibrationManager::getInstance();
-
- if(mysizeX == mysizeY){
- unsigned int mysize = mysizeX;
- // Energy
- for (UShort_t i = 0; i < mysize ; ++i) {
- if(m_EventData->GetXE_DetectorNbr(i) == m_EventData->GetYE_DetectorNbr(i)){
- if (m_EventData->GetX_Energy(i) > m_E_RAW_Threshold) {
- double x = (m_EventData->GetXE_StripNbr(i) - 16.5)*2, y = (m_EventData->GetYE_StripNbr(i) - 16.5)*2;
- //int dist = std::round(TMath::Sqrt(x*x+y*y+4)-2);
- int dist = std::round(TMath::Sqrt(y*y+4)-2);
- //if(dist > 38) dist = 38;
- string name = "SEASON/D";
- name += NPL::itoa(m_EventData->GetXE_DetectorNbr(i));
- name += "_DIST";
- name += NPL::itoa(dist);
- name += "_ENERGY";
- string namestrip = "SEASON/D";
- namestrip += NPL::itoa(m_EventData->GetXE_DetectorNbr(i));
- namestrip += "_STRIPX";
- namestrip += NPL::itoa(m_EventData->GetXE_StripNbr(i));
- namestrip += "_ENERGY";
- //cout << name << endl;
- Double_t Estrip = Cal->ApplyCalibration(namestrip,m_EventData->GetX_Energy(i));
- Double_t E = Cal->ApplyCalibration(name,Estrip);
- if (E > m_XE_Threshold) {
- m_PreTreatedData->SetEnergy(m_EventData->GetXE_DetectorNbr(i),m_EventData->GetXE_StripNbr(i),m_EventData->GetYE_StripNbr(i),E);
- }
- }
- }
- }
-
- // Time
- for (UShort_t i = 0; i < mysize; ++i) {
- if(m_EventData->GetXT_DetectorNbr(i) == m_EventData->GetYT_DetectorNbr(i)){
- Double_t T= Cal->ApplyCalibration("SEASON/D"+NPL::itoa(m_EventData->GetXT_DetectorNbr(i))+"_STRIPX"+NPL::itoa(m_EventData->GetXT_StripNbr(i))+"_STRIPY"+NPL::itoa(m_EventData->GetYT_StripNbr(i))+"_Time",m_EventData->GetX_Time(i));
- m_PreTreatedData->SetTime(m_EventData->GetXT_DetectorNbr(i),m_EventData->GetXT_StripNbr(i),m_EventData->GetYT_StripNbr(i),T);
- }
- }
- }
-
- // X Energy
- for (UShort_t i = 0; i < mysizeX ; ++i) {
- if (m_EventData->GetX_Energy(i) > m_E_RAW_Threshold) {
- Double_t EX = Cal->ApplyCalibration("SEASON/D"+NPL::itoa(m_EventData->GetXE_DetectorNbr(i))+"_STRIP"+NPL::itoa(m_EventData->GetXE_StripNbr(i))+"_Energy",m_EventData->GetX_Energy(i));
- if (EX > m_XE_Threshold) {
- m_PreTreatedData->SetXEnergy(m_EventData->GetXE_DetectorNbr(i),m_EventData->GetXE_StripNbr(i) ,EX);
- }
- }
- }
-
- // X Time
- for (UShort_t i = 0; i < mysizeX; ++i) {
- Double_t TX= Cal->ApplyCalibration("SEASON/D"+NPL::itoa(m_EventData->GetXT_DetectorNbr(i))+"_Strip"+NPL::itoa(m_EventData->GetXT_StripNbr(i))+"_Time",m_EventData->GetX_Time(i));
- m_PreTreatedData->SetXTime(m_EventData->GetXT_DetectorNbr(i),m_EventData->GetXT_StripNbr(i) ,TX);
- }
-
- // Y Energy
- for (UShort_t i = 0; i < mysizeY; ++i) {
- if (m_EventData->GetY_Energy(i) > m_E_RAW_Threshold) {
- Double_t EY = Cal->ApplyCalibration("SEASON/D"+NPL::itoa(m_EventData->GetYE_DetectorNbr(i))+"_STRIP"+NPL::itoa(m_EventData->GetYE_StripNbr(i))+"_Energy",m_EventData->GetY_Energy(i));
- if (EY > m_YE_Threshold) {
- m_PreTreatedData->SetYEnergy(m_EventData->GetYE_DetectorNbr(i),m_EventData->GetYE_StripNbr(i) ,EY);
- }
- }
- }
-
- // Y Time
- for (UShort_t i = 0; i < mysizeY; ++i) {
- Double_t TY= Cal->ApplyCalibration("SEASON/D"+NPL::itoa(m_EventData->GetYT_DetectorNbr(i))+"_Strip"+NPL::itoa(m_EventData->GetYT_StripNbr(i))+"_Time",m_EventData->GetY_Time(i));
- m_PreTreatedData->SetYTime(m_EventData->GetYT_DetectorNbr(i),m_EventData->GetYT_StripNbr(i) ,TY);
- }
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::ReadAnalysisConfig() {
- bool ReadingStatus = false;
-
- // path to file
- string FileName = "./configs/ConfigSEASON.dat";
-
- // open analysis config file
- ifstream AnalysisConfigFile;
- AnalysisConfigFile.open(FileName.c_str());
-
- if (!AnalysisConfigFile.is_open()) {
- cout << " No ConfigSEASON.dat found: Default parameter loaded for Analayis " << FileName << endl;
- return;
- }
- cout << " Loading user parameter for Analysis from ConfigSEASON.dat " << endl;
-
- // Save it in a TAsciiFile
- TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
- asciiConfig->AppendLine("%%% ConfigSEASON.dat %%%");
- asciiConfig->Append(FileName.c_str());
- asciiConfig->AppendLine("");
- // read analysis config file
- string LineBuffer,DataBuffer,whatToDo;
- while (!AnalysisConfigFile.eof()) {
- // Pick-up next line
- getline(AnalysisConfigFile, LineBuffer);
-
- // search for "header"
- string name = "ConfigSEASON";
- if (LineBuffer.compare(0, name.length(), name) == 0)
- ReadingStatus = true;
-
- // loop on tokens and data
- while (ReadingStatus ) {
- whatToDo="";
- AnalysisConfigFile >> whatToDo;
-
- // Search for comment symbol (%)
- if (whatToDo.compare(0, 1, "%") == 0) {
- AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
- }
- else {
- ReadingStatus = false;
- }
- }
- }
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::Clear() {
- DetectorNumber.clear();
- StripNumberX.clear();
- StripNumberY.clear();
- StripX.clear();
- StripY.clear();
- Energy.clear();
- Time.clear();
-
- DetectorNumberX.clear();
- DetectorNumberY.clear();
- EnergyX.clear();
- EnergyY.clear();
- TimeX.clear();
- TimeY.clear();
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::ReadConfiguration(NPL::InputParser parser) {
- vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SEASON");
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
-
- vector<string> cart = {"X1_Y1","X1_YMax","XMax_Y1","XMax_YMax","NStripsX","NStripsY","Group"};
-
- for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(cart)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// SEASON " << i+1 << endl;
-
- TVector3 A = blocks[i]->GetTVector3("X1_Y1","mm");
- TVector3 B = blocks[i]->GetTVector3("X1_YMax","mm");
- TVector3 C = blocks[i]->GetTVector3("XMax_Y1","mm");
- TVector3 D = blocks[i]->GetTVector3("XMax_YMax","mm");
- int NumberOfStripsX = blocks[i]->GetInt("NStripsX");
- int NumberOfStripsY = blocks[i]->GetInt("NStripsY");
- AddDetector(A,B,C,D, NumberOfStripsX, NumberOfStripsY);
- }
- else{
- cout << "ERROR: check your input file formatting " << endl;
- exit(1);
- }
- }
- }
- ///////////////////////////////////////////////////////////////////////////
- TVector3 TSEASONPhysics::GetPositionOfInteraction(const int i) const
- {
- TVector3 Position;
- Position = TVector3( GetStripPositionX( DetectorNumber[i], StripNumberX[i], StripNumberY[i]),
- GetStripPositionY( DetectorNumber[i], StripNumberX[i], StripNumberY[i]),
- GetStripPositionZ( DetectorNumber[i], StripNumberX[i], StripNumberY[i]));
- return Position;
- }
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::InitSpectra() {
- m_Spectra = new TSEASONSpectra(m_NumberOfDetectors, m_NumberOfStripsX, m_NumberOfStripsY);
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::FillSpectra() {
- m_Spectra -> FillRawSpectra(m_EventData);
- m_Spectra -> FillPreTreatedSpectra(m_PreTreatedData);
- m_Spectra -> FillPhysicsSpectra(m_EventPhysics);
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::CheckSpectra() {
- m_Spectra->CheckSpectra();
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::ClearSpectra() {
- // To be done
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- map< string , TH1*> TSEASONPhysics::GetSpectra() {
- if(m_Spectra)
- return m_Spectra->GetMapHisto();
- else{
- map< string , TH1*> empty;
- return empty;
- }
- }
-
- ///////////////////////////////////////////////////////////////////////////
- vector<TCanvas*> TSEASONPhysics::GetCanvas() {
- /*if(m_Spectra)
- return m_Spectra->GetCanvas();
- else{
- vector<TCanvas*> empty;
- return empty;
- }*/
- }
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::WriteSpectra() {
- m_Spectra->WriteSpectra();
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::AddParameterToCalibrationManager() {
- int NumberOfStrips = 32;
- CalibrationManager* Cal = CalibrationManager::getInstance();
- vector<double> standard = {0,1};
- for (int i = 0; i < m_NumberOfDetectors; ++i) {
- for(int k = 0;k<50;k++){
- Cal->AddParameter("SEASON", "D"+NPL::itoa(i+1)+"_DIST"+NPL::itoa(k)+"_ENERGY","SEASON_D"+ NPL::itoa(i+1)+"_DIST"+NPL::itoa(k)+"_ENERGY", standard);
- Cal->AddParameter("SEASON", "D"+NPL::itoa(i+1)+"_FOILDIST"+NPL::itoa(k)+"_ENERGY","SEASON_D"+ NPL::itoa(i+1)+"_FOILDIST"+NPL::itoa(k)+"_ENERGY", standard);
- }
-
- for(int j = 0; j < NumberOfStrips; ++j){
- Cal->AddParameter("SEASON", "D"+NPL::itoa(i+1)+"_STRIPX"+NPL::itoa(j+1)+"_ENERGY",
- "SEASON_D"+ NPL::itoa(i+1)+"_STRIPX"+NPL::itoa(j+1)+"_ENERGY");
- Cal->AddParameter("SEASON", "D"+ NPL::itoa(i+1)+"_STRIPX"+NPL::itoa(j+1)+"_TIME","SEASON_D"+ NPL::itoa(i+1)+"_STRIP"+NPL::itoa(j+1)+"_TIME");
- }
- }
- }
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::InitializeRootInputRaw() {
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchStatus("SEASON", true );
- inputChain->SetBranchAddress("SEASON", &m_EventData );
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::InitializeRootInputPhysics() {
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchAddress("SEASON", &m_EventPhysics);
- }
-
-
-
- ///////////////////////////////////////////////////////////////////////////
- void TSEASONPhysics::InitializeRootOutput() {
- TTree* outputTree = RootOutput::getInstance()->GetTree();
- outputTree->Branch("SEASON", "TSEASONPhysics", &m_EventPhysics);
- }
-
-
-
- ////////////////////////////////////////////////////////////////////////////////
- // Construct Method to be pass to the DetectorFactory //
- ////////////////////////////////////////////////////////////////////////////////
- NPL::VDetector* TSEASONPhysics::Construct() {
- return (NPL::VDetector*) new TSEASONPhysics();
- }
-
-
-
- ////////////////////////////////////////////////////////////////////////////////
- // Registering the construct method to the factory //
- ////////////////////////////////////////////////////////////////////////////////
- extern "C"{
- class proxy_SEASON{
- public:
- proxy_SEASON(){
- NPL::DetectorFactory::getInstance()->AddToken("SEASON","SEASON");
- NPL::DetectorFactory::getInstance()->AddDetector("SEASON",TSEASONPhysics::Construct);
- }
- };
-
- proxy_SEASON p_SEASON;
- }
+
+#include "TSEASONPhysics.h"
+
+// STL
+#include <cmath>
+#include <iostream>
+#include <limits>
+#include <sstream>
+#include <stdlib.h>
+using namespace std;
+
+// NPL
+#include "NPDetectorFactory.h"
+#include "NPOptionManager.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+
+// ROOT
+#include "TChain.h"
+
+ClassImp(TSEASONPhysics)
+
+ ///////////////////////////////////////////////////////////////////////////
+ TSEASONPhysics::TSEASONPhysics()
+ : m_EventData(new TSEASONData), m_PreTreatedData(new TSEASONData), m_EventPhysics(this), m_Spectra(0),
+ m_E_RAW_Threshold(0), // adc channels
+ m_E_Threshold(0), // MeV
+ m_NumberOfDetectors(0) {}
+
+///////////////////////////////////////////////////////////////////////////
+/// A usefull method to bundle all operation to add a detector
+void TSEASONPhysics::AddDetector(TVector3 X1_Y1, TVector3 X1_YMax, TVector3 XMax_Y1, TVector3 XMax_YMax,
+ int NumberOfStripsX, int NumberOfStripsY) {
+ // In That simple case nothing is done
+ // Typically for more complex detector one would calculate the relevant
+ // positions (stripped silicon) or angles (gamma array)
+ m_NumberOfDetectors++;
+ m_NumberOfStripsX = NumberOfStripsX;
+ m_NumberOfStripsY = NumberOfStripsY;
+ // Vector U on Telescope Face (parallele to 0X axis) (parallele to X Strip) (NB: remember that Y strip are allong X
+ // axis)
+ TVector3 U = XMax_Y1 - X1_Y1;
+ double FaceX = U.Mag();
+ double Ushift = (U.Mag() - FaceX) / 2.;
+ U = U.Unit();
+
+ // Vector V on Telescope Face (parallele to 0Y axis) (parallele to Y Strip)
+ TVector3 V = X1_YMax - X1_Y1;
+ double FaceY = V.Mag();
+ double Vshift = (V.Mag() - FaceY) / 2.;
+ V = V.Unit();
+
+ // Position Vector of Strip Center
+ TVector3 StripCenter = TVector3(0, 0, 0);
+ // Position Vector of X=1 Y=1 Strip
+ TVector3 Strip_1_1;
+
+ // Geometry Parameter
+ double StripPitchX = FaceX / NumberOfStripsX; // mm
+ double StripPitchY = FaceY / NumberOfStripsY; // mm
+
+ // buffer object to fill Position Array
+ vector<double> lineX;
+ vector<double> lineY;
+ vector<double> lineZ;
+
+ vector<vector<double>> OneTelescopeStripPositionX;
+ vector<vector<double>> OneTelescopeStripPositionY;
+ vector<vector<double>> OneTelescopeStripPositionZ;
+
+ // Moving StripCEnter to 1.1 corner
+
+ Strip_1_1 = X1_Y1 + U * (StripPitchX / 2.) + V * (StripPitchY / 2.);
+ Strip_1_1 += U * Ushift + V * Vshift;
+
+ for (int i = 0; i < NumberOfStripsX; ++i) {
+ lineX.clear();
+ lineY.clear();
+ lineZ.clear();
+
+ for (int j = 0; j < NumberOfStripsY; ++j) {
+ StripCenter = Strip_1_1 + StripPitchX * i * U + StripPitchY * i * V;
+ lineX.push_back(StripCenter.X());
+ lineY.push_back(StripCenter.Y());
+ lineZ.push_back(StripCenter.Z());
+ }
+
+ OneTelescopeStripPositionX.push_back(lineX);
+ OneTelescopeStripPositionY.push_back(lineY);
+ OneTelescopeStripPositionZ.push_back(lineZ);
+ }
+
+ m_StripPositionX.push_back(OneTelescopeStripPositionX);
+ m_StripPositionY.push_back(OneTelescopeStripPositionY);
+ m_StripPositionZ.push_back(OneTelescopeStripPositionZ);
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::BuildSimplePhysicalEvent() { BuildPhysicalEvent(); }
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::BuildPhysicalEvent() {
+ // apply thresholds and calibration
+ PreTreat();
+
+ // match energy and time together
+ unsigned int mysizeXE = m_PreTreatedData->GetXMultEnergy();
+ unsigned int mysizeYE = m_PreTreatedData->GetYMultEnergy();
+ unsigned int mysizeXT = m_PreTreatedData->GetXMultTime();
+ unsigned int mysizeYT = m_PreTreatedData->GetYMultTime();
+ unsigned int mysizeE = m_PreTreatedData->GetMultEnergy();
+ unsigned int mysizeT = m_PreTreatedData->GetMultTime();
+
+ for (UShort_t e = 0; e < mysizeXE; e++) {
+ for (UShort_t t = 0; t < mysizeXT; t++) {
+ if (m_PreTreatedData->GetXE_DetectorNbr(e) == m_PreTreatedData->GetXT_DetectorNbr(t) &&
+ m_PreTreatedData->GetXE_StripNbr(e) == m_PreTreatedData->GetXT_StripNbr(t)) {
+ DetectorNumberX.push_back(m_PreTreatedData->GetXE_DetectorNbr(e));
+ StripNumberX.push_back(m_PreTreatedData->GetXE_StripNbr(e));
+ EnergyX.push_back(m_PreTreatedData->GetX_Energy(e));
+ TimeX.push_back(m_PreTreatedData->GetX_Time(t));
+ }
+ }
+ }
+ for (UShort_t e = 0; e < mysizeYE; e++) {
+ for (UShort_t t = 0; t < mysizeYT; t++) {
+ if (m_PreTreatedData->GetYE_DetectorNbr(e) == m_PreTreatedData->GetYT_DetectorNbr(t) &&
+ m_PreTreatedData->GetYE_StripNbr(e) == m_PreTreatedData->GetYT_StripNbr(t)) {
+ DetectorNumberY.push_back(m_PreTreatedData->GetYE_DetectorNbr(e));
+ StripNumberY.push_back(m_PreTreatedData->GetYE_StripNbr(e));
+ EnergyY.push_back(m_PreTreatedData->GetY_Energy(e));
+ TimeY.push_back(m_PreTreatedData->GetY_Time(t));
+ }
+ }
+ }
+
+ 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) &&
+ m_PreTreatedData->GetE_StripNbrX(e) == m_PreTreatedData->GetT_StripNbrX(t) &&
+ m_PreTreatedData->GetE_StripNbrY(e) == m_PreTreatedData->GetT_StripNbrY(t)) {
+ DetectorNumber.push_back(m_PreTreatedData->GetE_DetectorNbr(e));
+ StripX.push_back(m_PreTreatedData->GetE_StripNbrX(e));
+ StripY.push_back(m_PreTreatedData->GetE_StripNbrY(e));
+ Energy.push_back(m_PreTreatedData->Get_Energy(e));
+ Time.push_back(m_PreTreatedData->Get_Time(t));
+ }
+ }
+ }
+ // DetectorNumber = DetectorNumberX;
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::PreTreat() {
+ // This method typically applies thresholds and calibrations
+ // Might test for disabled channels for more complex detector
+
+ // clear pre-treated object
+ ClearPreTreatedData();
+ unsigned int mysizeX = m_EventData->GetXMultEnergy();
+ unsigned int mysizeY = m_EventData->GetYMultEnergy();
+
+ // instantiate CalibrationManager
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+
+ if (mysizeX == mysizeY) {
+ unsigned int mysize = mysizeX;
+ // Energy
+ for (UShort_t i = 0; i < mysize; ++i) {
+ if (m_EventData->GetXE_DetectorNbr(i) == m_EventData->GetYE_DetectorNbr(i)) {
+ if (m_EventData->GetX_Energy(i) > m_E_RAW_Threshold) {
+ double x = (m_EventData->GetXE_StripNbr(i) - 16.5) * 2, y = (m_EventData->GetYE_StripNbr(i) - 16.5) * 2;
+ // int dist = std::round(TMath::Sqrt(x*x+y*y+4)-2);
+ int dist = std::round(TMath::Sqrt(y * y + 4) - 2);
+ // if(dist > 38) dist = 38;
+ string name = "SEASON/D";
+ name += NPL::itoa(m_EventData->GetXE_DetectorNbr(i));
+ name += "_DIST";
+ name += NPL::itoa(dist);
+ name += "_ENERGY";
+ string namestrip = "SEASON/D";
+ namestrip += NPL::itoa(m_EventData->GetXE_DetectorNbr(i));
+ namestrip += "_STRIPX";
+ namestrip += NPL::itoa(m_EventData->GetXE_StripNbr(i));
+ namestrip += "_ENERGY";
+ // cout << name << endl;
+ Double_t Estrip = Cal->ApplyCalibration(namestrip, m_EventData->GetX_Energy(i));
+ Double_t E = Cal->ApplyCalibration(name, Estrip);
+ if (E > m_XE_Threshold) {
+ m_PreTreatedData->SetEnergy(m_EventData->GetXE_DetectorNbr(i), m_EventData->GetXE_StripNbr(i),
+ m_EventData->GetYE_StripNbr(i), E);
+ }
+ }
+ }
+ }
+
+ // Time
+ for (UShort_t i = 0; i < mysize; ++i) {
+ if (m_EventData->GetXT_DetectorNbr(i) == m_EventData->GetYT_DetectorNbr(i)) {
+ Double_t T = Cal->ApplyCalibration("SEASON/D" + NPL::itoa(m_EventData->GetXT_DetectorNbr(i)) + "_STRIPX" +
+ NPL::itoa(m_EventData->GetXT_StripNbr(i)) + "_STRIPY" +
+ NPL::itoa(m_EventData->GetYT_StripNbr(i)) + "_Time",
+ m_EventData->GetX_Time(i));
+ m_PreTreatedData->SetTime(m_EventData->GetXT_DetectorNbr(i), m_EventData->GetXT_StripNbr(i),
+ m_EventData->GetYT_StripNbr(i), T);
+ }
+ }
+ }
+
+ // X Energy
+ for (UShort_t i = 0; i < mysizeX; ++i) {
+ if (m_EventData->GetX_Energy(i) > m_E_RAW_Threshold) {
+ Double_t EX = Cal->ApplyCalibration("SEASON/D" + NPL::itoa(m_EventData->GetXE_DetectorNbr(i)) + "_STRIP" +
+ NPL::itoa(m_EventData->GetXE_StripNbr(i)) + "_Energy",
+ m_EventData->GetX_Energy(i));
+ if (EX > m_XE_Threshold) {
+ m_PreTreatedData->SetXEnergy(m_EventData->GetXE_DetectorNbr(i), m_EventData->GetXE_StripNbr(i), EX);
+ }
+ }
+ }
+
+ // X Time
+ for (UShort_t i = 0; i < mysizeX; ++i) {
+ Double_t TX = Cal->ApplyCalibration("SEASON/D" + NPL::itoa(m_EventData->GetXT_DetectorNbr(i)) + "_Strip" +
+ NPL::itoa(m_EventData->GetXT_StripNbr(i)) + "_Time",
+ m_EventData->GetX_Time(i));
+ m_PreTreatedData->SetXTime(m_EventData->GetXT_DetectorNbr(i), m_EventData->GetXT_StripNbr(i), TX);
+ }
+
+ // Y Energy
+ for (UShort_t i = 0; i < mysizeY; ++i) {
+ if (m_EventData->GetY_Energy(i) > m_E_RAW_Threshold) {
+ Double_t EY = Cal->ApplyCalibration("SEASON/D" + NPL::itoa(m_EventData->GetYE_DetectorNbr(i)) + "_STRIP" +
+ NPL::itoa(m_EventData->GetYE_StripNbr(i)) + "_Energy",
+ m_EventData->GetY_Energy(i));
+ if (EY > m_YE_Threshold) {
+ m_PreTreatedData->SetYEnergy(m_EventData->GetYE_DetectorNbr(i), m_EventData->GetYE_StripNbr(i), EY);
+ }
+ }
+ }
+
+ // Y Time
+ for (UShort_t i = 0; i < mysizeY; ++i) {
+ Double_t TY = Cal->ApplyCalibration("SEASON/D" + NPL::itoa(m_EventData->GetYT_DetectorNbr(i)) + "_Strip" +
+ NPL::itoa(m_EventData->GetYT_StripNbr(i)) + "_Time",
+ m_EventData->GetY_Time(i));
+ m_PreTreatedData->SetYTime(m_EventData->GetYT_DetectorNbr(i), m_EventData->GetYT_StripNbr(i), TY);
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::ReadAnalysisConfig() {
+ bool ReadingStatus = false;
+
+ // path to file
+ string FileName = "./configs/ConfigSEASON.dat";
+
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigSEASON.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameter for Analysis from ConfigSEASON.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigSEASON.dat %%%");
+ asciiConfig->Append(FileName.c_str());
+ asciiConfig->AppendLine("");
+ // read analysis config file
+ string LineBuffer, DataBuffer, whatToDo;
+ while (!AnalysisConfigFile.eof()) {
+ // Pick-up next line
+ getline(AnalysisConfigFile, LineBuffer);
+
+ // search for "header"
+ string name = "ConfigSEASON";
+ if (LineBuffer.compare(0, name.length(), name) == 0)
+ ReadingStatus = true;
+
+ // loop on tokens and data
+ while (ReadingStatus) {
+ whatToDo = "";
+ AnalysisConfigFile >> whatToDo;
+
+ // Search for comment symbol (%)
+ if (whatToDo.compare(0, 1, "%") == 0) {
+ AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n');
+ }
+ else {
+ ReadingStatus = false;
+ }
+ }
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::Clear() {
+ DetectorNumber.clear();
+ StripNumberX.clear();
+ StripNumberY.clear();
+ StripX.clear();
+ StripY.clear();
+ Energy.clear();
+ Time.clear();
+
+ DetectorNumberX.clear();
+ DetectorNumberY.clear();
+ EnergyX.clear();
+ EnergyY.clear();
+ TimeX.clear();
+ TimeY.clear();
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::ReadConfiguration(NPL::InputParser parser) {
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SEASON");
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks.size() << " detectors found " << endl;
+
+ vector<string> cart = {"X1_Y1", "X1_YMax", "XMax_Y1", "XMax_YMax", "NStripsX", "NStripsY", "Group"};
+
+ for (unsigned int i = 0; i < blocks.size(); i++) {
+ if (blocks[i]->HasTokenList(cart)) {
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// SEASON " << i + 1 << endl;
+
+ TVector3 A = blocks[i]->GetTVector3("X1_Y1", "mm");
+ TVector3 B = blocks[i]->GetTVector3("X1_YMax", "mm");
+ TVector3 C = blocks[i]->GetTVector3("XMax_Y1", "mm");
+ TVector3 D = blocks[i]->GetTVector3("XMax_YMax", "mm");
+ int NumberOfStripsX = blocks[i]->GetInt("NStripsX");
+ int NumberOfStripsY = blocks[i]->GetInt("NStripsY");
+ AddDetector(A, B, C, D, NumberOfStripsX, NumberOfStripsY);
+ }
+ else {
+ cout << "ERROR: check your input file formatting " << endl;
+ exit(1);
+ }
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+TVector3 TSEASONPhysics::GetPositionOfInteraction(const int i) const {
+ TVector3 Position;
+ Position = TVector3(GetStripPositionX(DetectorNumber[i], StripNumberX[i], StripNumberY[i]),
+ GetStripPositionY(DetectorNumber[i], StripNumberX[i], StripNumberY[i]),
+ GetStripPositionZ(DetectorNumber[i], StripNumberX[i], StripNumberY[i]));
+ return Position;
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::InitSpectra() {
+ m_Spectra = new TSEASONSpectra(m_NumberOfDetectors, m_NumberOfStripsX, m_NumberOfStripsY);
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::FillSpectra() {
+ m_Spectra->FillRawSpectra(m_EventData);
+ m_Spectra->FillPreTreatedSpectra(m_PreTreatedData);
+ m_Spectra->FillPhysicsSpectra(m_EventPhysics);
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::CheckSpectra() { m_Spectra->CheckSpectra(); }
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::ClearSpectra() {
+ // To be done
+}
+
+///////////////////////////////////////////////////////////////////////////
+map<string, TH1*> TSEASONPhysics::GetSpectra() {
+ if (m_Spectra)
+ return m_Spectra->GetMapHisto();
+ else {
+ map<string, TH1*> empty;
+ return empty;
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+vector<TCanvas*> TSEASONPhysics::GetCanvas() {
+ /* if (m_Spectra)
+ return m_Spectra->GetCanvas();
+ else {*/
+ vector<TCanvas*> empty;
+ return empty;
+ // }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::WriteSpectra() { m_Spectra->WriteSpectra(); }
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::AddParameterToCalibrationManager() {
+ int NumberOfStrips = 32;
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+ vector<double> standard = {0, 1};
+ for (int i = 0; i < m_NumberOfDetectors; ++i) {
+ for (int k = 0; k < 50; k++) {
+ Cal->AddParameter("SEASON", "D" + NPL::itoa(i + 1) + "_DIST" + NPL::itoa(k) + "_ENERGY",
+ "SEASON_D" + NPL::itoa(i + 1) + "_DIST" + NPL::itoa(k) + "_ENERGY", standard);
+ Cal->AddParameter("SEASON", "D" + NPL::itoa(i + 1) + "_FOILDIST" + NPL::itoa(k) + "_ENERGY",
+ "SEASON_D" + NPL::itoa(i + 1) + "_FOILDIST" + NPL::itoa(k) + "_ENERGY", standard);
+ }
+
+ for (int j = 0; j < NumberOfStrips; ++j) {
+ Cal->AddParameter("SEASON", "D" + NPL::itoa(i + 1) + "_STRIPX" + NPL::itoa(j + 1) + "_ENERGY",
+ "SEASON_D" + NPL::itoa(i + 1) + "_STRIPX" + NPL::itoa(j + 1) + "_ENERGY");
+ Cal->AddParameter("SEASON", "D" + NPL::itoa(i + 1) + "_STRIPX" + NPL::itoa(j + 1) + "_TIME",
+ "SEASON_D" + NPL::itoa(i + 1) + "_STRIP" + NPL::itoa(j + 1) + "_TIME");
+ }
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("SEASON", true);
+ inputChain->SetBranchAddress("SEASON", &m_EventData);
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::InitializeRootInputPhysics() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchAddress("SEASON", &m_EventPhysics);
+}
+
+///////////////////////////////////////////////////////////////////////////
+void TSEASONPhysics::InitializeRootOutput() {
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("SEASON", "TSEASONPhysics", &m_EventPhysics);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VDetector* TSEASONPhysics::Construct() { return (NPL::VDetector*)new TSEASONPhysics(); }
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C" {
+class proxy_SEASON {
+ public:
+ proxy_SEASON() {
+ NPL::DetectorFactory::getInstance()->AddToken("SEASON", "SEASON");
+ NPL::DetectorFactory::getInstance()->AddDetector("SEASON", TSEASONPhysics::Construct);
+ }
+};
+
+proxy_SEASON p_SEASON;
+}
diff --git a/NPLib/Utility/npreader.cxx b/NPLib/Utility/npreader.cxx
index 6a7fd709300569945c17f421aea9faf83884d8a5..670091b247a57be18838d506e5b4f3db5b8c494b 100644
--- a/NPLib/Utility/npreader.cxx
+++ b/NPLib/Utility/npreader.cxx
@@ -198,15 +198,10 @@ int main(int argc , char** argv){
if(!IsPhysics){
std::cout << "test npa branch 2" << std::endl;
while(inputTreeReader->Next()){
- //for (unsigned long i = first_entry ; i < nentries + first_entry; i++) {
- // Get the raw Data
- // Chain -> GetEntry(i);
+
// Build the current event
-
- //myDetector->GetRawData();
if(UserAnalysis->UnallocateBeforeBuild()){
myDetector->BuildPhysicalEvent();
- //std::cout << "test" << std::endl;
// User Analysis
if(UserAnalysis->UnallocateBeforeTreat()){
UserAnalysis->TreatEvent();
@@ -216,28 +211,7 @@ int main(int argc , char** argv){
}
}
current_tree = Chain->GetTreeNumber()+1;
- // std::cout << "Display :" << treated << " " << nentries << " " << mean_rate << std::endl;
- //ProgressDisplay(begin,end,treated,inter,nentries,mean_rate,displayed,current_tree,total_tree);
ProgressDisplay(tv_begin,tv_end,treated,inter,nentries,mean_rate,displayed,current_tree,total_tree);
- //if(myOptionManager->GetOnline() && i%10000==0){
- // myDetector->CheckSpectraServer();
- // bool first = true;
- // while(!Chain || first){
- // first = false;
- // RootInput::getInstance()->GetFile()->ReadKeys(kTRUE);
-
- // Chain = (TChain*) RootInput::getInstance()->GetFile()->FindKeyAny(ChainName)->ReadObj();
- // new_nentries = Chain->GetEntries();
- // if(new_nentries!=nentries){
- // RootInput::getInstance()->SetChain(Chain);
- // myDetector->InitializeRootInput();
- // nentries = Chain->GetEntries();
- // }
- // else{
- // first = true;
- // }
- // }
- //}
}
}
diff --git a/NPSimulation/Detectors/Actar/Actar.cc b/NPSimulation/Detectors/Actar/Actar.cc
index 3691d550d311eea3036da6ab4e4e122a8ce8e284..f41979c9a448f080c95826f49dc0e4c2abb7cf21 100644
--- a/NPSimulation/Detectors/Actar/Actar.cc
+++ b/NPSimulation/Detectors/Actar/Actar.cc
@@ -20,493 +20,499 @@
*****************************************************************************/
// C++ headers
-#include <sstream>
#include <cmath>
#include <limits>
-//G4 Geometry object
-#include "G4Tubs.hh"
+#include <sstream>
+// G4 Geometry object
#include "G4Box.hh"
+#include "G4Tubs.hh"
-//G4 sensitive
-#include "G4SDManager.hh"
+// G4 sensitive
#include "G4MultiFunctionalDetector.hh"
+#include "G4SDManager.hh"
-//G4 various object
+// G4 various object
+#include "G4Colour.hh"
#include "G4Material.hh"
-#include "G4Transform3D.hh"
#include "G4PVPlacement.hh"
+#include "G4Transform3D.hh"
#include "G4VisAttributes.hh"
-#include "G4Colour.hh"
// G4 Field
-#include "G4FieldManager.hh"
+#include "G4ChordFinder.hh"
+#include "G4ClassicalRK4.hh"
#include "G4ElectricField.hh"
-#include "G4UniformElectricField.hh"
-#include "G4TransportationManager.hh"
#include "G4EqMagElectricField.hh"
-#include "G4MagIntegratorStepper.hh"
-#include "G4ClassicalRK4.hh"
+#include "G4FieldManager.hh"
#include "G4MagIntegratorDriver.hh"
-#include "G4ChordFinder.hh"
+#include "G4MagIntegratorStepper.hh"
#include "G4MaterialPropertiesTable.hh"
+#include "G4TransportationManager.hh"
+#include "G4UniformElectricField.hh"
// NPTool header
#include "Actar.hh"
+#include "CalorimeterScorers.hh"
#include "DSSDScorers.hh"
-#include "TPCScorers.hh"
-#include "RootOutput.h"
+#include "FastDriftElectron.hh"
#include "MaterialManager.hh"
-#include "NPSDetectorFactory.hh"
#include "NPOptionManager.h"
-#include "FastDriftElectron.hh"
+#include "NPSDetectorFactory.hh"
#include "NPSHitsMap.hh"
-#include "CalorimeterScorers.hh"
-
+#include "RootOutput.h"
+#include "TPCScorers.hh"
// CLHEP header
#include "CLHEP/Random/RandGauss.h"
// ROOT
-#include "TH1D.h"
#include "TF1.h"
+#include "TH1D.h"
using namespace std;
using namespace CLHEP;
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-namespace Actar_NS{
- // Energy and time Resolution
- const double ChargeThreshold = 0;
- //const double ResoTime = 0.1*ns ;
- const double ResoCharge = 5./100 ;
- const double ChamberThickness = 376*mm ;
- const double ChamberWidth = 376*mm ;
- const double ChamberHeight = 40*cm ;
- //const int NumberOfPads = 16384;
- //const int PadX = 128;
- //const int PadZ = 128;
-
- const double Nose_Rmin = 2.5*cm;
- const double Nose_Rmax = 3.5*cm;
- const double Nose_Length = 12*cm;
-
- const double Mylar_Rmax = 3.5*cm;
- const double Mylar_Thickness = 7*micrometer;
-
- const double XGazVolume = 256.*mm;
- const double YGazVolume = 256.*mm;
- const double ZGazVolume = 256.*mm;
-
- const double SiliconHeight = 53.*mm;
- const double SiliconWidth = 53.*mm;
- const double SiliconThickness = 0.7*mm;
- const double DistInterSi = 1.*mm;
- const double Si_PosZ=175.*mm;
- const double ResoSilicon = 0.60/2.35;
- const double EnergyThreshold = 0.1;
-
- const double VamosSiliconHeight = 70*mm;
- const double VamosSiliconWidth = 50*mm;
- const double VamosSiliconThickness = 0.5*mm;
- const double VamosSiliconDistanInterSi = 1*mm;
- const double VamosSilicon_PosZ = -160*mm;
-
- const double CsIThickness = 1.*cm;
- const double CsIHeight = 2.5*cm;
- const double CsIWidth = 2.5*cm;
- const double DistInterCsI = 1.*mm;
- const double CsI_PosZ = 16.*cm;
- //const double ResoCsI = 0.200/2.35;
-
- const double BeamDumpRadius = 30*mm;
- const double BeamDumpThickness = 5*mm;
- const double BeamDump_PosZ = 160*mm;
-
-}
+namespace Actar_NS {
+ // Energy and time Resolution
+ const double ChargeThreshold = 0;
+ // const double ResoTime = 0.1*ns ;
+ const double ResoCharge = 5. / 100;
+ const double ChamberThickness = 376 * mm;
+ const double ChamberWidth = 376 * mm;
+ const double ChamberHeight = 40 * cm;
+ // const int NumberOfPads = 16384;
+ // const int PadX = 128;
+ // const int PadZ = 128;
+
+ const double Nose_Rmin = 2.5 * cm;
+ const double Nose_Rmax = 3.5 * cm;
+ const double Nose_Length = 12 * cm;
+
+ const double Mylar_Rmax = 3.5 * cm;
+ const double Mylar_Thickness = 7 * micrometer;
+
+ const double XGazVolume = 256. * mm;
+ const double YGazVolume = 256. * mm;
+ const double ZGazVolume = 256. * mm;
+
+ const double SiliconHeight = 53. * mm;
+ const double SiliconWidth = 53. * mm;
+ const double SiliconThickness = 0.7 * mm;
+ const double DistInterSi = 1. * mm;
+ const double Si_PosZ = 175. * mm;
+ const double ResoSilicon = 0.60 / 2.35;
+ const double EnergyThreshold = 0.1;
+
+ const double VamosSiliconHeight = 70 * mm;
+ const double VamosSiliconWidth = 50 * mm;
+ const double VamosSiliconThickness = 0.5 * mm;
+ const double VamosSiliconDistanInterSi = 1 * mm;
+ const double VamosSilicon_PosZ = -160 * mm;
+
+ const double CsIThickness = 1. * cm;
+ const double CsIHeight = 2.5 * cm;
+ const double CsIWidth = 2.5 * cm;
+ const double DistInterCsI = 1. * mm;
+ const double CsI_PosZ = 16. * cm;
+ // const double ResoCsI = 0.200/2.35;
+
+ const double BeamDumpRadius = 30 * mm;
+ const double BeamDumpThickness = 5 * mm;
+ const double BeamDump_PosZ = 160 * mm;
+
+} // namespace Actar_NS
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Actar Specific Method
-Actar::Actar(){
- m_Event = new TActarData() ;
- m_EventReduced = new MEventReduced();
- m_ActarScorer = 0;
- m_SquareDetector = 0;
-
- // RGB Color + Transparency
- m_VisChamber = new G4VisAttributes(G4Colour(0.7, 0.7, 0.7, 0.3));
- m_VisWindows = new G4VisAttributes(G4Colour(1, 0, 0, 0.25));
- m_VisGas = new G4VisAttributes(G4Colour(0, 0.5, 0.5, 0.3));
- m_VisPads = new G4VisAttributes(G4Colour(255, 223, 50, 0.8));
- m_VisMicromegas = new G4VisAttributes(G4Colour(100, 100, 100, 0.4));
- m_SiliconVisAtt = new G4VisAttributes(G4Colour(0.529412, 0.807843, 0.980392, 0.95)) ;
- m_CsIVisAtt = new G4VisAttributes(G4Colour(0.429412, 0.607843, 0.780392, 0.95));
- m_BeamDumpVisAtt = new G4VisAttributes(G4Colour(0.9, 0.5, 0.5));
- m_VisPads->SetForceWireframe(true);
-
- m_build_BeamDump= 0;
- m_build_Silicon=1;
- m_build_Vamos_Silicon=0;
- m_build_CsI=1;
- m_ReactionRegion=NULL;
- // Lookup table //
- // bool ReadingLookupTable = false;
- // Opening the LookUp Table LT.dat
- G4String GlobalPath = getenv("NPTOOL");
- G4String LT_FileName = GlobalPath + "/NPSimulation/Detectors/Actar/LT.dat";
- //string LT_FileName = "./Detectors/Actar/LT.dat";
- //string LT_FileName = "./configs/LT.dat";
- ifstream LTConfigFile;
- LTConfigFile.open(LT_FileName.c_str());
- if(!LTConfigFile.is_open()){
- cout << "No Lookup Table in " << LT_FileName << " found!" << endl;
- return;
+Actar::Actar() {
+ m_Event = new TActarData();
+ m_EventReduced = new MEventReduced();
+ m_ActarScorer = 0;
+ m_SquareDetector = 0;
+
+ // RGB Color + Transparency
+ m_VisChamber = new G4VisAttributes(G4Colour(0.7, 0.7, 0.7, 0.3));
+ m_VisWindows = new G4VisAttributes(G4Colour(1, 0, 0, 0.25));
+ m_VisGas = new G4VisAttributes(G4Colour(0, 0.5, 0.5, 0.3));
+ m_VisPads = new G4VisAttributes(G4Colour(255, 223, 50, 0.8));
+ m_VisMicromegas = new G4VisAttributes(G4Colour(100, 100, 100, 0.4));
+ m_SiliconVisAtt = new G4VisAttributes(G4Colour(0.529412, 0.807843, 0.980392, 0.95));
+ m_CsIVisAtt = new G4VisAttributes(G4Colour(0.429412, 0.607843, 0.780392, 0.95));
+ m_BeamDumpVisAtt = new G4VisAttributes(G4Colour(0.9, 0.5, 0.5));
+ m_VisPads->SetForceWireframe(true);
+
+ m_build_BeamDump = 0;
+ m_build_Silicon = 1;
+ m_build_Vamos_Silicon = 0;
+ m_build_CsI = 1;
+ m_ReactionRegion = NULL;
+ // Lookup table //
+ // bool ReadingLookupTable = false;
+ // Opening the LookUp Table LT.dat
+ G4String GlobalPath = getenv("NPTOOL");
+ G4String LT_FileName = GlobalPath + "/NPSimulation/Detectors/Actar/LT.dat";
+ // string LT_FileName = "./Detectors/Actar/LT.dat";
+ // string LT_FileName = "./configs/LT.dat";
+ ifstream LTConfigFile;
+ LTConfigFile.open(LT_FileName.c_str());
+ if (!LTConfigFile.is_open()) {
+ cout << "No Lookup Table in " << LT_FileName << " found!" << endl;
+ return;
+ }
+ else {
+ cout << "/// Using LookupTable from: " << LT_FileName << " ///" << endl;
+ int co, as, ag, ch;
+ int pX, pY;
+ for (int i = 0; i < NumberOfCobo * NumberOfASAD * NumberOfAGET * NumberOfChannel; i++) {
+ LTConfigFile >> co >> as >> ag >> ch >> pX >> pY;
+ if (pX != -1 && pY != -1) {
+ m_PadToCobo[pX][pY] = co;
+ m_PadToAsad[pX][pY] = as;
+ m_PadToAGET[pX][pY] = ag;
+ m_PadToChannel[pX][pY] = ch;
+ }
}
- else{
- cout << "/// Using LookupTable from: " << LT_FileName << " ///" << endl;
- int co, as, ag, ch;
- int pX, pY;
- for(int i=0;i<NumberOfCobo*NumberOfASAD*NumberOfAGET*NumberOfChannel;i++){
- LTConfigFile >> co >> as >> ag >> ch >> pX >> pY;
- if(pX!=-1 && pY !=-1){
- m_PadToCobo[pX][pY] = co;
- m_PadToAsad[pX][pY] = as;
- m_PadToAGET[pX][pY] = ag;
- m_PadToChannel[pX][pY] = ch;
- }
- }
- //ReadingLookupTable = true;
- }
- LTConfigFile.close();
-
+ // ReadingLookupTable = true;
+ }
+ LTConfigFile.close();
}
-Actar::~Actar(){
-}
+Actar::~Actar() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Actar::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);
+void Actar::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 Actar::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);
+void Actar::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* Actar::BuildDetector(){
+G4LogicalVolume* Actar::BuildDetector() {
- G4Material* Cu= MaterialManager::getInstance()->GetMaterialFromLibrary("Cu");
- G4Material* Si= MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
- G4Material* Al= MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
- G4Material* Mylar= MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
- G4Material* MaterialCsI = MaterialManager::getInstance()->GetMaterialFromLibrary("CsI");
+ G4Material* Cu = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu");
+ G4Material* Si = MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
+ G4Material* Al = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ G4Material* Mylar = MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
+ G4Material* MaterialCsI = MaterialManager::getInstance()->GetMaterialFromLibrary("CsI");
- if(!m_SquareDetector){
- // Main volume
- G4Box* sChamber = new G4Box("Actar_Box",Actar_NS::ChamberWidth*0.5,
- Actar_NS::ChamberHeight*0.5,Actar_NS::ChamberThickness*0.5);
+ if (!m_SquareDetector) {
+ // Main volume
+ G4Box* sChamber = new G4Box("Actar_Box", Actar_NS::ChamberWidth * 0.5, Actar_NS::ChamberHeight * 0.5,
+ Actar_NS::ChamberThickness * 0.5);
- //Nose volume
- G4Tubs* sNose = new G4Tubs("Actar_Nose",Actar_NS::Nose_Rmin, Actar_NS::Nose_Rmax,Actar_NS::Nose_Length*0.5,
- 0*deg, 360*deg);
+ // Nose volume
+ G4Tubs* sNose = new G4Tubs("Actar_Nose", Actar_NS::Nose_Rmin, Actar_NS::Nose_Rmax, Actar_NS::Nose_Length * 0.5,
+ 0 * deg, 360 * deg);
- //Mylar volume
- G4Tubs* sWindows = new G4Tubs("Actar_Windows",0, Actar_NS::Mylar_Rmax,Actar_NS::Mylar_Thickness*0.5,
- 0*deg, 360*deg);
+ // Mylar volume
+ G4Tubs* sWindows =
+ new G4Tubs("Actar_Windows", 0, Actar_NS::Mylar_Rmax, Actar_NS::Mylar_Thickness * 0.5, 0 * deg, 360 * deg);
- // Cage volume
- G4Box* sCage = new G4Box("Actar_Gas",Actar_NS::XGazVolume*0.5,
- Actar_NS::YGazVolume*0.5,Actar_NS::ZGazVolume*0.5);
+ // Cage volume
+ G4Box* sCage =
+ new G4Box("Actar_Gas", Actar_NS::XGazVolume * 0.5, Actar_NS::YGazVolume * 0.5, Actar_NS::ZGazVolume * 0.5);
- // Pad
- G4Box* sPad = new G4Box("Actar_Pad",256*mm*0.5,
- 2*mm*0.5,256*mm*0.5);
+ // Pad
+ G4Box* sPad = new G4Box("Actar_Pad", 256 * mm * 0.5, 2 * mm * 0.5, 256 * mm * 0.5);
- // Micromegas
- G4Box* sMicromegas = new G4Box("Actar_Micromegas",256*mm*0.5,
- 220*micrometer*0.5,256*mm*0.5);
+ // Micromegas
+ G4Box* sMicromegas = new G4Box("Actar_Micromegas", 256 * mm * 0.5, 220 * micrometer * 0.5, 256 * mm * 0.5);
- // Cathode
- G4Box* sCathode = new G4Box("Actar_Cathode",26.5*cm*0.5,
- 1*um*0.5,25.6*cm*0.5);
+ // Cathode
+ G4Box* sCathode = new G4Box("Actar_Cathode", 26.5 * cm * 0.5, 1 * um * 0.5, 25.6 * cm * 0.5);
+ unsigned const int NumberOfGasMix = m_GasMaterial.size();
+ double density = 0;
+ double density_sum = 0;
+ vector<G4Material*> GasComponent;
+ vector<double> FractionMass;
- unsigned const int NumberOfGasMix = m_GasMaterial.size();
-
- double density=0;
- double density_sum=0;
- vector<G4Material*> GasComponent;
- vector<double> FractionMass;
-
- for(unsigned int i=0; i<NumberOfGasMix; i++){
- GasComponent.push_back(MaterialManager::getInstance()->GetGasFromLibrary(m_GasMaterial[i],m_Pressure,m_Temperature) );
- }
- for(unsigned int i=0; i<NumberOfGasMix; i++){
- density += ((double)m_GasFraction[i]/100)*GasComponent[i]->GetDensity();
- density_sum += GasComponent[i]->GetDensity();
- }
- //cout << "density = " << density*cm3/g << endl;
-
- for(unsigned int i=0; i<NumberOfGasMix; i++){
- FractionMass.push_back(GasComponent[i]->GetDensity()/density_sum);
- }
-
- G4Material* GasMaterial = new G4Material("GasMix", density, NumberOfGasMix, kStateGas, m_Temperature, m_Pressure);
- G4Material* DriftGasMaterial = new G4Material("DriftGasMix", density, NumberOfGasMix, kStateGas, m_Temperature, m_Pressure);
-
- for(unsigned int i=0; i<NumberOfGasMix; i++){
- GasMaterial->AddMaterial(GasComponent[i], FractionMass[i]);
- DriftGasMaterial->AddMaterial(GasComponent[i], FractionMass[i]);
- cout << GasComponent[i] << endl;
- }
-
- G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
- MPT->AddConstProperty("DE_PAIRENERGY",20*eV);
- MPT->AddConstProperty("DE_YIELD",3e-1);
- //MPT->AddConstProperty("DE_AMPLIFICATION",2);
- MPT->AddConstProperty("DE_ABSLENGTH",1*pc);
- MPT->AddConstProperty("DE_DRIFTSPEED",0.8*cm/microsecond);
- MPT->AddConstProperty("DE_TRANSVERSALSPREAD",2e-5*mm2/ns);
- MPT->AddConstProperty("DE_LONGITUDINALSPREAD",4e-6*mm2/ns);
-
- DriftGasMaterial->SetMaterialPropertiesTable(MPT);
-
- G4MaterialPropertiesTable* MPT2 = new G4MaterialPropertiesTable();
- MPT2->AddConstProperty("DE_AMPLIFICATION",1000);
- MPT2->AddConstProperty("DE_ABSLENGTH",1*pc);
-
- Al->SetMaterialPropertiesTable(MPT2);
-
- m_SquareDetector = new G4LogicalVolume(sChamber,GasMaterial,"logic_Actar_Box",0,0,0);
- m_logicGas = new G4LogicalVolume(sCage,DriftGasMaterial,"logic_Gas",0,0,0);
- G4LogicalVolume* logicPad = new G4LogicalVolume(sPad,Cu,"logic_Pad",0,0,0);
- G4LogicalVolume* logicMicromegas = new G4LogicalVolume(sMicromegas,Al,"logic_Micromegas",0,0,0);
-
- G4LogicalVolume* logicNose = new G4LogicalVolume(sNose,Cu,"logic_Nose",0,0,0);
- /*G4LogicalVolume* logicCathode = */ new G4LogicalVolume(sCathode,Cu,"logic_Cathode",0,0,0);
- G4LogicalVolume* logicWindows = new G4LogicalVolume(sWindows,Mylar,"logic_Windows",0,0,0);
-
- G4RotationMatrix* Rot = new G4RotationMatrix();
- //new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Nose_Length*0.5)),
- // logicNose,
- // "ActarNose",m_SquareDetector,false,0);
-
- //new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Mylar_Thickness*0.5+Actar_NS::Nose_Length)),
- // logicWindows,
- // "ActarEntranceWindows",m_SquareDetector,false,0);
-
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,0)),
- m_logicGas,
- "ActarGas",m_SquareDetector,false,0);
-
- //new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,Actar_NS::YGazVolume*0.5-0.3*cm,0)),
- // logicMicromegas,
- // "ActarMicromegas",m_logicGas,false,0);
-
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,Actar_NS::YGazVolume*0.5,0)),
- logicPad,
- "ActarPad",m_logicGas,false,0);
-
- /*int pad=0;
- m_PadToXRow.clear();
- m_PadToZColumn.clear();
- for(int i=0; i<Actar_NS::PadX; i++){
- for(int j=0; j<Actar_NS::PadZ; j++){
- m_PadToXRow[pad] = i;
- m_PadToZColumn[pad] = j;
- double X=(i-64)*2*mm;
- double Z=(j-64)*2*mm;
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(X,Actar_NS::YGazVolume*0.5,Z)),
- logicPad,
- "ActarPad",m_logicGas,false,pad+1);
-
- pad++;
- }
- }*/
-
- /*new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(3*cm-0.5*1*um,0,0)),
- logicCathode,
- "ActarCathode",m_logicGas,false,0);
-
-
-
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-6*cm+6*micrometer)),
- logicWindows,
- "ActarEntranceWindows",m_SquareDetector,false,0);*/
-
- G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0,-400*volt/cm,0.0));
- // Create an equation of motion for this field
- G4EqMagElectricField* Equation = new G4EqMagElectricField(field);
- G4MagIntegratorStepper* Stepper = new G4ClassicalRK4( Equation, 8 );
-
- // Get the global field manager
- G4FieldManager* FieldManager= new G4FieldManager();
- // Set this field to the global field manager
- FieldManager->SetDetectorField(field);
- m_logicGas->SetFieldManager(FieldManager,true);
-
- G4MagInt_Driver* IntgrDriver = new G4MagInt_Driver(0.1*mm,
- Stepper,
- Stepper->GetNumberOfVariables() );
-
- G4ChordFinder* ChordFinder = new G4ChordFinder(IntgrDriver);
- FieldManager->SetChordFinder( ChordFinder );
-
-
- logicPad->SetSensitiveDetector(m_ActarScorer);
- logicNose->SetVisAttributes(m_VisChamber);
- m_SquareDetector->SetVisAttributes(m_VisChamber);
- m_logicGas->SetVisAttributes(m_VisGas);
- logicWindows->SetVisAttributes(m_VisWindows);
- logicPad->SetVisAttributes(m_VisPads);
- logicMicromegas->SetVisAttributes(m_VisMicromegas);
- //m_SquareDetector->SetSensitiveDetector(m_ActarScorer);
+ for (unsigned int i = 0; i < NumberOfGasMix; i++) {
+ GasComponent.push_back(
+ MaterialManager::getInstance()->GetGasFromLibrary(m_GasMaterial[i], m_Pressure, m_Temperature));
}
+ for (unsigned int i = 0; i < NumberOfGasMix; i++) {
+ density += ((double)m_GasFraction[i] / 100) * GasComponent[i]->GetDensity();
+ density_sum += GasComponent[i]->GetDensity();
+ }
+ // cout << "density = " << density*cm3/g << endl;
-
- ////////////////////////////////////////////////////
- ///////////////////// Beam Dump ////////////////////
- ////////////////////////////////////////////////////
- if(m_build_BeamDump){
- G4Tubs* sBeamDump = new G4Tubs("Actar_BeamDump",0*mm, Actar_NS::BeamDumpRadius,Actar_NS::BeamDumpThickness*0.5,
- 0*deg, 360*deg);
- m_LogicBeamDump = new G4LogicalVolume(sBeamDump, Cu, "logicBeamDump",0,0,0);
-
- G4ThreeVector positionBeamDump = G4ThreeVector(0, 0, Actar_NS::BeamDump_PosZ);
- new G4PVPlacement(new G4RotationMatrix(0,0,0),
- positionBeamDump,
- m_LogicBeamDump,"BeamDump",
- m_SquareDetector,false,0);
-
- m_LogicBeamDump->SetVisAttributes(m_BeamDumpVisAtt) ;
-
+ for (unsigned int i = 0; i < NumberOfGasMix; i++) {
+ FractionMass.push_back(GasComponent[i]->GetDensity() / density_sum);
}
- ///////////////////////////////////////////////////
- ///////////////////// Thin Si /////////////////////
- ///////////////////////////////////////////////////
- int SiliconNumber=0;
- if(m_build_Silicon){
- G4Box* solidSi = new G4Box("Si", 0.5*Actar_NS::SiliconWidth, 0.5*Actar_NS::SiliconHeight, 0.5*Actar_NS::SiliconThickness); ;
- m_LogicSilicon = new G4LogicalVolume(solidSi, Si, "logicSi", 0, 0, 0);
-
- for(int k=0;k<4; k++){
- for(int p=0; p<5; p++){
- double PosX;
- double PosY;
- if(k==0) PosY= -1.5*Actar_NS::SiliconHeight-2*Actar_NS::DistInterSi;
- if(k==1) PosY= -0.5*Actar_NS::SiliconHeight-1*Actar_NS::DistInterSi;
- if(k==2) PosY= 0.5*Actar_NS::SiliconHeight+1*Actar_NS::DistInterSi;
- if(k==3) PosY= 1.5*Actar_NS::SiliconHeight+2*Actar_NS::DistInterSi;
- if(p==0) PosX= -2*Actar_NS::SiliconWidth-2*Actar_NS::DistInterSi;
- if(p==1) PosX= -1*Actar_NS::SiliconWidth-1*Actar_NS::DistInterSi;
- if(p==2) PosX= 0;
- if(p==3) PosX= 1*Actar_NS::SiliconWidth+2*Actar_NS::DistInterSi;
- if(p==4) PosX= 2*Actar_NS::SiliconWidth+1*Actar_NS::DistInterSi;
-
- G4ThreeVector positionSi = G4ThreeVector(PosX, PosY, Actar_NS::Si_PosZ);
- new G4PVPlacement(new G4RotationMatrix(0,0,0),
- positionSi,
- m_LogicSilicon,"Si",
- m_SquareDetector,false,SiliconNumber);
- SiliconNumber++;
- }
- }
- // Set Si sensible
- m_LogicSilicon->SetSensitiveDetector(m_SiliconScorer);
+ G4Material* GasMaterial = new G4Material("GasMix", density, NumberOfGasMix, kStateGas, m_Temperature, m_Pressure);
+ G4Material* DriftGasMaterial =
+ new G4Material("DriftGasMix", density, NumberOfGasMix, kStateGas, m_Temperature, m_Pressure);
- // Visualisation of ThinSi
- m_LogicSilicon->SetVisAttributes(m_SiliconVisAtt);
+ for (unsigned int i = 0; i < NumberOfGasMix; i++) {
+ GasMaterial->AddMaterial(GasComponent[i], FractionMass[i]);
+ DriftGasMaterial->AddMaterial(GasComponent[i], FractionMass[i]);
+ cout << GasComponent[i] << endl;
}
- ////////////////////////////////////////////////////
- ///////////////////// Vamos Si /////////////////////
- ////////////////////////////////////////////////////
- if(m_build_Vamos_Silicon){
- G4Box* solidSi = new G4Box("Si", 0.5*Actar_NS::VamosSiliconWidth, 0.5*Actar_NS::VamosSiliconHeight, 0.5*Actar_NS::VamosSiliconThickness); ;
- m_LogicVamosSilicon = new G4LogicalVolume(solidSi, Si, "logicVamosSi", 0, 0, 0);
-
-
- int VamosSiliconNumber=0;
- for(int k=0;k<4; k++){
- for(int p=0; p<3; p++){
- double PosX;
- double PosY;
- if(k==0) PosX= -35*mm -0.5*Actar_NS::VamosSiliconWidth - Actar_NS::VamosSiliconDistanInterSi;
- if(k==1) PosX= -35*mm -1.5*Actar_NS::VamosSiliconWidth - Actar_NS::VamosSiliconDistanInterSi;
- if(k==2) PosX= 35*mm + 0.5*Actar_NS::VamosSiliconWidth + Actar_NS::VamosSiliconDistanInterSi;
- if(k==3) PosX= 35*mm + 1.5*Actar_NS::VamosSiliconWidth + Actar_NS::VamosSiliconDistanInterSi;
-
- if(p==0) PosY= -Actar_NS::VamosSiliconHeight-Actar_NS::VamosSiliconDistanInterSi;
- if(p==1) PosY= 0;
- if(p==2) PosY= Actar_NS::VamosSiliconHeight+Actar_NS::VamosSiliconDistanInterSi;
-
- G4ThreeVector positionSi = G4ThreeVector(PosX, PosY, Actar_NS::VamosSilicon_PosZ);
- new G4PVPlacement(new G4RotationMatrix(0,0,0),
- positionSi,
- m_LogicVamosSilicon,"Si",
- m_SquareDetector,false,SiliconNumber);
- VamosSiliconNumber++;
- SiliconNumber++;
- }
- }
+ G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
+ MPT->AddConstProperty("DE_PAIRENERGY", 20 * eV, true);
+ MPT->AddConstProperty("DE_YIELD", 3e-1, true);
+ // MPT->AddConstProperty("DE_AMPLIFICATION",2,true);
+ MPT->AddConstProperty("DE_ABSLENGTH", 1 * pc, true);
+ MPT->AddConstProperty("DE_DRIFTSPEED", 0.8 * cm / microsecond, true);
+ MPT->AddConstProperty("DE_TRANSVERSALSPREAD", 2e-5 * mm2 / ns, true);
+ MPT->AddConstProperty("DE_LONGITUDINALSPREAD", 4e-6 * mm2 / ns, true);
+
+ DriftGasMaterial->SetMaterialPropertiesTable(MPT);
+
+ G4MaterialPropertiesTable* MPT2 = new G4MaterialPropertiesTable();
+ MPT2->AddConstProperty("DE_AMPLIFICATION", 1000, true);
+ MPT2->AddConstProperty("DE_ABSLENGTH", 1 * pc, true);
+
+ Al->SetMaterialPropertiesTable(MPT2);
+
+ m_SquareDetector = new G4LogicalVolume(sChamber, GasMaterial, "logic_Actar_Box", 0, 0, 0);
+ m_logicGas = new G4LogicalVolume(sCage, DriftGasMaterial, "logic_Gas", 0, 0, 0);
+ G4LogicalVolume* logicPad = new G4LogicalVolume(sPad, Cu, "logic_Pad", 0, 0, 0);
+ G4LogicalVolume* logicMicromegas = new G4LogicalVolume(sMicromegas, Al, "logic_Micromegas", 0, 0, 0);
+
+ G4LogicalVolume* logicNose = new G4LogicalVolume(sNose, Cu, "logic_Nose", 0, 0, 0);
+ /*G4LogicalVolume* logicCathode = */ new G4LogicalVolume(sCathode, Cu, "logic_Cathode", 0, 0, 0);
+ G4LogicalVolume* logicWindows = new G4LogicalVolume(sWindows, Mylar, "logic_Windows", 0, 0, 0);
+
+ G4RotationMatrix* Rot = new G4RotationMatrix();
+ // new
+ // G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Nose_Length*0.5)),
+ // logicNose,
+ // "ActarNose",m_SquareDetector,false,0);
+
+ // new
+ // G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-Actar_NS::ChamberThickness*0.5+Actar_NS::Mylar_Thickness*0.5+Actar_NS::Nose_Length)),
+ // logicWindows,
+ // "ActarEntranceWindows",m_SquareDetector,false,0);
+
+ new G4PVPlacement(G4Transform3D(*Rot, G4ThreeVector(0, 0, 0)), m_logicGas, "ActarGas", m_SquareDetector, false, 0);
+
+ // new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,Actar_NS::YGazVolume*0.5-0.3*cm,0)),
+ // logicMicromegas,
+ // "ActarMicromegas",m_logicGas,false,0);
+
+ new G4PVPlacement(G4Transform3D(*Rot, G4ThreeVector(0, Actar_NS::YGazVolume * 0.5, 0)), logicPad, "ActarPad",
+ m_logicGas, false, 0);
+
+ /*int pad=0;
+ m_PadToXRow.clear();
+ m_PadToZColumn.clear();
+ for(int i=0; i<Actar_NS::PadX; i++){
+ for(int j=0; j<Actar_NS::PadZ; j++){
+ m_PadToXRow[pad] = i;
+ m_PadToZColumn[pad] = j;
+ double X=(i-64)*2*mm;
+ double Z=(j-64)*2*mm;
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(X,Actar_NS::YGazVolume*0.5,Z)),
+ logicPad,
+ "ActarPad",m_logicGas,false,pad+1);
+
+ pad++;
+ }
+ }*/
+
+ /*new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(3*cm-0.5*1*um,0,0)),
+ logicCathode,
+ "ActarCathode",m_logicGas,false,0);
+
+
+
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-6*cm+6*micrometer)),
+ logicWindows,
+ "ActarEntranceWindows",m_SquareDetector,false,0);*/
+
+ G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0, -400 * volt / cm, 0.0));
+ // Create an equation of motion for this field
+ G4EqMagElectricField* Equation = new G4EqMagElectricField(field);
+ G4MagIntegratorStepper* Stepper = new G4ClassicalRK4(Equation, 8);
+
+ // Get the global field manager
+ G4FieldManager* FieldManager = new G4FieldManager();
+ // Set this field to the global field manager
+ FieldManager->SetDetectorField(field);
+ m_logicGas->SetFieldManager(FieldManager, true);
+
+ G4MagInt_Driver* IntgrDriver = new G4MagInt_Driver(0.1 * mm, Stepper, Stepper->GetNumberOfVariables());
+
+ G4ChordFinder* ChordFinder = new G4ChordFinder(IntgrDriver);
+ FieldManager->SetChordFinder(ChordFinder);
+
+ logicPad->SetSensitiveDetector(m_ActarScorer);
+ logicNose->SetVisAttributes(m_VisChamber);
+ m_SquareDetector->SetVisAttributes(m_VisChamber);
+ m_logicGas->SetVisAttributes(m_VisGas);
+ logicWindows->SetVisAttributes(m_VisWindows);
+ logicPad->SetVisAttributes(m_VisPads);
+ logicMicromegas->SetVisAttributes(m_VisMicromegas);
+ // m_SquareDetector->SetSensitiveDetector(m_ActarScorer);
+ }
+
+ ////////////////////////////////////////////////////
+ ///////////////////// Beam Dump ////////////////////
+ ////////////////////////////////////////////////////
+ if (m_build_BeamDump) {
+ G4Tubs* sBeamDump = new G4Tubs("Actar_BeamDump", 0 * mm, Actar_NS::BeamDumpRadius,
+ Actar_NS::BeamDumpThickness * 0.5, 0 * deg, 360 * deg);
+ m_LogicBeamDump = new G4LogicalVolume(sBeamDump, Cu, "logicBeamDump", 0, 0, 0);
+
+ G4ThreeVector positionBeamDump = G4ThreeVector(0, 0, Actar_NS::BeamDump_PosZ);
+ new G4PVPlacement(new G4RotationMatrix(0, 0, 0), positionBeamDump, m_LogicBeamDump, "BeamDump", m_SquareDetector,
+ false, 0);
+
+ m_LogicBeamDump->SetVisAttributes(m_BeamDumpVisAtt);
+ }
+ ///////////////////////////////////////////////////
+ ///////////////////// Thin Si /////////////////////
+ ///////////////////////////////////////////////////
+ int SiliconNumber = 0;
+ if (m_build_Silicon) {
+ G4Box* solidSi =
+ new G4Box("Si", 0.5 * Actar_NS::SiliconWidth, 0.5 * Actar_NS::SiliconHeight, 0.5 * Actar_NS::SiliconThickness);
+ ;
+ m_LogicSilicon = new G4LogicalVolume(solidSi, Si, "logicSi", 0, 0, 0);
+
+ for (int k = 0; k < 4; k++) {
+ for (int p = 0; p < 5; p++) {
+ double PosX;
+ double PosY;
+ if (k == 0)
+ PosY = -1.5 * Actar_NS::SiliconHeight - 2 * Actar_NS::DistInterSi;
+ if (k == 1)
+ PosY = -0.5 * Actar_NS::SiliconHeight - 1 * Actar_NS::DistInterSi;
+ if (k == 2)
+ PosY = 0.5 * Actar_NS::SiliconHeight + 1 * Actar_NS::DistInterSi;
+ if (k == 3)
+ PosY = 1.5 * Actar_NS::SiliconHeight + 2 * Actar_NS::DistInterSi;
+ if (p == 0)
+ PosX = -2 * Actar_NS::SiliconWidth - 2 * Actar_NS::DistInterSi;
+ if (p == 1)
+ PosX = -1 * Actar_NS::SiliconWidth - 1 * Actar_NS::DistInterSi;
+ if (p == 2)
+ PosX = 0;
+ if (p == 3)
+ PosX = 1 * Actar_NS::SiliconWidth + 2 * Actar_NS::DistInterSi;
+ if (p == 4)
+ PosX = 2 * Actar_NS::SiliconWidth + 1 * Actar_NS::DistInterSi;
+
+ G4ThreeVector positionSi = G4ThreeVector(PosX, PosY, Actar_NS::Si_PosZ);
+ new G4PVPlacement(new G4RotationMatrix(0, 0, 0), positionSi, m_LogicSilicon, "Si", m_SquareDetector, false,
+ SiliconNumber);
+ SiliconNumber++;
+ }
+ }
- // Set Si sensible
- m_LogicVamosSilicon->SetSensitiveDetector(m_SiliconScorer);
+ // Set Si sensible
+ m_LogicSilicon->SetSensitiveDetector(m_SiliconScorer);
+
+ // Visualisation of ThinSi
+ m_LogicSilicon->SetVisAttributes(m_SiliconVisAtt);
+ }
+
+ ////////////////////////////////////////////////////
+ ///////////////////// Vamos Si /////////////////////
+ ////////////////////////////////////////////////////
+ if (m_build_Vamos_Silicon) {
+ G4Box* solidSi = new G4Box("Si", 0.5 * Actar_NS::VamosSiliconWidth, 0.5 * Actar_NS::VamosSiliconHeight,
+ 0.5 * Actar_NS::VamosSiliconThickness);
+ ;
+ m_LogicVamosSilicon = new G4LogicalVolume(solidSi, Si, "logicVamosSi", 0, 0, 0);
+
+ int VamosSiliconNumber = 0;
+ for (int k = 0; k < 4; k++) {
+ for (int p = 0; p < 3; p++) {
+ double PosX;
+ double PosY;
+ if (k == 0)
+ PosX = -35 * mm - 0.5 * Actar_NS::VamosSiliconWidth - Actar_NS::VamosSiliconDistanInterSi;
+ if (k == 1)
+ PosX = -35 * mm - 1.5 * Actar_NS::VamosSiliconWidth - Actar_NS::VamosSiliconDistanInterSi;
+ if (k == 2)
+ PosX = 35 * mm + 0.5 * Actar_NS::VamosSiliconWidth + Actar_NS::VamosSiliconDistanInterSi;
+ if (k == 3)
+ PosX = 35 * mm + 1.5 * Actar_NS::VamosSiliconWidth + Actar_NS::VamosSiliconDistanInterSi;
+
+ if (p == 0)
+ PosY = -Actar_NS::VamosSiliconHeight - Actar_NS::VamosSiliconDistanInterSi;
+ if (p == 1)
+ PosY = 0;
+ if (p == 2)
+ PosY = Actar_NS::VamosSiliconHeight + Actar_NS::VamosSiliconDistanInterSi;
+
+ G4ThreeVector positionSi = G4ThreeVector(PosX, PosY, Actar_NS::VamosSilicon_PosZ);
+ new G4PVPlacement(new G4RotationMatrix(0, 0, 0), positionSi, m_LogicVamosSilicon, "Si", m_SquareDetector, false,
+ SiliconNumber);
+ VamosSiliconNumber++;
+ SiliconNumber++;
+ }
+ }
- // Visualisation of ThinSi
- m_LogicVamosSilicon->SetVisAttributes(m_SiliconVisAtt);
+ // Set Si sensible
+ m_LogicVamosSilicon->SetSensitiveDetector(m_SiliconScorer);
+
+ // Visualisation of ThinSi
+ m_LogicVamosSilicon->SetVisAttributes(m_SiliconVisAtt);
+ }
+
+ ///////////////////////////////////////////////
+ ///////////////////// CsI /////////////////////
+ ///////////////////////////////////////////////
+ if (m_build_CsI) {
+ G4Box* solidCsI =
+ new G4Box("Si", 0.5 * Actar_NS::CsIWidth, 0.5 * Actar_NS::CsIHeight, 0.5 * Actar_NS::CsIThickness);
+ ;
+ m_LogicCsICrystal = new G4LogicalVolume(solidCsI, MaterialCsI, "logicCsI", 0, 0, 0);
+
+ int CsINumber = 0;
+ for (int k = 0; k < 8; k++) {
+ for (int p = 0; p < 10; p++) {
+ double PosX;
+ double PosY;
+ if (k < 4)
+ PosY = -0.5 * Actar_NS::CsIHeight - 0.5 * Actar_NS::DistInterCsI +
+ (k - 3) * (Actar_NS::CsIHeight + Actar_NS::DistInterCsI);
+ if (k > 3)
+ PosY = 0.5 * Actar_NS::CsIHeight + 0.5 * Actar_NS::DistInterCsI +
+ (k - 4) * (Actar_NS::CsIHeight + Actar_NS::DistInterCsI);
+
+ if (p < 5)
+ PosX = 0.5 * Actar_NS::CsIWidth + 0.5 * Actar_NS::DistInterCsI +
+ (4 - p) * (Actar_NS::CsIWidth + Actar_NS::DistInterCsI);
+ if (p > 4)
+ PosX = -0.5 * Actar_NS::CsIWidth - 0.5 * Actar_NS::DistInterCsI +
+ (5 - p) * (Actar_NS::CsIWidth + Actar_NS::DistInterCsI);
+
+ G4ThreeVector positionCsI = G4ThreeVector(PosX, PosY, Actar_NS::CsI_PosZ);
+ new G4PVPlacement(new G4RotationMatrix(0, 0, 0), positionCsI, m_LogicCsICrystal, "CsI", m_SquareDetector, false,
+ CsINumber);
+ CsINumber++;
+ }
}
- ///////////////////////////////////////////////
- ///////////////////// CsI /////////////////////
- ///////////////////////////////////////////////
- if(m_build_CsI){
- G4Box* solidCsI = new G4Box("Si", 0.5*Actar_NS::CsIWidth, 0.5*Actar_NS::CsIHeight, 0.5*Actar_NS::CsIThickness); ;
- m_LogicCsICrystal = new G4LogicalVolume(solidCsI, MaterialCsI, "logicCsI", 0, 0, 0);
-
- int CsINumber=0;
- for(int k=0;k<8; k++){
- for(int p=0; p<10; p++){
- double PosX;
- double PosY;
- if(k<4) PosY= -0.5*Actar_NS::CsIHeight-0.5*Actar_NS::DistInterCsI+(k-3)*(Actar_NS::CsIHeight+Actar_NS::DistInterCsI);
- if(k>3) PosY= 0.5*Actar_NS::CsIHeight+0.5*Actar_NS::DistInterCsI+(k-4)*(Actar_NS::CsIHeight+Actar_NS::DistInterCsI);
-
- if(p<5) PosX= 0.5*Actar_NS::CsIWidth+0.5*Actar_NS::DistInterCsI+(4-p)*(Actar_NS::CsIWidth+Actar_NS::DistInterCsI);
- if(p>4) PosX= -0.5*Actar_NS::CsIWidth-0.5*Actar_NS::DistInterCsI+(5-p)*(Actar_NS::CsIWidth+Actar_NS::DistInterCsI);
-
- G4ThreeVector positionCsI = G4ThreeVector(PosX, PosY, Actar_NS::CsI_PosZ);
- new G4PVPlacement(new G4RotationMatrix(0,0,0),
- positionCsI,
- m_LogicCsICrystal,"CsI",
- m_SquareDetector,false,CsINumber);
- CsINumber++;
- }
- }
+ // Set Si sensible
+ m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer);
- // Set Si sensible
- m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer);
+ // Visualisation of ThinSi
+ m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt);
+ }
- // Visualisation of ThinSi
- m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt) ;
- }
-
- return m_SquareDetector;
+ return m_SquareDetector;
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -514,232 +520,229 @@ G4LogicalVolume* Actar::BuildDetector(){
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void Actar::ReadConfiguration(NPL::InputParser parser){
- vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Actar");
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
-
- vector<string> cart = {"POS","Shape","GasMaterial","GasFraction","Temperature","Pressure","Si","VamosSi","CsI","BeamDump"};
- //vector<string> sphe = {"R","Theta","Phi","Shape","GasMaterial","GasFraction","Temperature","Pressure","Si","CsI","BeamDump"};
-
- for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(cart)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Actar " << i+1 << endl;
-
- G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
- string Shape = blocks[i]->GetString("Shape");
- vector<string> GasName = blocks[i]->GetVectorString("GasMaterial");
- vector<int> GasFraction = blocks[i]->GetVectorInt("GasFraction");
- for(unsigned int j=0; j<GasName.size(); j++){
- m_GasMaterial.push_back(GasName[j]);
- m_GasFraction.push_back(GasFraction[j]);
- }
- m_Temperature = blocks[i]->GetDouble("Temperature","kelvin");
- m_Pressure = blocks[i]->GetDouble("Pressure","bar");
- m_build_Silicon = blocks[i]->GetInt("Si");
- m_build_Vamos_Silicon = blocks[i]->GetInt("VamosSi");
- m_build_CsI = blocks[i]->GetInt("CsI");
- m_build_BeamDump = blocks[i]->GetInt("BeamDump");
-
- AddDetector(Pos,Shape);
- }
- /*else if(blocks[i]->HasTokenList(sphe)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Actar " << 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");
- vector<string> GasName = blocks[i]->GetVectorString("GasMaterial");
- vector<int> GasFraction = blocks[i]->GetVectorInt("GasFraction");
- for(unsigned int j=0; j<GasName.size(); j++){
- m_GasMaterial.push_back(GasName[j]);
- m_GasFraction.push_back(GasFraction[j]);
- }
- m_Temperature = blocks[i]->GetDouble("Temperature","kelvin");
- m_Pressure = blocks[i]->GetDouble("Pressure","bar");
- m_build_Silicon = blocks[i]->GetInt("Si");
- m_build_CsI = blocks[i]->GetInt("CsI");
-
- AddDetector(R,Theta,Phi,Shape);
- }*/
- else{
- cout << "ERROR: check your input file formatting " << endl;
- exit(1);
+void Actar::ReadConfiguration(NPL::InputParser parser) {
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Actar");
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks.size() << " detectors found " << endl;
+
+ vector<string> cart = {"POS", "Shape", "GasMaterial", "GasFraction", "Temperature",
+ "Pressure", "Si", "VamosSi", "CsI", "BeamDump"};
+ // vector<string> sphe =
+ // {"R","Theta","Phi","Shape","GasMaterial","GasFraction","Temperature","Pressure","Si","CsI","BeamDump"};
+
+ for (unsigned int i = 0; i < blocks.size(); i++) {
+ if (blocks[i]->HasTokenList(cart)) {
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Actar " << i + 1 << endl;
+
+ G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS", "mm"));
+ string Shape = blocks[i]->GetString("Shape");
+ vector<string> GasName = blocks[i]->GetVectorString("GasMaterial");
+ vector<int> GasFraction = blocks[i]->GetVectorInt("GasFraction");
+ for (unsigned int j = 0; j < GasName.size(); j++) {
+ m_GasMaterial.push_back(GasName[j]);
+ m_GasFraction.push_back(GasFraction[j]);
+ }
+ m_Temperature = blocks[i]->GetDouble("Temperature", "kelvin");
+ m_Pressure = blocks[i]->GetDouble("Pressure", "bar");
+ m_build_Silicon = blocks[i]->GetInt("Si");
+ m_build_Vamos_Silicon = blocks[i]->GetInt("VamosSi");
+ m_build_CsI = blocks[i]->GetInt("CsI");
+ m_build_BeamDump = blocks[i]->GetInt("BeamDump");
+
+ AddDetector(Pos, Shape);
+ }
+ /*else if(blocks[i]->HasTokenList(sphe)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Actar " << 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");
+ vector<string> GasName = blocks[i]->GetVectorString("GasMaterial");
+ vector<int> GasFraction = blocks[i]->GetVectorInt("GasFraction");
+ for(unsigned int j=0; j<GasName.size(); j++){
+ m_GasMaterial.push_back(GasName[j]);
+ m_GasFraction.push_back(GasFraction[j]);
}
+ m_Temperature = blocks[i]->GetDouble("Temperature","kelvin");
+ m_Pressure = blocks[i]->GetDouble("Pressure","bar");
+ m_build_Silicon = blocks[i]->GetInt("Si");
+ m_build_CsI = blocks[i]->GetInt("CsI");
+
+ AddDetector(R,Theta,Phi,Shape);
+ }*/
+ else {
+ cout << "ERROR: check your input file formatting " << endl;
+ exit(1);
}
+ }
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
-void Actar::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()*Actar_NS::ChamberThickness*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);
-
- new G4PVPlacement(G4Transform3D(*Rot,Det_pos),
- BuildDetector(),
- "Actar",world,false,i+1);
- }
- if(!m_ReactionRegion){
- G4ProductionCuts* ecut = new G4ProductionCuts();
- ecut->SetProductionCut(1000,"e-");
-
- m_ReactionRegion= new G4Region("NPSimulationProcess");
- m_ReactionRegion->SetProductionCuts(ecut);
- m_ReactionRegion->AddRootLogicalVolume(m_logicGas);
- m_ReactionRegion->SetUserLimits(new G4UserLimits(1.2*mm));
-
- G4Region* Region_cut = new G4Region("RegionCut");
- Region_cut->SetProductionCuts(ecut);
- Region_cut->AddRootLogicalVolume(m_SquareDetector);
- }
- G4FastSimulationManager* mng = m_ReactionRegion->GetFastSimulationManager();
- unsigned int size = m_ReactionModel.size();
- for(unsigned int i = 0 ; i < size ; i++){
- mng->RemoveFastSimulationModel(m_ReactionModel[i]);
- }
- m_ReactionModel.clear();
- G4VFastSimulationModel* fsm;
- fsm = new NPS::BeamReaction("BeamReaction",m_ReactionRegion);
- m_ReactionModel.push_back(fsm);
- fsm = new NPS::Decay("Decay",m_ReactionRegion);
- m_ReactionModel.push_back(fsm);
+void Actar::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() * Actar_NS::ChamberThickness * 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);
+
+ new G4PVPlacement(G4Transform3D(*Rot, Det_pos), BuildDetector(), "Actar", world, false, i + 1);
+ }
+ if (!m_ReactionRegion) {
+ G4ProductionCuts* ecut = new G4ProductionCuts();
+ ecut->SetProductionCut(1000, "e-");
+
+ m_ReactionRegion = new G4Region("NPSimulationProcess");
+ m_ReactionRegion->SetProductionCuts(ecut);
+ m_ReactionRegion->AddRootLogicalVolume(m_logicGas);
+ m_ReactionRegion->SetUserLimits(new G4UserLimits(1.2 * mm));
+
+ G4Region* Region_cut = new G4Region("RegionCut");
+ Region_cut->SetProductionCuts(ecut);
+ Region_cut->AddRootLogicalVolume(m_SquareDetector);
+ }
+ G4FastSimulationManager* mng = m_ReactionRegion->GetFastSimulationManager();
+ unsigned int size = m_ReactionModel.size();
+ for (unsigned int i = 0; i < size; i++) {
+ mng->RemoveFastSimulationModel(m_ReactionModel[i]);
+ }
+ m_ReactionModel.clear();
+ G4VFastSimulationModel* fsm;
+ fsm = new NPS::BeamReaction("BeamReaction", m_ReactionRegion);
+ m_ReactionModel.push_back(fsm);
+ fsm = new NPS::Decay("Decay", m_ReactionRegion);
+ m_ReactionModel.push_back(fsm);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method
-void Actar::InitializeRootOutput(){
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- if(!pTree->FindBranch("data")){
- //pTree->Branch("Actar", "TActarData", &m_Event) ;
- pTree->Branch("data", "MEventReduced", &m_EventReduced) ;
- }
- //pTree->SetBranchAddress("Actar", &m_Event) ;
- pTree->SetBranchAddress("data", &m_EventReduced) ;
+void Actar::InitializeRootOutput() {
+ RootOutput* pAnalysis = RootOutput::getInstance();
+ TTree* pTree = pAnalysis->GetTree();
+ if (!pTree->FindBranch("data")) {
+ // pTree->Branch("Actar", "TActarData", &m_Event) ;
+ pTree->Branch("data", "MEventReduced", &m_EventReduced);
+ }
+ // pTree->SetBranchAddress("Actar", &m_Event) ;
+ pTree->SetBranchAddress("data", &m_EventReduced);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
-void Actar::ReadSensitive(const G4Event*){
- m_EventReduced->CoboAsad.clear();
- static ReducedData DataReduced;
-
-
- //////////////
- // Pad scorer
- TPCScorers::PS_TPCCathode* PadScorer= (TPCScorers::PS_TPCCathode*) m_ActarScorer->GetPrimitive(0);
- unsigned int size = PadScorer->GetMult();
- // Loop on the Pad map
- //TH1D* h = new TH1D("h","h",25000,0,25000);
- for (unsigned int i = 0 ; i < size ; i++){
- DataReduced.clear();
- double Count = RandGauss::shoot(PadScorer->GetCharge(i),Actar_NS::ResoCharge*PadScorer->GetCharge(i));
- double Time = PadScorer->GetTime(i);//RandGauss::shoot(Info[1],Actar_NS::ResoTime);
- //cout << "Time= " << Time << endl;
- //int iTime = ((int) Time*20/512)+1;
- //int PadNbr = Info[7];
- int Pad_X = PadScorer->GetPadX(i);//m_PadToXRow[PadNbr];
- int Pad_Y = PadScorer->GetPadY(i);//m_PadToZColumn[PadNbr];
-
- int co = m_PadToCobo[Pad_X][Pad_Y];
- int as = m_PadToAsad[Pad_X][Pad_Y];
- int ag = m_PadToAGET[Pad_X][Pad_Y];
- int ch = m_PadToChannel[Pad_X][Pad_Y];
-
- if(Count>Actar_NS::ChargeThreshold){
- DataReduced.globalchannelid = ch+(ag<<7)+(as<<9)+(co<<11);
- DataReduced.peakheight.push_back(Count);
- DataReduced.peaktime.push_back(Time);
- }
- m_EventReduced->CoboAsad.push_back(DataReduced);
+void Actar::ReadSensitive(const G4Event*) {
+ m_EventReduced->CoboAsad.clear();
+ static ReducedData DataReduced;
+
+ //////////////
+ // Pad scorer
+ TPCScorers::PS_TPCCathode* PadScorer = (TPCScorers::PS_TPCCathode*)m_ActarScorer->GetPrimitive(0);
+ unsigned int size = PadScorer->GetMult();
+ // Loop on the Pad map
+ // TH1D* h = new TH1D("h","h",25000,0,25000);
+ for (unsigned int i = 0; i < size; i++) {
+ DataReduced.clear();
+ double Count = RandGauss::shoot(PadScorer->GetCharge(i), Actar_NS::ResoCharge * PadScorer->GetCharge(i));
+ double Time = PadScorer->GetTime(i); // RandGauss::shoot(Info[1],Actar_NS::ResoTime);
+ // cout << "Time= " << Time << endl;
+ // int iTime = ((int) Time*20/512)+1;
+ // int PadNbr = Info[7];
+ int Pad_X = PadScorer->GetPadX(i); // m_PadToXRow[PadNbr];
+ int Pad_Y = PadScorer->GetPadY(i); // m_PadToZColumn[PadNbr];
+
+ int co = m_PadToCobo[Pad_X][Pad_Y];
+ int as = m_PadToAsad[Pad_X][Pad_Y];
+ int ag = m_PadToAGET[Pad_X][Pad_Y];
+ int ch = m_PadToChannel[Pad_X][Pad_Y];
+
+ if (Count > Actar_NS::ChargeThreshold) {
+ DataReduced.globalchannelid = ch + (ag << 7) + (as << 9) + (co << 11);
+ DataReduced.peakheight.push_back(Count);
+ DataReduced.peaktime.push_back(Time);
}
-
- /*vector<double> Q, T;
- for(int i=0; i<h->GetNbinsX(); i++){
- double count = h->GetBinContent(i);
- double time = h->GetBinCenter(i);
- if(count){
- Q.push_back(count);
- T.push_back(time+500);
-
- }
- }
- // clear map for next event
- SimulateDigitizer(Q,T,1.40*microsecond,0,8750,25,5);
- delete h;*/
-
- // Silicon //
- if(m_build_Silicon){
- DSSDScorers::PS_Rectangle* SiScorer= (DSSDScorers::PS_Rectangle*) m_SiliconScorer->GetPrimitive(0);
- unsigned int sizeSi = SiScorer->GetLengthMult();
- // Loop on the ThinSi map
- for(unsigned int i = 0 ; i < sizeSi ; i++){
- DataReduced.clear();
- double E_Si = RandGauss::shoot(SiScorer->GetEnergyLength(i),Actar_NS::ResoSilicon);
-
- int co = 31;
- int as = 0;
- int ag = 0;
- int ch = SiScorer->GetDetectorLength(i);
-
- if(E_Si>Actar_NS::EnergyThreshold){
- DataReduced.globalchannelid = ch+(ag<<7)+(as<<9)+(co<<11);
- DataReduced.peaktime.push_back(ch);
- DataReduced.peakheight.push_back(E_Si);
- }
- m_EventReduced->CoboAsad.push_back(DataReduced);
- }
- // Clear Map for next event
- SiScorer->clear();
+ m_EventReduced->CoboAsad.push_back(DataReduced);
+ }
+
+ /*vector<double> Q, T;
+ for(int i=0; i<h->GetNbinsX(); i++){
+ double count = h->GetBinContent(i);
+ double time = h->GetBinCenter(i);
+ if(count){
+ Q.push_back(count);
+ T.push_back(time+500);
+
+ }
+ }
+ // clear map for next event
+ SimulateDigitizer(Q,T,1.40*microsecond,0,8750,25,5);
+ delete h;*/
+
+ // Silicon //
+ if (m_build_Silicon) {
+ DSSDScorers::PS_Rectangle* SiScorer = (DSSDScorers::PS_Rectangle*)m_SiliconScorer->GetPrimitive(0);
+ unsigned int sizeSi = SiScorer->GetLengthMult();
+ // Loop on the ThinSi map
+ for (unsigned int i = 0; i < sizeSi; i++) {
+ DataReduced.clear();
+ double E_Si = RandGauss::shoot(SiScorer->GetEnergyLength(i), Actar_NS::ResoSilicon);
+
+ int co = 31;
+ int as = 0;
+ int ag = 0;
+ int ch = SiScorer->GetDetectorLength(i);
+
+ if (E_Si > Actar_NS::EnergyThreshold) {
+ DataReduced.globalchannelid = ch + (ag << 7) + (as << 9) + (co << 11);
+ DataReduced.peaktime.push_back(ch);
+ DataReduced.peakheight.push_back(E_Si);
+ }
+ m_EventReduced->CoboAsad.push_back(DataReduced);
}
-
- // CsI //
- if(m_build_CsI){
- CalorimeterScorers::PS_Calorimeter* CsIScorer= (CalorimeterScorers::PS_Calorimeter*) m_CsIScorer->GetPrimitive(0);
- unsigned int sizeCsI = CsIScorer->GetMult();
- // Loop on the ThinSi map
- for(unsigned int i = 0 ; i < sizeCsI ; i++){
- DataReduced.clear();
- vector<unsigned int> level = CsIScorer->GetLevel(i);
- //double E_CsI = RandGauss::shoot(CsIScorer->GetEnergy(i),Actar_NS::ResoCsI);
-
- /*if(E_CsI>Actar_NS::EnergyThreshold){
- m_Event->SetCsIEnergy(E_CsI);
- m_Event->SetCsICrystalNumber(level[0]);
- }
- m_EventReduced->CoboAsad.push_back(DataReduced);*/
- }
+ // Clear Map for next event
+ SiScorer->clear();
+ }
+
+ // CsI //
+ if (m_build_CsI) {
+ CalorimeterScorers::PS_Calorimeter* CsIScorer = (CalorimeterScorers::PS_Calorimeter*)m_CsIScorer->GetPrimitive(0);
+ unsigned int sizeCsI = CsIScorer->GetMult();
+ // Loop on the ThinSi map
+ for (unsigned int i = 0; i < sizeCsI; i++) {
+ DataReduced.clear();
+ vector<unsigned int> level = CsIScorer->GetLevel(i);
+ // double E_CsI = RandGauss::shoot(CsIScorer->GetEnergy(i),Actar_NS::ResoCsI);
+
+ /*if(E_CsI>Actar_NS::EnergyThreshold){
+ m_Event->SetCsIEnergy(E_CsI);
+ m_Event->SetCsICrystalNumber(level[0]);
+ }
+ m_EventReduced->CoboAsad.push_back(DataReduced);*/
}
-
-
+ }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-/*void Actar::SimulateDigitizer(vector<double> E, vector<double> T, double fallTime,double start,double stop, double step,double noise){
+/*void Actar::SimulateDigitizer(vector<double> E, vector<double> T, double fallTime,double start,double stop, double
+ step,double noise){
static string formula;
formula= "";
@@ -771,37 +774,41 @@ void Actar::ReadSensitive(const G4Event*){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////
void Actar::InitializeScorers() {
- // This check is necessary in case the geometry is reloaded
+ // This check is necessary in case the geometry is reloaded
- bool already_exist = false;
- vector<G4int> NestingLevel;
- NestingLevel.push_back(0);
- NestingLevel.push_back(2);
+ bool already_exist = false;
+ vector<G4int> NestingLevel;
+ NestingLevel.push_back(0);
+ NestingLevel.push_back(2);
- m_ActarScorer = CheckScorer("ActarScorer",already_exist) ;
- m_SiliconScorer = CheckScorer("SiliconScorer",already_exist);
- //m_VamosSiliconScorer = CheckScorer("VamosSiliconScorer",already_exist);
- m_CsIScorer = CheckScorer("CsIScorer",already_exist);
+ m_ActarScorer = CheckScorer("ActarScorer", already_exist);
+ m_SiliconScorer = CheckScorer("SiliconScorer", already_exist);
+ // m_VamosSiliconScorer = CheckScorer("VamosSiliconScorer",already_exist);
+ m_CsIScorer = CheckScorer("CsIScorer", already_exist);
- if(already_exist) return;
+ if (already_exist)
+ return;
- G4VPrimitiveScorer* SiScorer = new DSSDScorers::PS_Rectangle("SiliconScorer",0,Actar_NS::SiliconHeight,Actar_NS::SiliconWidth,1,1);
- m_SiliconScorer->RegisterPrimitive(SiScorer);
+ G4VPrimitiveScorer* SiScorer =
+ new DSSDScorers::PS_Rectangle("SiliconScorer", 0, Actar_NS::SiliconHeight, Actar_NS::SiliconWidth, 1, 1);
+ m_SiliconScorer->RegisterPrimitive(SiScorer);
- /*G4VPrimitiveScorer* VamosSiScorer = new SILICONSCORERS::PS_Silicon_Rectangle("VamosSiliconScorer",0,Actar_NS::VamosSiliconHeight,Actar_NS::VamosSiliconWidth,1,1);
- m_VamosSiliconScorer->RegisterPrimitive(VamosSiScorer);*/
+ /*G4VPrimitiveScorer* VamosSiScorer = new
+ SILICONSCORERS::PS_Silicon_Rectangle("VamosSiliconScorer",0,Actar_NS::VamosSiliconHeight,Actar_NS::VamosSiliconWidth,1,1);
+ m_VamosSiliconScorer->RegisterPrimitive(VamosSiScorer);*/
- G4VPrimitiveScorer* CsIScorer= new CalorimeterScorers::PS_Calorimeter("CsI",NestingLevel);
- m_CsIScorer->RegisterPrimitive(CsIScorer);
+ G4VPrimitiveScorer* CsIScorer = new CalorimeterScorers::PS_Calorimeter("CsI", NestingLevel);
+ m_CsIScorer->RegisterPrimitive(CsIScorer);
- vector<int> level; level.push_back(0);
- G4VPrimitiveScorer* Actar_dig= new TPCScorers::PS_TPCCathode("Actar_dig",0) ;
- m_ActarScorer->RegisterPrimitive(Actar_dig);
+ vector<int> level;
+ level.push_back(0);
+ G4VPrimitiveScorer* Actar_dig = new TPCScorers::PS_TPCCathode("Actar_dig", 0);
+ m_ActarScorer->RegisterPrimitive(Actar_dig);
- G4SDManager::GetSDMpointer()->AddNewDetector(m_ActarScorer);
- G4SDManager::GetSDMpointer()->AddNewDetector(m_SiliconScorer);
- //G4SDManager::GetSDMpointer()->AddNewDetector(m_VamosSiliconScorer);
- G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer) ;
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_ActarScorer);
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_SiliconScorer);
+ // G4SDManager::GetSDMpointer()->AddNewDetector(m_VamosSiliconScorer);
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -810,22 +817,20 @@ void Actar::InitializeScorers() {
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
-NPS::VDetector* Actar::Construct(){
- return (NPS::VDetector*) new Actar();
-}
+NPS::VDetector* Actar::Construct() { return (NPS::VDetector*)new Actar(); }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
-extern"C" {
- class proxy_nps_Actar{
- public:
- proxy_nps_Actar(){
- NPS::DetectorFactory::getInstance()->AddToken("Actar","Actar");
- NPS::DetectorFactory::getInstance()->AddDetector("Actar",Actar::Construct);
- }
- };
-
- proxy_nps_Actar p_nps_Actar;
+extern "C" {
+class proxy_nps_Actar {
+ public:
+ proxy_nps_Actar() {
+ NPS::DetectorFactory::getInstance()->AddToken("Actar", "Actar");
+ NPS::DetectorFactory::getInstance()->AddDetector("Actar", Actar::Construct);
+ }
+};
+
+proxy_nps_Actar p_nps_Actar;
}
diff --git a/NPSimulation/Detectors/Chio/Chio.cc b/NPSimulation/Detectors/Chio/Chio.cc
index 4d6f446f893258209ff32d737562bcf3e30b8575..b1d3c77db7eb6ee0795928e906544bfeb622d4a7 100644
--- a/NPSimulation/Detectors/Chio/Chio.cc
+++ b/NPSimulation/Detectors/Chio/Chio.cc
@@ -20,100 +20,95 @@
*****************************************************************************/
// C++ headers
-#include <sstream>
#include <cmath>
#include <limits>
-//G4 Geometry object
-#include "G4Tubs.hh"
+#include <sstream>
+// G4 Geometry object
#include "G4Box.hh"
+#include "G4Tubs.hh"
-//G4 sensitive
-#include "G4SDManager.hh"
+// G4 sensitive
#include "G4MultiFunctionalDetector.hh"
+#include "G4SDManager.hh"
-//G4 various object
+// G4 various object
+#include "G4Colour.hh"
#include "G4Material.hh"
-#include "G4Transform3D.hh"
#include "G4PVPlacement.hh"
+#include "G4Transform3D.hh"
#include "G4VisAttributes.hh"
-#include "G4Colour.hh"
// G4 Field
-#include "G4FieldManager.hh"
+#include "G4ChordFinder.hh"
+#include "G4ClassicalRK4.hh"
#include "G4ElectricField.hh"
-#include "G4UniformElectricField.hh"
-#include "G4TransportationManager.hh"
#include "G4EqMagElectricField.hh"
-#include "G4MagIntegratorStepper.hh"
-#include "G4ClassicalRK4.hh"
+#include "G4FieldManager.hh"
#include "G4MagIntegratorDriver.hh"
-#include "G4ChordFinder.hh"
+#include "G4MagIntegratorStepper.hh"
#include "G4MaterialPropertiesTable.hh"
+#include "G4TransportationManager.hh"
+#include "G4UniformElectricField.hh"
// NPTool header
#include "Chio.hh"
#include "DriftElectronScorers.hh"
-#include "RootOutput.h"
+#include "FastDriftElectron.hh"
#include "MaterialManager.hh"
-#include "NPSDetectorFactory.hh"
#include "NPOptionManager.h"
-#include "FastDriftElectron.hh"
+#include "NPSDetectorFactory.hh"
+#include "RootOutput.h"
// CLHEP header
#include "CLHEP/Random/RandGauss.h"
// ROOT
-#include "TH1D.h"
#include "TF1.h"
+#include "TH1D.h"
using namespace std;
using namespace CLHEP;
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-namespace Chio_NS{
+namespace Chio_NS {
// Energy and time Resolution
const double ChargeThreshold = 1;
- const double ResoTime = 4.5*ns ;
+ const double ResoTime = 4.5 * ns;
// const double ResoEnergy = 1.0*MeV ;
- //const double Radius = 50*mm ;
+ // const double Radius = 50*mm ;
// const double Width = 100*mm ;
- const double Thickness = 300*mm ;
+ const double Thickness = 300 * mm;
const string Material = "BC400";
-}
+} // namespace Chio_NS
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Chio Specific Method
-Chio::Chio(){
- m_Event_an = new TChio_anData() ;
- m_Event_dig = new TChio_digData() ;
+Chio::Chio() {
+ m_Event_an = new TChio_anData();
+ m_Event_dig = new TChio_digData();
m_ChioScorer = 0;
m_SquareDetector = 0;
m_CylindricalDetector = 0;
-
// RGB Color + Transparency
- m_VisChamber = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.25));
- m_VisWindows = new G4VisAttributes(G4Colour(1, 0, 0, 0.25));
- m_VisGas = new G4VisAttributes(G4Colour(0, 0.5, 0.5, 0.5));
-
+ m_VisChamber = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.25));
+ m_VisWindows = new G4VisAttributes(G4Colour(1, 0, 0, 0.25));
+ m_VisGas = new G4VisAttributes(G4Colour(0, 0.5, 0.5, 0.5));
}
-Chio::~Chio(){
-}
+Chio::~Chio() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Chio::AddDetector(G4ThreeVector POS, string Shape){
- // Convert the POS value to R theta Phi as Spherical coordinate is easier in G4
+void Chio::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 Chio::AddDetector(double R, double Theta, double Phi, string Shape){
+void Chio::AddDetector(double R, double Theta, double Phi, string Shape) {
m_R.push_back(R);
m_Theta.push_back(Theta);
m_Phi.push_back(Phi);
@@ -121,102 +116,86 @@ void Chio::AddDetector(double R, double Theta, double Phi, string Shape){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* Chio::BuildDetector(){
- if(!m_SquareDetector){
+G4LogicalVolume* Chio::BuildDetector() {
+ if (!m_SquareDetector) {
// Main volume
- G4Box* sChamber = new G4Box("Chio_Box",7*cm*0.5,
- 7*cm*0.5,12*cm*0.5+2*12*micrometer*0.5);
+ G4Box* sChamber = new G4Box("Chio_Box", 7 * cm * 0.5, 7 * cm * 0.5, 12 * cm * 0.5 + 2 * 12 * micrometer * 0.5);
// Gaz volume
- G4Box* sGas = new G4Box("Chio_Gas",6*cm*0.5,
- 6*cm*0.5,12*cm*0.5-2*12*micrometer*0.5);
+ G4Box* sGas = new G4Box("Chio_Gas", 6 * cm * 0.5, 6 * cm * 0.5, 12 * cm * 0.5 - 2 * 12 * micrometer * 0.5);
// Frish grid
- G4Box* sGrid = new G4Box("Chio_Grid",1*um*0.5,
- 6*cm*0.5,12*cm*0.5-2*12*micrometer*0.5);
+ G4Box* sGrid = new G4Box("Chio_Grid", 1 * um * 0.5, 6 * cm * 0.5, 12 * cm * 0.5 - 2 * 12 * micrometer * 0.5);
// Cathode
- G4Box* sCathode = new G4Box("Chio_Cathode",1*um*0.5,
- 6*cm*0.5,12*cm*0.5-2*12*micrometer*0.5);
-
+ G4Box* sCathode = new G4Box("Chio_Cathode", 1 * um * 0.5, 6 * cm * 0.5, 12 * cm * 0.5 - 2 * 12 * micrometer * 0.5);
// Entrance/Exit windows
- G4Box* sWindows = new G4Box("Chio_Windows",7*cm*0.5,
- 7*cm*0.5,12*micrometer*0.5);
+ G4Box* sWindows = new G4Box("Chio_Windows", 7 * cm * 0.5, 7 * cm * 0.5, 12 * micrometer * 0.5);
+ G4Material* Fe = MaterialManager::getInstance()->GetMaterialFromLibrary("Fe");
+ G4Material* Al = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
- G4Material* Fe= MaterialManager::getInstance()->GetMaterialFromLibrary("Fe");
- G4Material* Al= MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ // G4Material* CF4= MaterialManager::getInstance()->GetGasFromLibrary("CF4",0.0693276*bar,273.15*kelvin);
+ G4Material* CF4 = MaterialManager::getInstance()->GetGasFromLibrary("iC4H10", 8. * bar / 1000., 273.15 * kelvin);
- //G4Material* CF4= MaterialManager::getInstance()->GetGasFromLibrary("CF4",0.0693276*bar,273.15*kelvin);
- G4Material* CF4= MaterialManager::getInstance()->GetGasFromLibrary("iC4H10",8.*bar/1000.,273.15*kelvin);
+ G4Material* Mylar = MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
- G4Material* Mylar= MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
-
- G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
- MPT->AddConstProperty("DE_PAIRENERGY",30*eV);
- // MPT->AddConstProperty("DE_AMPLIFICATION",1e4);
- MPT->AddConstProperty("DE_ABSLENGTH",1*pc);
- MPT->AddConstProperty("DE_DRIFTSPEED",11*cm/microsecond);
- MPT->AddConstProperty("DE_TRANSVERSALSPREAD",6e-5*mm2/ns);
- MPT->AddConstProperty("DE_LONGITUDINALSPREAD",4e-5*mm2/ns);
+ G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
+ MPT->AddConstProperty("DE_PAIRENERGY", 30 * eV, true);
+ // MPT->AddConstProperty("DE_AMPLIFICATION",1e4,true);
+ MPT->AddConstProperty("DE_ABSLENGTH", 1 * pc, true);
+ MPT->AddConstProperty("DE_DRIFTSPEED", 11 * cm / microsecond, true);
+ MPT->AddConstProperty("DE_TRANSVERSALSPREAD", 6e-5 * mm2 / ns, true);
+ MPT->AddConstProperty("DE_LONGITUDINALSPREAD", 4e-5 * mm2 / ns, true);
CF4->SetMaterialPropertiesTable(MPT);
- G4MaterialPropertiesTable* MPT2 = new G4MaterialPropertiesTable();
- MPT2->AddConstProperty("DE_YIELD",1);
- MPT2->AddConstProperty("DE_AMPLIFICATION",2);
- MPT2->AddConstProperty("DE_ABSLENGTH",1*pc);
+ G4MaterialPropertiesTable* MPT2 = new G4MaterialPropertiesTable();
+ MPT2->AddConstProperty("DE_YIELD", 1, true);
+ MPT2->AddConstProperty("DE_AMPLIFICATION", 2, true);
+ MPT2->AddConstProperty("DE_ABSLENGTH", 1 * pc, true);
Al->SetMaterialPropertiesTable(MPT2);
- m_SquareDetector = new G4LogicalVolume(sChamber,Fe,"logic_Chio_Box",0,0,0);
- G4LogicalVolume* logicGas = new G4LogicalVolume(sGas,CF4,"logic_Gas",0,0,0);
- G4LogicalVolume* logicGrid = new G4LogicalVolume(sGrid,Al,"logic_Grid",0,0,0);
- G4LogicalVolume* logicCathode = new G4LogicalVolume(sCathode,Fe,"logic_Cathode",0,0,0);
- G4LogicalVolume* logicWindows = new G4LogicalVolume(sWindows,Mylar,"logic_Windows",0,0,0);
+ m_SquareDetector = new G4LogicalVolume(sChamber, Fe, "logic_Chio_Box", 0, 0, 0);
+ G4LogicalVolume* logicGas = new G4LogicalVolume(sGas, CF4, "logic_Gas", 0, 0, 0);
+ G4LogicalVolume* logicGrid = new G4LogicalVolume(sGrid, Al, "logic_Grid", 0, 0, 0);
+ G4LogicalVolume* logicCathode = new G4LogicalVolume(sCathode, Fe, "logic_Cathode", 0, 0, 0);
+ G4LogicalVolume* logicWindows = new G4LogicalVolume(sWindows, Mylar, "logic_Windows", 0, 0, 0);
G4RotationMatrix* Rot = new G4RotationMatrix();
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,0)),
- logicGas,
- "ChioGas",m_SquareDetector,false,0);
+ new G4PVPlacement(G4Transform3D(*Rot, G4ThreeVector(0, 0, 0)), logicGas, "ChioGas", m_SquareDetector, false, 0);
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(2.5*cm,0,0)),
- logicGrid,
- "ChioGrid",logicGas,false,0);
+ new G4PVPlacement(G4Transform3D(*Rot, G4ThreeVector(2.5 * cm, 0, 0)), logicGrid, "ChioGrid", logicGas, false, 0);
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(3*cm-0.5*1*um,0,0)),
- logicCathode,
- "ChioCathode",logicGas,false,0);
+ new G4PVPlacement(G4Transform3D(*Rot, G4ThreeVector(3 * cm - 0.5 * 1 * um, 0, 0)), logicCathode, "ChioCathode",
+ logicGas, false, 0);
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,6*cm-6*micrometer)),
- logicWindows,
- "ChioExitWindows",m_SquareDetector,false,0);
+ new G4PVPlacement(G4Transform3D(*Rot, G4ThreeVector(0, 0, 6 * cm - 6 * micrometer)), logicWindows,
+ "ChioExitWindows", m_SquareDetector, false, 0);
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-6*cm+6*micrometer)),
- logicWindows,
- "ChioEntranceWindows",m_SquareDetector,false,0);
+ new G4PVPlacement(G4Transform3D(*Rot, G4ThreeVector(0, 0, -6 * cm + 6 * micrometer)), logicWindows,
+ "ChioEntranceWindows", m_SquareDetector, false, 0);
- G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0,-1000*volt/cm,0.0));
+ G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0, -1000 * volt / cm, 0.0));
// Create an equation of motion for this field
- G4EqMagElectricField* Equation = new G4EqMagElectricField(field);
- G4MagIntegratorStepper* Stepper = new G4ClassicalRK4( Equation, 8 );
+ G4EqMagElectricField* Equation = new G4EqMagElectricField(field);
+ G4MagIntegratorStepper* Stepper = new G4ClassicalRK4(Equation, 8);
- // Get the global field manager
- G4FieldManager* FieldManager= new G4FieldManager();
- // Set this field to the global field manager
- FieldManager->SetDetectorField(field );
- logicGas->SetFieldManager(FieldManager,true);
+ // Get the global field manager
+ G4FieldManager* FieldManager = new G4FieldManager();
+ // Set this field to the global field manager
+ FieldManager->SetDetectorField(field);
+ logicGas->SetFieldManager(FieldManager, true);
- G4MagInt_Driver* IntgrDriver = new G4MagInt_Driver(0.1*mm,
- Stepper,
- Stepper->GetNumberOfVariables() );
+ G4MagInt_Driver* IntgrDriver = new G4MagInt_Driver(0.1 * mm, Stepper, Stepper->GetNumberOfVariables());
G4ChordFinder* ChordFinder = new G4ChordFinder(IntgrDriver);
- FieldManager->SetChordFinder( ChordFinder );
+ FieldManager->SetChordFinder(ChordFinder);
- logicCathode->SetSensitiveDetector(m_ChioScorer);
+ logicCathode->SetSensitiveDetector(m_ChioScorer);
m_SquareDetector->SetVisAttributes(m_VisChamber);
logicGas->SetVisAttributes(m_VisGas);
logicWindows->SetVisAttributes(m_VisWindows);
@@ -232,96 +211,92 @@ G4LogicalVolume* Chio::BuildDetector(){
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void Chio::ReadConfiguration(NPL::InputParser parser){
+void Chio::ReadConfiguration(NPL::InputParser parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Chio");
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks.size() << " detectors found " << endl;
- vector<string> cart = {"POS","Shape"};
- vector<string> sphe = {"R","Theta","Phi","Shape"};
+ vector<string> cart = {"POS", "Shape"};
+ vector<string> sphe = {"R", "Theta", "Phi", "Shape"};
- for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(cart)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Chio " << i+1 << endl;
+ for (unsigned int i = 0; i < blocks.size(); i++) {
+ if (blocks[i]->HasTokenList(cart)) {
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Chio " << i + 1 << endl;
- G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
+ G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS", "mm"));
string Shape = blocks[i]->GetString("Shape");
- AddDetector(Pos,Shape);
+ AddDetector(Pos, Shape);
}
- else if(blocks[i]->HasTokenList(sphe)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Chio " << i+1 << endl;
- double R = blocks[i]->GetDouble("R","mm");
- double Theta = blocks[i]->GetDouble("Theta","deg");
- double Phi = blocks[i]->GetDouble("Phi","deg");
+ else if (blocks[i]->HasTokenList(sphe)) {
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Chio " << 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);
+ AddDetector(R, Theta, Phi, Shape);
}
- else{
+ else {
cout << "ERROR: check your input file formatting " << endl;
exit(1);
}
}
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
-void Chio::ConstructDetector(G4LogicalVolume* world){
- for (unsigned short i = 0 ; i < m_R.size() ; i++) {
+void Chio::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) ;
+ 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()*Chio_NS::Thickness*0.5;
+ Det_pos += Det_pos.unit() * Chio_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 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);
+ G4RotationMatrix* Rot = new G4RotationMatrix(u, v, w);
- if(m_Shape[i] == "Square"){
- new G4PVPlacement(G4Transform3D(*Rot,Det_pos),
- BuildDetector(),
- "Chio",world,false,i+1);
+ if (m_Shape[i] == "Square") {
+ new G4PVPlacement(G4Transform3D(*Rot, Det_pos), BuildDetector(), "Chio", world, false, i + 1);
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method
-void Chio::InitializeRootOutput(){
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- if(!pTree->FindBranch("ChioAn")){
- pTree->Branch("ChioAn", "TChio_anData", &m_Event_an) ;
+void Chio::InitializeRootOutput() {
+ RootOutput* pAnalysis = RootOutput::getInstance();
+ TTree* pTree = pAnalysis->GetTree();
+ if (!pTree->FindBranch("ChioAn")) {
+ pTree->Branch("ChioAn", "TChio_anData", &m_Event_an);
}
- pTree->SetBranchAddress("ChioAn", &m_Event_an) ;
+ pTree->SetBranchAddress("ChioAn", &m_Event_an);
- ///////////////
- if(!pTree->FindBranch("ChioDig")){
- pTree->Branch("ChioDig", "TChio_digData", &m_Event_dig) ;
+ ///////////////
+ if (!pTree->FindBranch("ChioDig")) {
+ pTree->Branch("ChioDig", "TChio_digData", &m_Event_dig);
}
- pTree->SetBranchAddress("ChioDig", &m_Event_dig) ;
-
+ pTree->SetBranchAddress("ChioDig", &m_Event_dig);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
-void Chio::ReadSensitive(const G4Event* event){
+void Chio::ReadSensitive(const G4Event* event) {
m_Event_an->Clear();
///////////
@@ -333,12 +308,12 @@ void Chio::ReadSensitive(const G4Event* event){
CathodeHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CathodeCollectionID));
// Loop on the Cathode map
- for (Cathode_itr = CathodeHitMap->GetMap()->begin() ; Cathode_itr != CathodeHitMap->GetMap()->end() ; Cathode_itr++){
+ for (Cathode_itr = CathodeHitMap->GetMap()->begin(); Cathode_itr != CathodeHitMap->GetMap()->end(); Cathode_itr++) {
G4double* Info = *(Cathode_itr->second);
- double Count= Info[0];
- if(Count>Chio_NS::ChargeThreshold-1){
- double Time = RandGauss::shoot(Info[1],Chio_NS::ResoTime);
- m_Event_an->SetEnergyAndTime(Count,Time,Info[2]);
+ double Count = Info[0];
+ if (Count > Chio_NS::ChargeThreshold - 1) {
+ double Time = RandGauss::shoot(Info[1], Chio_NS::ResoTime);
+ m_Event_an->SetEnergyAndTime(Count, Time, Info[2]);
}
}
// clear map for next event
@@ -351,81 +326,81 @@ void Chio::ReadSensitive(const G4Event* event){
CathodeHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CathodeCollectionID));
// Loop on the Cathode map
- TH1D* h = new TH1D("h","h",25000,0,25000);
- for (Cathode_itr = CathodeHitMap->GetMap()->begin() ; Cathode_itr != CathodeHitMap->GetMap()->end() ; Cathode_itr++){
+ TH1D* h = new TH1D("h", "h", 25000, 0, 25000);
+ for (Cathode_itr = CathodeHitMap->GetMap()->begin(); Cathode_itr != CathodeHitMap->GetMap()->end(); Cathode_itr++) {
G4double* Info = *(Cathode_itr->second);
- if(Info[0]){
- h->Fill(Info[1],Info[0]);
+ if (Info[0]) {
+ h->Fill(Info[1], Info[0]);
}
}
- vector<double> E,T;
- for(int i = 0 ; i < h->GetNbinsX() ; i++){
+ vector<double> E, T;
+ for (int i = 0; i < h->GetNbinsX(); i++) {
double count = h->GetBinContent(i);
- double time = h->GetBinCenter(i);
- if(count)
- // m_Event_dig->AddEnergyPoint(count,time);
- E.push_back(count);
- T.push_back(time+500);
+ double time = h->GetBinCenter(i);
+ if (count)
+ // m_Event_dig->AddEnergyPoint(count,time);
+ E.push_back(count);
+ T.push_back(time + 500);
}
- SimulateDigitizer(E,T,1.40*microsecond,0,8750,25,5);
+ SimulateDigitizer(E, T, 1.40 * microsecond, 0, 8750, 25, 5);
delete h;
// clear map for next event
CathodeHitMap->clear();
-
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Chio::SimulateDigitizer(vector<double> E, vector<double> T, double fallTime,double start,double stop, double step,double noise){
+void Chio::SimulateDigitizer(vector<double> E, vector<double> T, double fallTime, double start, double stop,
+ double step, double noise) {
- static string formula;
- formula= "";
- static string Es,Ts,var,cond;
+ static string formula;
+ formula = "";
+ static string Es, Ts, var, cond;
static string fall;
- fall=std::to_string(fallTime);
+ fall = std::to_string(fallTime);
- for(unsigned int i = 0 ; i < E.size() ; i++){
- if(E[i]!=0 && T[i]!=0){
+ for (unsigned int i = 0; i < E.size(); i++) {
+ if (E[i] != 0 && T[i] != 0) {
Es = std::to_string(E[i]);
Ts = std::to_string(T[i]);
- cond = ")*(x>"+Ts+")+";
- var = "(x-"+Ts+")";
- formula += Es+"*-1*exp(-"+var+"/"+fall+cond;
+ cond = ")*(x>" + Ts + ")+";
+ var = "(x-" + Ts + ")";
+ formula += Es + "*-1*exp(-" + var + "/" + fall + cond;
}
}
- formula+="0";
- TF1* f = new TF1("f",formula.c_str(),start,stop);
- unsigned int size = (stop-start)/step;
- for(unsigned int i = 0 ; i < size ; i++){
- double time = start+i*step;
- double energy = f->Eval(time)+noise*(1-2*G4UniformRand());
- m_Event_dig->AddEnergyPoint(energy,time);
+ formula += "0";
+ TF1* f = new TF1("f", formula.c_str(), start, stop);
+ unsigned int size = (stop - start) / step;
+ for (unsigned int i = 0; i < size; i++) {
+ double time = start + i * step;
+ double energy = f->Eval(time) + noise * (1 - 2 * G4UniformRand());
+ m_Event_dig->AddEnergyPoint(energy, time);
}
-
- delete f;
+
+ delete f;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-////////////////////////////////////////////////////////////////
-void Chio::InitializeScorers() {
+////////////////////////////////////////////////////////////////
+void Chio::InitializeScorers() {
// This check is necessary in case the geometry is reloaded
- bool already_exist = false;
- m_ChioScorer = CheckScorer("ChioScorer",already_exist) ;
+ bool already_exist = false;
+ m_ChioScorer = CheckScorer("ChioScorer", already_exist);
- if(already_exist)
- return ;
+ if (already_exist)
+ return;
// Otherwise the scorer is initialised
- G4VPrimitiveScorer* Cathode_an= new DRIFTELECTRONSCORERS::PS_DECathode("Cathode_an",0) ;
- G4VPrimitiveScorer* Cathode_dig= new DRIFTELECTRONSCORERS::PS_DEDigitizer("Cathode_dig",0) ;
+ G4VPrimitiveScorer* Cathode_an = new DRIFTELECTRONSCORERS::PS_DECathode("Cathode_an", 0);
+ G4VPrimitiveScorer* Cathode_dig = new DRIFTELECTRONSCORERS::PS_DEDigitizer("Cathode_dig", 0);
- //and register it to the multifunctionnal detector
+ // and register it to the multifunctionnal detector
m_ChioScorer->RegisterPrimitive(Cathode_an);
m_ChioScorer->RegisterPrimitive(Cathode_dig);
- G4SDManager::GetSDMpointer()->AddNewDetector(m_ChioScorer) ;
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_ChioScorer);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -434,22 +409,20 @@ void Chio::InitializeScorers() {
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
-NPS::VDetector* Chio::Construct(){
- return (NPS::VDetector*) new Chio();
-}
+NPS::VDetector* Chio::Construct() { return (NPS::VDetector*)new Chio(); }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
-extern"C" {
- class proxy_nps_Chio{
- public:
- proxy_nps_Chio(){
- NPS::DetectorFactory::getInstance()->AddToken("Chio","Chio");
- NPS::DetectorFactory::getInstance()->AddDetector("Chio",Chio::Construct);
- }
- };
-
- proxy_nps_Chio p_nps_Chio;
+extern "C" {
+class proxy_nps_Chio {
+ public:
+ proxy_nps_Chio() {
+ NPS::DetectorFactory::getInstance()->AddToken("Chio", "Chio");
+ NPS::DetectorFactory::getInstance()->AddDetector("Chio", Chio::Construct);
+ }
+};
+
+proxy_nps_Chio p_nps_Chio;
}
diff --git a/NPSimulation/Detectors/SEASON/SEASON.cc b/NPSimulation/Detectors/SEASON/SEASON.cc
index 4e785717dabf624b3aeb9beb391e4b28adecda7b..4fc96a0e1cd86757834a8d30d5288a705e2ee983 100644
--- a/NPSimulation/Detectors/SEASON/SEASON.cc
+++ b/NPSimulation/Detectors/SEASON/SEASON.cc
@@ -1,59 +1,59 @@
/*****************************************************************************
-* Copyright (C) 2009-2023 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 *
-*****************************************************************************/
+ * Copyright (C) 2009-2023 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: Emmanuel Rey-herme *
-* contact address: marine.vandebrouck@cea.fr *
-* *
-* Creation Date : septembre 2023 *
-* Last update : *
-*---------------------------------------------------------------------------*
-* Decription: *
-* This class describe SEASON simulation *
-* *
-*---------------------------------------------------------------------------*
-* Comment: *
-* *
-*****************************************************************************/
+ * Original Author: Emmanuel Rey-herme *
+ * contact address: marine.vandebrouck@cea.fr *
+ * *
+ * Creation Date : septembre 2023 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class describe SEASON simulation *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
// C++ headers
-#include <sstream>
#include <cmath>
#include <limits>
-//G4 Geometry object
-#include "G4Tubs.hh"
+#include <sstream>
+// G4 Geometry object
#include "G4Box.hh"
-#include "G4Trd.hh"
+#include "G4MultiUnion.hh"
#include "G4SubtractionSolid.hh"
+#include "G4Trd.hh"
+#include "G4Tubs.hh"
#include "G4UnionSolid.hh"
-#include "G4MultiUnion.hh"
#include "G4VSolid.hh"
-//G4 sensitive
-#include "G4SDManager.hh"
+// G4 sensitive
#include "G4MultiFunctionalDetector.hh"
+#include "G4SDManager.hh"
-//G4 various object
+// G4 various object
+#include "G4Colour.hh"
#include "G4Material.hh"
-#include "G4Transform3D.hh"
#include "G4PVPlacement.hh"
+#include "G4Transform3D.hh"
#include "G4VisAttributes.hh"
-#include "G4Colour.hh"
#include "Randomize.hh"
// NPTool header
-#include "SEASON.hh"
#include "DSSDScorers.hh"
#include "InteractionScorers.hh"
-#include "RootOutput.h"
#include "MaterialManager.hh"
-#include "NPSDetectorFactory.hh"
#include "NPOptionManager.h"
+#include "NPSDetectorFactory.hh"
#include "NPSHitsMap.hh"
+#include "RootOutput.h"
+#include "SEASON.hh"
// CLHEP header
#include "CLHEP/Random/RandGauss.h"
@@ -62,67 +62,65 @@
using namespace std;
using namespace CLHEP;
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-namespace SEASON_NS{
+namespace SEASON_NS {
// Energy and time Resolution
- const double ResoTime = 8.5*ns;
- const double ResoAlpha = 6.37*keV; // Equivalent to 15 keV FWHM
- const double ResoElectron = 2.97*keV; // Equivalent to 7 keV FWHM;
-
+ const double ResoTime = 8.5 * ns;
+ const double ResoAlpha = 6.37 * keV; // Equivalent to 15 keV FWHM
+ const double ResoElectron = 2.97 * keV; // Equivalent to 7 keV FWHM;
+
//
- const double Thickness = 1*mm; //Thickness of DSSDs
- const double FoilWheelCenterDist = 16.*cm; // distance from foil center to wheel center
- const double WheelRadius = 11.5*cm;
- const double WheelThickness = 2*mm;
- const double FinThickness = 0.5*mm;
- const double FinLenght = 6.*cm;
- const double FinWidth = 3.*cm;
+ const double Thickness = 1 * mm; // Thickness of DSSDs
+ const double FoilWheelCenterDist = 16. * cm; // distance from foil center to wheel center
+ const double WheelRadius = 11.5 * cm;
+ const double WheelThickness = 2 * mm;
+ const double FinThickness = 0.5 * mm;
+ const double FinLenght = 6. * cm;
+ const double FinWidth = 3. * cm;
const int FoilNbr = 11;
- const double DistDSSD1_AlWindow = 9*mm;
-
- const double StripWidth = 2*mm;
- const double StripPitch = 1.925*mm;
- const double StripLength = 63.96*mm;
- const double DetSize = 67.975*mm;
+ const double DistDSSD1_AlWindow = 9 * mm;
+
+ const double StripWidth = 2 * mm;
+ const double StripPitch = 1.925 * mm;
+ const double StripLength = 63.96 * mm;
+ const double DetSize = 67.975 * mm;
const int NumberOfStripsX = 32;
const int NumberOfStripsY = 32;
-
+
const double DeadLayerThickness = 50 * nm;
const double GridThickness = 300 * nm;
const double GridWidth = 30 * um;
-
- G4Material * FoilMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("C");
+
+ G4Material* FoilMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("C");
G4Material* DeadLayerMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
G4Material* Grid_Material = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
-}
+} // namespace SEASON_NS
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// SEASON Specific Method
-SEASON::SEASON(){
- m_Event = new TSEASONData() ;
+SEASON::SEASON() {
+ m_Event = new TSEASONData();
m_SEASONScorer = 0;
-
- m_DistDSSD1 = 2*mm;
- m_DistTunnel = 2*mm;
-
+
+ m_DistDSSD1 = 2 * mm;
+ m_DistTunnel = 2 * mm;
+
// RGB Color + Transparency
m_VisSquare = new G4VisAttributes(G4Colour(0, 1, 0.5, 1));
m_VisWheel = new G4VisAttributes(G4Colour(1, 1, 0, 1));
m_VisFoil = new G4VisAttributes(G4Colour(1, 0, 0, 1));
m_VisGe = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.5));
- m_VisWindow = new G4VisAttributes(G4VisAttributes::Invisible);
- m_VisGrid = new G4VisAttributes(G4VisAttributes::Invisible);
- m_VisChamber = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.5));
+ m_VisWindow = new G4VisAttributes(G4VisAttributes::GetInvisible());
+ m_VisGrid = new G4VisAttributes(G4VisAttributes::GetInvisible());
+ m_VisChamber = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.5));
}
-SEASON::~SEASON(){
-}
+SEASON::~SEASON() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void SEASON::AddDetector(G4ThreeVector X1_Y1, G4ThreeVector X1_YMax, G4ThreeVector XMax_Y1, G4ThreeVector XMax_YMax){
+void SEASON::AddDetector(G4ThreeVector X1_Y1, G4ThreeVector X1_YMax, G4ThreeVector XMax_Y1, G4ThreeVector XMax_YMax) {
m_X1_Y1.push_back(X1_Y1);
m_X1_YMax.push_back(X1_YMax);
m_XMax_Y1.push_back(XMax_Y1);
@@ -131,95 +129,117 @@ void SEASON::AddDetector(G4ThreeVector X1_Y1, G4ThreeVector X1_YMax, G4ThreeVect
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void SEASON::AddDetector(double DistDSSD1, double DistTunnel, bool Deloc, bool EnableWheel, bool RotateWheel, bool EnableChamber){
-
+void SEASON::AddDetector(double DistDSSD1, double DistTunnel, bool Deloc, bool EnableWheel, bool RotateWheel,
+ bool EnableChamber) {
+
m_EnableWheel = EnableWheel;
m_RotateWheel = RotateWheel;
m_EnableChamber = EnableChamber;
-
+
m_DistDSSD1 = DistDSSD1;
m_DistTunnel = DistTunnel;
-
- double pos = SEASON_NS::DetSize*0.5;
- double interX = 0.2*mm;
- double interY = 0.2*mm;
-
+
+ double pos = SEASON_NS::DetSize * 0.5;
+ double interX = 0.2 * mm;
+ double interY = 0.2 * mm;
+
// Add DSSD1
- AddDetector(G4ThreeVector(-pos,-pos,DistDSSD1), G4ThreeVector(-pos,pos,DistDSSD1), G4ThreeVector(pos,-pos,DistDSSD1), G4ThreeVector(pos,pos,DistDSSD1));
-
+ AddDetector(G4ThreeVector(-pos, -pos, DistDSSD1), G4ThreeVector(-pos, pos, DistDSSD1),
+ G4ThreeVector(pos, -pos, DistDSSD1), G4ThreeVector(pos, pos, DistDSSD1));
+
// Add Tunnel detectors
- AddDetector(G4ThreeVector(pos+interX,-pos,-DistTunnel), G4ThreeVector(pos+interX,pos,-DistTunnel), G4ThreeVector(pos+interX,-pos,-DistTunnel-SEASON_NS::DetSize), G4ThreeVector(pos+interX,pos,-DistTunnel-SEASON_NS::DetSize));
- AddDetector(G4ThreeVector(pos,pos+interY,-DistTunnel), G4ThreeVector(-pos,pos+interY,-DistTunnel), G4ThreeVector(pos,pos+interY,-DistTunnel-SEASON_NS::DetSize), G4ThreeVector(-pos,pos+interY,-DistTunnel-SEASON_NS::DetSize));
- AddDetector(G4ThreeVector(-pos-interX,pos,-DistTunnel), G4ThreeVector(-pos-interX,-pos,-DistTunnel), G4ThreeVector(-pos-interX,pos,-DistTunnel-SEASON_NS::DetSize), G4ThreeVector(-pos,-pos-interX,-DistTunnel-SEASON_NS::DetSize));
- AddDetector(G4ThreeVector(-pos,-pos-interY,-DistTunnel), G4ThreeVector(pos,-pos-interY,-DistTunnel), G4ThreeVector(-pos,-pos-interY,-DistTunnel-SEASON_NS::DetSize), G4ThreeVector(pos,-pos-interY,-DistTunnel-SEASON_NS::DetSize));
-
+ AddDetector(G4ThreeVector(pos + interX, -pos, -DistTunnel), G4ThreeVector(pos + interX, pos, -DistTunnel),
+ G4ThreeVector(pos + interX, -pos, -DistTunnel - SEASON_NS::DetSize),
+ G4ThreeVector(pos + interX, pos, -DistTunnel - SEASON_NS::DetSize));
+ AddDetector(G4ThreeVector(pos, pos + interY, -DistTunnel), G4ThreeVector(-pos, pos + interY, -DistTunnel),
+ G4ThreeVector(pos, pos + interY, -DistTunnel - SEASON_NS::DetSize),
+ G4ThreeVector(-pos, pos + interY, -DistTunnel - SEASON_NS::DetSize));
+ AddDetector(G4ThreeVector(-pos - interX, pos, -DistTunnel), G4ThreeVector(-pos - interX, -pos, -DistTunnel),
+ G4ThreeVector(-pos - interX, pos, -DistTunnel - SEASON_NS::DetSize),
+ G4ThreeVector(-pos, -pos - interX, -DistTunnel - SEASON_NS::DetSize));
+ AddDetector(G4ThreeVector(-pos, -pos - interY, -DistTunnel), G4ThreeVector(pos, -pos - interY, -DistTunnel),
+ G4ThreeVector(-pos, -pos - interY, -DistTunnel - SEASON_NS::DetSize),
+ G4ThreeVector(pos, -pos - interY, -DistTunnel - SEASON_NS::DetSize));
+
// Add Delocalized station detectors if activated
- if(Deloc){
- G4double angle = -4.*2.*TMath::Pi()/SEASON_NS::FoilNbr - TMath::Pi()*0.5 ;
- G4ThreeVector * DelocFoilPos = new G4ThreeVector(TMath::Cos(angle)*SEASON_NS::FoilWheelCenterDist, (TMath::Sin(angle)+1)*SEASON_NS::FoilWheelCenterDist, 0*cm);
-
+ if (Deloc) {
+ G4double angle = -4. * 2. * TMath::Pi() / SEASON_NS::FoilNbr - TMath::Pi() * 0.5;
+ G4ThreeVector* DelocFoilPos = new G4ThreeVector(TMath::Cos(angle) * SEASON_NS::FoilWheelCenterDist,
+ (TMath::Sin(angle) + 1) * SEASON_NS::FoilWheelCenterDist, 0 * cm);
+
double DelocFoilPosX = DelocFoilPos->getX();
double DelocFoilPosY = DelocFoilPos->getY();
-
- AddDetector(G4ThreeVector(DelocFoilPosX-pos,DelocFoilPosY-pos,DistDSSD1), G4ThreeVector(DelocFoilPosX-pos,DelocFoilPosY+pos,DistDSSD1), G4ThreeVector(DelocFoilPosX+pos,DelocFoilPosY-pos,DistDSSD1), G4ThreeVector(DelocFoilPosX+pos,DelocFoilPosY+pos,DistDSSD1));
- AddDetector(G4ThreeVector(DelocFoilPosX-pos,DelocFoilPosY-pos,-DistTunnel), G4ThreeVector(DelocFoilPosX-pos,DelocFoilPosY+pos,-DistTunnel), G4ThreeVector(DelocFoilPosX+pos,DelocFoilPosY-pos,-DistTunnel), G4ThreeVector(DelocFoilPosX+pos,DelocFoilPosY+pos,-DistTunnel));
+
+ AddDetector(G4ThreeVector(DelocFoilPosX - pos, DelocFoilPosY - pos, DistDSSD1),
+ G4ThreeVector(DelocFoilPosX - pos, DelocFoilPosY + pos, DistDSSD1),
+ G4ThreeVector(DelocFoilPosX + pos, DelocFoilPosY - pos, DistDSSD1),
+ G4ThreeVector(DelocFoilPosX + pos, DelocFoilPosY + pos, DistDSSD1));
+ AddDetector(G4ThreeVector(DelocFoilPosX - pos, DelocFoilPosY - pos, -DistTunnel),
+ G4ThreeVector(DelocFoilPosX - pos, DelocFoilPosY + pos, -DistTunnel),
+ G4ThreeVector(DelocFoilPosX + pos, DelocFoilPosY - pos, -DistTunnel),
+ G4ThreeVector(DelocFoilPosX + pos, DelocFoilPosY + pos, -DistTunnel));
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* SEASON::BuildSquareDetector(){
-
- G4Box* box = new G4Box("SEASON_Box",SEASON_NS::DetSize*0.5, SEASON_NS::DetSize*0.5,SEASON_NS::Thickness*0.5);
+G4LogicalVolume* SEASON::BuildSquareDetector() {
+
+ G4Box* box = new G4Box("SEASON_Box", SEASON_NS::DetSize * 0.5, SEASON_NS::DetSize * 0.5, SEASON_NS::Thickness * 0.5);
G4Material* DetectorMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Si");
- G4LogicalVolume* m_SquareDetector = new G4LogicalVolume(box,DetectorMaterial,"logic_SEASON_Box",0,0,0);
+ G4LogicalVolume* m_SquareDetector = new G4LogicalVolume(box, DetectorMaterial, "logic_SEASON_Box", 0, 0, 0);
m_SquareDetector->SetVisAttributes(m_VisSquare);
m_SquareDetector->SetSensitiveDetector(m_SEASONScorer);
return m_SquareDetector;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* SEASON::BuildWindow(){
-
- G4Box* solidWindow = new G4Box("WindowBox",SEASON_NS::DetSize*0.5, SEASON_NS::DetSize*0.5,SEASON_NS::DeadLayerThickness*0.5);
- G4LogicalVolume* logicWindow = new G4LogicalVolume(solidWindow,SEASON_NS::DeadLayerMaterial,"logicWindow",0,0,0);
- //logicWindow->SetVisAttributes(new G4VisAttributes(G4Colour(0, 1, 1, 0.5)));
+G4LogicalVolume* SEASON::BuildWindow() {
+
+ G4Box* solidWindow =
+ new G4Box("WindowBox", SEASON_NS::DetSize * 0.5, SEASON_NS::DetSize * 0.5, SEASON_NS::DeadLayerThickness * 0.5);
+ G4LogicalVolume* logicWindow = new G4LogicalVolume(solidWindow, SEASON_NS::DeadLayerMaterial, "logicWindow", 0, 0, 0);
+ // logicWindow->SetVisAttributes(new G4VisAttributes(G4Colour(0, 1, 1, 0.5)));
logicWindow->SetVisAttributes(m_VisWindow);
return logicWindow;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* SEASON::BuildChamber(){
- G4Box* plainChamber = new G4Box("PlainChamber",300*mm, 372.75*mm,4*mm);
-
+G4LogicalVolume* SEASON::BuildChamber() {
+ G4Box* plainChamber = new G4Box("PlainChamber", 300 * mm, 372.75 * mm, 4 * mm);
+
// make holes
- G4Tubs* hole = new G4Tubs("Hole", 0, 141.5*mm+1*micrometer, 4*mm+100*micrometer, 0*degree, 360*degree);
-
- G4SubtractionSolid * DrilledOnceChamber = new G4SubtractionSolid("DrilledOnceChamber",plainChamber,hole,NULL,G4ThreeVector(70.9*mm,-159.05*mm,1.5*mm));
- G4SubtractionSolid * DrilledChamber = new G4SubtractionSolid("DrilledChamber",DrilledOnceChamber,hole,NULL,G4ThreeVector(-50*mm,105.75*mm,1.5*mm));
-
-
- G4Material * ChamberMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
- G4LogicalVolume * m_Chamber = new G4LogicalVolume(DrilledChamber, ChamberMaterial, "logic_SEASON_Chamber");
+ G4Tubs* hole =
+ new G4Tubs("Hole", 0, 141.5 * mm + 1 * micrometer, 4 * mm + 100 * micrometer, 0 * degree, 360 * degree);
+
+ G4SubtractionSolid* DrilledOnceChamber = new G4SubtractionSolid("DrilledOnceChamber", plainChamber, hole, NULL,
+ G4ThreeVector(70.9 * mm, -159.05 * mm, 1.5 * mm));
+ G4SubtractionSolid* DrilledChamber = new G4SubtractionSolid("DrilledChamber", DrilledOnceChamber, hole, NULL,
+ G4ThreeVector(-50 * mm, 105.75 * mm, 1.5 * mm));
+
+ G4Material* ChamberMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ G4LogicalVolume* m_Chamber = new G4LogicalVolume(DrilledChamber, ChamberMaterial, "logic_SEASON_Chamber");
m_Chamber->SetVisAttributes(m_VisChamber);
-
+
return m_Chamber;
}
-G4LogicalVolume* SEASON::BuildGrid(){
- G4Box* plainStrip = new G4Box("plainStrip",SEASON_NS::StripPitch * 0.5, SEASON_NS::StripLength * 0.5, SEASON_NS::GridThickness * 0.5);
- G4Box* hole = new G4Box("hole", SEASON_NS::StripPitch * 0.5 - SEASON_NS::GridWidth, SEASON_NS::StripLength * 0.5 - SEASON_NS::GridWidth, SEASON_NS::GridThickness * 0.5 + 10 * um);
- G4SubtractionSolid * solidStrip = new G4SubtractionSolid("SEASON_Strip", plainStrip, hole, 0, G4ThreeVector(0,0,0));
-
- G4LogicalVolume * logicStrip = new G4LogicalVolume(solidStrip,SEASON_NS::Grid_Material, "logicStrip");
- //logicStrip->SetVisAttributes(new G4VisAttributes(G4Colour(1, 0, 0, 1)));
+G4LogicalVolume* SEASON::BuildGrid() {
+ G4Box* plainStrip = new G4Box("plainStrip", SEASON_NS::StripPitch * 0.5, SEASON_NS::StripLength * 0.5,
+ SEASON_NS::GridThickness * 0.5);
+ G4Box* hole =
+ new G4Box("hole", SEASON_NS::StripPitch * 0.5 - SEASON_NS::GridWidth,
+ SEASON_NS::StripLength * 0.5 - SEASON_NS::GridWidth, SEASON_NS::GridThickness * 0.5 + 10 * um);
+ G4SubtractionSolid* solidStrip = new G4SubtractionSolid("SEASON_Strip", plainStrip, hole, 0, G4ThreeVector(0, 0, 0));
+
+ G4LogicalVolume* logicStrip = new G4LogicalVolume(solidStrip, SEASON_NS::Grid_Material, "logicStrip");
+ // logicStrip->SetVisAttributes(new G4VisAttributes(G4Colour(1, 0, 0, 1)));
logicStrip->SetVisAttributes(m_VisGrid);
return logicStrip;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Wheel with the mounting implantation foils
-G4LogicalVolume* SEASON::BuildWheel()
-{
+G4LogicalVolume* SEASON::BuildWheel() {
G4double WheelRadius = SEASON_NS::WheelRadius;
G4double WheelThickness = SEASON_NS::WheelThickness;
G4double FinLenght = SEASON_NS::FinLenght;
@@ -227,43 +247,47 @@ G4LogicalVolume* SEASON::BuildWheel()
G4double FinThickness = SEASON_NS::FinThickness;
G4double FoilNbr = SEASON_NS::FoilNbr;
// Make wheel
- G4Tubs * plainWheel = new G4Tubs("PlainWheel", 0, WheelRadius, WheelThickness*0.5, 0*degree, 360*degree);
-
+ G4Tubs* plainWheel = new G4Tubs("PlainWheel", 0, WheelRadius, WheelThickness * 0.5, 0 * degree, 360 * degree);
+
// Make fins
- G4Box * plainFin = new G4Box("PlainFin", FinWidth*0.5, FinLenght*0.5+0.5*cm, FinThickness*0.5);
-
+ G4Box* plainFin = new G4Box("PlainFin", FinWidth * 0.5, FinLenght * 0.5 + 0.5 * cm, FinThickness * 0.5);
+
// make holes
- G4Tubs * hole = new G4Tubs("Hole", 0, m_FoilRadius+1*micrometer, FinThickness*0.5+100*micrometer, 0*degree, 360*degree);
-
- G4UnionSolid * wheel = new G4UnionSolid("SEASON_Wheel",plainWheel,plainFin,0,G4ThreeVector(0, WheelRadius+FinLenght*0.5-0.5*cm, -FinThickness*0.5));
- G4SubtractionSolid * wheel2 = new G4SubtractionSolid("SEASON_Wheel",wheel,hole,NULL,G4ThreeVector(0,SEASON_NS::FoilWheelCenterDist, -FinThickness*0.5));
-
- for ( int i = 1; i < FoilNbr; i++ )
- {
+ G4Tubs* hole = new G4Tubs("Hole", 0, m_FoilRadius + 1 * micrometer, FinThickness * 0.5 + 100 * micrometer, 0 * degree,
+ 360 * degree);
+
+ G4UnionSolid* wheel =
+ new G4UnionSolid("SEASON_Wheel", plainWheel, plainFin, 0,
+ G4ThreeVector(0, WheelRadius + FinLenght * 0.5 - 0.5 * cm, -FinThickness * 0.5));
+ G4SubtractionSolid* wheel2 = new G4SubtractionSolid(
+ "SEASON_Wheel", wheel, hole, NULL, G4ThreeVector(0, SEASON_NS::FoilWheelCenterDist, -FinThickness * 0.5));
+
+ for (int i = 1; i < FoilNbr; i++) {
G4RotationMatrix zRot;
- zRot.rotateZ(i*2*TMath::Pi()/FoilNbr*rad);
- G4UnionSolid * wheelbis = new G4UnionSolid("SEASON_Wheel", wheel2, plainFin,
- &zRot, G4ThreeVector(TMath::Sin(i*2*TMath::Pi()/FoilNbr)*(WheelRadius+FinLenght*0.5-0.5*cm),
- TMath::Cos(i*2*TMath::Pi()/FoilNbr)*(WheelRadius+FinLenght*0.5-0.5*cm), -FinThickness*0.5));
+ zRot.rotateZ(i * 2 * TMath::Pi() / FoilNbr * rad);
+ G4UnionSolid* wheelbis = new G4UnionSolid(
+ "SEASON_Wheel", wheel2, plainFin, &zRot,
+ G4ThreeVector(TMath::Sin(i * 2 * TMath::Pi() / FoilNbr) * (WheelRadius + FinLenght * 0.5 - 0.5 * cm),
+ TMath::Cos(i * 2 * TMath::Pi() / FoilNbr) * (WheelRadius + FinLenght * 0.5 - 0.5 * cm),
+ -FinThickness * 0.5));
wheel = wheelbis;
- G4SubtractionSolid * wheelbis2 = new G4SubtractionSolid("SEASON_Wheel", wheel, hole,
- NULL, G4ThreeVector(TMath::Sin(i*2*TMath::Pi()/FoilNbr)*SEASON_NS::FoilWheelCenterDist,
- TMath::Cos(i*2*TMath::Pi()/FoilNbr)*SEASON_NS::FoilWheelCenterDist, -FinThickness*0.5));
+ G4SubtractionSolid* wheelbis2 = new G4SubtractionSolid(
+ "SEASON_Wheel", wheel, hole, NULL,
+ G4ThreeVector(TMath::Sin(i * 2 * TMath::Pi() / FoilNbr) * SEASON_NS::FoilWheelCenterDist,
+ TMath::Cos(i * 2 * TMath::Pi() / FoilNbr) * SEASON_NS::FoilWheelCenterDist, -FinThickness * 0.5));
wheel2 = wheelbis2;
}
-
- G4Material * WheelMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
- G4LogicalVolume * m_Wheel = new G4LogicalVolume(wheel2, WheelMaterial, "logic_SEASON_Wheel");
+
+ G4Material* WheelMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ G4LogicalVolume* m_Wheel = new G4LogicalVolume(wheel2, WheelMaterial, "logic_SEASON_Wheel");
m_Wheel->SetVisAttributes(m_VisWheel);
-
+
return m_Wheel;
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Wheel with the mounting implantation foils
-G4LogicalVolume* SEASON::BuildRotatedWheel()
-{
+G4LogicalVolume* SEASON::BuildRotatedWheel() {
G4double WheelRadius = SEASON_NS::WheelRadius;
G4double WheelThickness = SEASON_NS::WheelThickness;
G4double FinLenght = SEASON_NS::FinLenght;
@@ -271,53 +295,64 @@ G4LogicalVolume* SEASON::BuildRotatedWheel()
G4double FinThickness = SEASON_NS::FinThickness;
G4double FoilNbr = SEASON_NS::FoilNbr;
// Make wheel
- G4Tubs * plainWheel = new G4Tubs("PlainWheel", 0, WheelRadius, WheelThickness*0.5, 0*degree, 360*degree);
-
+ G4Tubs* plainWheel = new G4Tubs("PlainWheel", 0, WheelRadius, WheelThickness * 0.5, 0 * degree, 360 * degree);
+
// Make fins
- G4Box * plainFin = new G4Box("PlainFin", FinWidth*0.5, FinLenght*0.5+0.5*cm, FinThickness*0.5);
-
+ G4Box* plainFin = new G4Box("PlainFin", FinWidth * 0.5, FinLenght * 0.5 + 0.5 * cm, FinThickness * 0.5);
+
// make holes
- G4Tubs * hole = new G4Tubs("Hole", 0, m_FoilRadius+1*micrometer, FinThickness*0.5+100*micrometer, 0*degree, 360*degree);
-
+ G4Tubs* hole = new G4Tubs("Hole", 0, m_FoilRadius + 1 * micrometer, FinThickness * 0.5 + 100 * micrometer, 0 * degree,
+ 360 * degree);
+
G4RotationMatrix zRot0;
- zRot0.rotateZ(TMath::Pi()/FoilNbr*rad);
-
- G4UnionSolid * wheel = new G4UnionSolid("SEASON_Wheel", plainWheel, plainFin, &zRot0, G4ThreeVector(TMath::Sin(TMath::Pi()/FoilNbr)*(WheelRadius+FinLenght*0.5-0.5*cm), TMath::Cos(TMath::Pi()/FoilNbr)*(WheelRadius+FinLenght*0.5-0.5*cm), -FinThickness));
-
- G4SubtractionSolid * wheel2 = new G4SubtractionSolid("SEASON_Wheel", wheel, hole, NULL, G4ThreeVector(TMath::Sin(TMath::Pi()/FoilNbr)*SEASON_NS::FoilWheelCenterDist, TMath::Cos(TMath::Pi()/FoilNbr)*SEASON_NS::FoilWheelCenterDist, -FinThickness));
-
- for ( int i = 1; i < FoilNbr; i++ )
- {
+ zRot0.rotateZ(TMath::Pi() / FoilNbr * rad);
+
+ G4UnionSolid* wheel = new G4UnionSolid(
+ "SEASON_Wheel", plainWheel, plainFin, &zRot0,
+ G4ThreeVector(TMath::Sin(TMath::Pi() / FoilNbr) * (WheelRadius + FinLenght * 0.5 - 0.5 * cm),
+ TMath::Cos(TMath::Pi() / FoilNbr) * (WheelRadius + FinLenght * 0.5 - 0.5 * cm), -FinThickness));
+
+ G4SubtractionSolid* wheel2 = new G4SubtractionSolid(
+ "SEASON_Wheel", wheel, hole, NULL,
+ G4ThreeVector(TMath::Sin(TMath::Pi() / FoilNbr) * SEASON_NS::FoilWheelCenterDist,
+ TMath::Cos(TMath::Pi() / FoilNbr) * SEASON_NS::FoilWheelCenterDist, -FinThickness));
+
+ for (int i = 1; i < FoilNbr; i++) {
G4RotationMatrix zRot;
- zRot.rotateZ((i+1*0.5)*2*TMath::Pi()/FoilNbr*rad);
- G4UnionSolid * wheelbis = new G4UnionSolid("SEASON_Wheel", wheel2, plainFin, &zRot, G4ThreeVector(TMath::Sin((i+1*0.5)*2*TMath::Pi()/FoilNbr)*(WheelRadius+FinLenght*0.5-0.5*cm), TMath::Cos((i+1*0.5)*2*TMath::Pi()/FoilNbr)*(WheelRadius+FinLenght*0.5-0.5*cm), -FinThickness));
+ zRot.rotateZ((i + 1 * 0.5) * 2 * TMath::Pi() / FoilNbr * rad);
+ G4UnionSolid* wheelbis = new G4UnionSolid(
+ "SEASON_Wheel", wheel2, plainFin, &zRot,
+ G4ThreeVector(
+ TMath::Sin((i + 1 * 0.5) * 2 * TMath::Pi() / FoilNbr) * (WheelRadius + FinLenght * 0.5 - 0.5 * cm),
+ TMath::Cos((i + 1 * 0.5) * 2 * TMath::Pi() / FoilNbr) * (WheelRadius + FinLenght * 0.5 - 0.5 * cm),
+ -FinThickness));
wheel = wheelbis;
- G4SubtractionSolid * wheelbis2 = new G4SubtractionSolid("SEASON_Wheel", wheel, hole, NULL, G4ThreeVector(TMath::Sin((i+1*0.5)*2*TMath::Pi()/FoilNbr)*SEASON_NS::FoilWheelCenterDist, TMath::Cos((i+1*0.5)*2*TMath::Pi()/FoilNbr)*SEASON_NS::FoilWheelCenterDist, -FinThickness));
+ G4SubtractionSolid* wheelbis2 = new G4SubtractionSolid(
+ "SEASON_Wheel", wheel, hole, NULL,
+ G4ThreeVector(TMath::Sin((i + 1 * 0.5) * 2 * TMath::Pi() / FoilNbr) * SEASON_NS::FoilWheelCenterDist,
+ TMath::Cos((i + 1 * 0.5) * 2 * TMath::Pi() / FoilNbr) * SEASON_NS::FoilWheelCenterDist,
+ -FinThickness));
wheel2 = wheelbis2;
}
-
- G4Material * WheelMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
- G4LogicalVolume * m_Wheel = new G4LogicalVolume(wheel2, WheelMaterial, "logic_SEASON_Wheel");
+
+ G4Material* WheelMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ G4LogicalVolume* m_Wheel = new G4LogicalVolume(wheel2, WheelMaterial, "logic_SEASON_Wheel");
m_Wheel->SetVisAttributes(m_VisWheel);
-
+
return m_Wheel;
}
-
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// The implantation foil (Carbon)
-G4LogicalVolume* SEASON::BuildFoil()
-{
- G4Tubs * foil = new G4Tubs("ImplantationFoil", 0, m_FoilRadius, m_FoilThickness * 0.5, 0*degree, 360*degree);
- G4LogicalVolume * m_Foil = new G4LogicalVolume(foil, SEASON_NS::FoilMaterial, "logic_SEASON_Foil");
-
- m_Foil->SetVisAttributes(m_VisFoil); //red G4Colour(0.5, 0.5, 0.5, 1)); //grey
-
+G4LogicalVolume* SEASON::BuildFoil() {
+ G4Tubs* foil = new G4Tubs("ImplantationFoil", 0, m_FoilRadius, m_FoilThickness * 0.5, 0 * degree, 360 * degree);
+ G4LogicalVolume* m_Foil = new G4LogicalVolume(foil, SEASON_NS::FoilMaterial, "logic_SEASON_Foil");
+
+ m_Foil->SetVisAttributes(m_VisFoil); // red G4Colour(0.5, 0.5, 0.5, 1)); //grey
+
return m_Foil;
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -325,112 +360,123 @@ G4LogicalVolume* SEASON::BuildFoil()
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void SEASON::ReadConfiguration(NPL::InputParser parser){
+void SEASON::ReadConfiguration(NPL::InputParser parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SEASON");
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
-
- vector<string> cart = {"X1_Y1","X1_YMax","XMax_Y1","XMax_YMax"};
- vector<string> cartfirst = {"X1_Y1","X1_YMax","XMax_Y1","XMax_YMax", "EnergyThreshold"};
- vector<string> full = {"DSSD1Dist", "TunnelDist", "UseDelocStation", "UseWheel", "UseChamber", "EnergyThreshold", "FoilRadius", "FoilThickness"};
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks.size() << " detectors found " << endl;
+
+ vector<string> cart = {"X1_Y1", "X1_YMax", "XMax_Y1", "XMax_YMax"};
+ vector<string> cartfirst = {"X1_Y1", "X1_YMax", "XMax_Y1", "XMax_YMax", "EnergyThreshold"};
+ vector<string> full = {"DSSD1Dist", "TunnelDist", "UseDelocStation", "UseWheel",
+ "UseChamber", "EnergyThreshold", "FoilRadius", "FoilThickness"};
double FoilThickness;
- for(unsigned int i = 0 ; i < blocks.size() ; i++){
- //cout << "!!! !!! !!! OK1 " << endl;
- if(blocks[i]->HasTokenList(cartfirst)){
- //cout << "!!! !!! !!! Cart ok" << endl;
- if(NPOptionManager::getInstance()->GetVerboseLevel()) cout << endl << "//// SEASON " << i+1 << endl;
-
- G4ThreeVector A = NPS::ConvertVector( blocks[i]->GetTVector3("X1_Y1","mm"));
- G4ThreeVector B = NPS::ConvertVector( blocks[i]->GetTVector3("X1_YMax","mm"));
- G4ThreeVector C = NPS::ConvertVector( blocks[i]->GetTVector3("XMax_Y1","mm"));
- G4ThreeVector D = NPS::ConvertVector( blocks[i]->GetTVector3("XMax_YMax","mm"));
- m_EnergyThreshold = blocks[i]->GetDouble("EnergyThreshold","keV");
- m_FoilThickness = (FoilThickness/0.226) * nm;
- AddDetector(A,B,C,D);
+ for (unsigned int i = 0; i < blocks.size(); i++) {
+ // cout << "!!! !!! !!! OK1 " << endl;
+ if (blocks[i]->HasTokenList(cartfirst)) {
+ // cout << "!!! !!! !!! Cart ok" << endl;
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// SEASON " << i + 1 << endl;
+
+ G4ThreeVector A = NPS::ConvertVector(blocks[i]->GetTVector3("X1_Y1", "mm"));
+ G4ThreeVector B = NPS::ConvertVector(blocks[i]->GetTVector3("X1_YMax", "mm"));
+ G4ThreeVector C = NPS::ConvertVector(blocks[i]->GetTVector3("XMax_Y1", "mm"));
+ G4ThreeVector D = NPS::ConvertVector(blocks[i]->GetTVector3("XMax_YMax", "mm"));
+ m_EnergyThreshold = blocks[i]->GetDouble("EnergyThreshold", "keV");
+ m_FoilThickness = (FoilThickness / 0.226) * nm;
+ AddDetector(A, B, C, D);
}
- else if(blocks[i]->HasTokenList(cart)){
- //cout << "!!! !!! !!! Cart ok" << endl;
- if(NPOptionManager::getInstance()->GetVerboseLevel()) cout << endl << "//// SEASON " << i+1 << endl;
-
- G4ThreeVector A = NPS::ConvertVector( blocks[i]->GetTVector3("X1_Y1","mm"));
- G4ThreeVector B = NPS::ConvertVector( blocks[i]->GetTVector3("X1_YMax","mm"));
- G4ThreeVector C = NPS::ConvertVector( blocks[i]->GetTVector3("XMax_Y1","mm"));
- G4ThreeVector D = NPS::ConvertVector( blocks[i]->GetTVector3("XMax_YMax","mm"));
- AddDetector(A,B,C,D);
+ else if (blocks[i]->HasTokenList(cart)) {
+ // cout << "!!! !!! !!! Cart ok" << endl;
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// SEASON " << i + 1 << endl;
+
+ G4ThreeVector A = NPS::ConvertVector(blocks[i]->GetTVector3("X1_Y1", "mm"));
+ G4ThreeVector B = NPS::ConvertVector(blocks[i]->GetTVector3("X1_YMax", "mm"));
+ G4ThreeVector C = NPS::ConvertVector(blocks[i]->GetTVector3("XMax_Y1", "mm"));
+ G4ThreeVector D = NPS::ConvertVector(blocks[i]->GetTVector3("XMax_YMax", "mm"));
+ AddDetector(A, B, C, D);
}
- else if(blocks[i]->HasTokenList(full)){
- if(NPOptionManager::getInstance()->GetVerboseLevel()) cout << endl << "//// SEASON " << i+1 << endl;
-
- double DistDSSD1 = blocks[i]->GetDouble("DSSD1Dist","mm");
- double DistTunnel = blocks[i]->GetDouble("TunnelDist","mm");
-
- m_EnergyThreshold = blocks[i]->GetDouble("EnergyThreshold","keV");
- m_FoilRadius = blocks[i]->GetDouble("FoilRadius","mm");
- FoilThickness = blocks[i]->GetDouble("FoilThickness","void");
- m_FoilThickness = (FoilThickness/0.226) * nm;
-
+ else if (blocks[i]->HasTokenList(full)) {
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// SEASON " << i + 1 << endl;
+
+ double DistDSSD1 = blocks[i]->GetDouble("DSSD1Dist", "mm");
+ double DistTunnel = blocks[i]->GetDouble("TunnelDist", "mm");
+
+ m_EnergyThreshold = blocks[i]->GetDouble("EnergyThreshold", "keV");
+ m_FoilRadius = blocks[i]->GetDouble("FoilRadius", "mm");
+ FoilThickness = blocks[i]->GetDouble("FoilThickness", "void");
+ m_FoilThickness = (FoilThickness / 0.226) * nm;
+
string Deloc = blocks[i]->GetString("UseDelocStation");
bool EnableDeloc;
- if(Deloc=="True") EnableDeloc = true;
- else if(Deloc=="False") EnableDeloc = false;
- else{
+ if (Deloc == "True")
+ EnableDeloc = true;
+ else if (Deloc == "False")
+ EnableDeloc = false;
+ else {
cout << "ERROR: UseDelocStation must be either True or False, please check input file formatting" << endl;
exit(1);
}
-
+
string Wheel = blocks[i]->GetString("UseWheel");
bool EnableWheel = false;
bool RotateWheel = false;
- if(Wheel=="Rotated") RotateWheel = true;
- else if(Wheel=="True") EnableWheel = true;
- else if(Wheel=="False") EnableWheel = false;
- else{
+ if (Wheel == "Rotated")
+ RotateWheel = true;
+ else if (Wheel == "True")
+ EnableWheel = true;
+ else if (Wheel == "False")
+ EnableWheel = false;
+ else {
cout << "ERROR: UseWheel must be either True, False or Rotated, please check input file formatting" << endl;
exit(1);
}
-
+
string Chamber = blocks[i]->GetString("UseChamber");
bool EnableChamber;
- if(Chamber=="True") EnableChamber = true;
- else if(Chamber=="False") EnableChamber = false;
- else{
+ if (Chamber == "True")
+ EnableChamber = true;
+ else if (Chamber == "False")
+ EnableChamber = false;
+ else {
cout << "ERROR: UseChamber must be either True or False, please check input file formatting" << endl;
exit(1);
}
-
+
AddDetector(DistDSSD1, DistTunnel, EnableDeloc, EnableWheel, RotateWheel, EnableChamber);
}
- else{
+ else {
cout << "ERROR: check your input file formatting " << endl;
exit(1);
}
}
- cout << " Foil Thickness : " << FoilThickness << " -> " << m_FoilThickness/nm << " nm with a carbon density rho = 2.26 g/cm3" << endl;
+ cout << " Foil Thickness : " << FoilThickness << " -> " << m_FoilThickness / nm
+ << " nm with a carbon density rho = 2.26 g/cm3" << endl;
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
-void SEASON::ConstructDetector(G4LogicalVolume* world){
-
+void SEASON::ConstructDetector(G4LogicalVolume* world) {
+
G4bool checkOverlaps = true;
-
+
//
// Construct DSSDs
-
- for (unsigned short i = 0 ; i < m_X1_Y1.size() ; i++) {
-
- G4RotationMatrix* Rot = NULL;
- G4ThreeVector Det_pos = G4ThreeVector(0, 0, 0);
- G4ThreeVector u = G4ThreeVector(0, 0, 0);
- G4ThreeVector v = G4ThreeVector(0, 0, 0);
- G4ThreeVector w = G4ThreeVector(0, 0, 0);
- G4ThreeVector Center = G4ThreeVector(0, 0, 0);
-
- G4ThreeVector Front_Window_pos = G4ThreeVector(0, 0, 0);
- G4ThreeVector Grid_pos = G4ThreeVector(0, 0, 0);
+
+ for (unsigned short i = 0; i < m_X1_Y1.size(); i++) {
+
+ G4RotationMatrix* Rot = NULL;
+ G4ThreeVector Det_pos = G4ThreeVector(0, 0, 0);
+ G4ThreeVector u = G4ThreeVector(0, 0, 0);
+ G4ThreeVector v = G4ThreeVector(0, 0, 0);
+ G4ThreeVector w = G4ThreeVector(0, 0, 0);
+ G4ThreeVector Center = G4ThreeVector(0, 0, 0);
+
+ G4ThreeVector Front_Window_pos = G4ThreeVector(0, 0, 0);
+ G4ThreeVector Grid_pos = G4ThreeVector(0, 0, 0);
// (u,v,w) unitary vector associated to telescope referencial
// (u,v) // to silicon plan
// w perpendicular to (u,v) plan and pointing CsI
@@ -440,191 +486,212 @@ void SEASON::ConstructDetector(G4LogicalVolume* world){
v = v.unit();
w = u.cross(v);
w = w.unit();
-
+
Center = (m_X1_Y1[i] + m_X1_YMax[i] + m_XMax_Y1[i] + m_XMax_YMax[i]) / 4;
-
+
// Passage Matrix from Lab Referential to Telescope Referential
Rot = new G4RotationMatrix(u, v, w);
- if( w == G4ThreeVector(0,0,-1) && Center.getZ() > 0) w = -w;
- Det_pos = w * SEASON_NS::Thickness * 0.5 + Center ;
-
- Front_Window_pos = - w * SEASON_NS::DeadLayerThickness * 0.5 + Center;
- Grid_pos = - w * (SEASON_NS::DeadLayerThickness + SEASON_NS::GridThickness * 0.5) + Center;
-
- if(SEASON_NS::DeadLayerThickness>0){
- new G4PVPlacement(G4Transform3D(*Rot,Front_Window_pos), BuildWindow(), "SEASON",world,false, m_X1_Y1.size()+i+1, FALSE);
+ if (w == G4ThreeVector(0, 0, -1) && Center.getZ() > 0)
+ w = -w;
+ Det_pos = w * SEASON_NS::Thickness * 0.5 + Center;
+
+ Front_Window_pos = -w * SEASON_NS::DeadLayerThickness * 0.5 + Center;
+ Grid_pos = -w * (SEASON_NS::DeadLayerThickness + SEASON_NS::GridThickness * 0.5) + Center;
+
+ if (SEASON_NS::DeadLayerThickness > 0) {
+ new G4PVPlacement(G4Transform3D(*Rot, Front_Window_pos), BuildWindow(), "SEASON", world, false,
+ m_X1_Y1.size() + i + 1, FALSE);
}
- else if(i==0) cout << "No entrance window " << endl;
-
-
- if(SEASON_NS::GridThickness>0 and SEASON_NS::GridWidth>0){
- for(int strip=0;strip<SEASON_NS::NumberOfStripsX;strip++){
- G4ThreeVector Strip_pos = u*((-strip+(SEASON_NS::NumberOfStripsX-1)/2)*SEASON_NS::StripWidth);
- new G4PVPlacement(G4Transform3D(*Rot,Grid_pos+Strip_pos), BuildGrid(), "SEASON",world,false, 2*m_X1_Y1.size()+i+1, FALSE);
+ else if (i == 0)
+ cout << "No entrance window " << endl;
+
+ if (SEASON_NS::GridThickness > 0 and SEASON_NS::GridWidth > 0) {
+ for (int strip = 0; strip < SEASON_NS::NumberOfStripsX; strip++) {
+ G4ThreeVector Strip_pos = u * ((-strip + (SEASON_NS::NumberOfStripsX - 1) / 2) * SEASON_NS::StripWidth);
+ new G4PVPlacement(G4Transform3D(*Rot, Grid_pos + Strip_pos), BuildGrid(), "SEASON", world, false,
+ 2 * m_X1_Y1.size() + i + 1, FALSE);
}
}
- else if (i==0) cout << "No aluminium grid" << endl;
-
- new G4PVPlacement(G4Transform3D(*Rot,Det_pos), BuildSquareDetector(),"SEASON",world,false,i+1, checkOverlaps);
-
+ else if (i == 0)
+ cout << "No aluminium grid" << endl;
+
+ new G4PVPlacement(G4Transform3D(*Rot, Det_pos), BuildSquareDetector(), "SEASON", world, false, i + 1,
+ checkOverlaps);
}
-
-
- if(m_RotateWheel){
- //
+
+ if (m_RotateWheel) {
+ //
// Build Wheel
G4RotationMatrix zRot;
- zRot.rotateZ(TMath::Pi()*rad);
-
- new G4PVPlacement(G4Transform3D(zRot, G4ThreeVector(0, SEASON_NS::FoilWheelCenterDist, 0)), BuildRotatedWheel(), "Wheel", world, false, 0, checkOverlaps);
-
+ zRot.rotateZ(TMath::Pi() * rad);
+
+ new G4PVPlacement(G4Transform3D(zRot, G4ThreeVector(0, SEASON_NS::FoilWheelCenterDist, 0)), BuildRotatedWheel(),
+ "Wheel", world, false, 0, checkOverlaps);
+
// Build foils in the holes
- for ( int i = 0; i < SEASON_NS::FoilNbr; i++ )
- {
- new G4PVPlacement(G4Transform3D(G4RotationMatrix(0, 0, 0), G4ThreeVector(TMath::Sin((i+1*0.5)*2*TMath::Pi()/SEASON_NS::FoilNbr + TMath::Pi()) * SEASON_NS::FoilWheelCenterDist, TMath::Cos((i+1*0.5)*2*TMath::Pi() / SEASON_NS::FoilNbr + TMath::Pi()) * SEASON_NS::FoilWheelCenterDist + SEASON_NS::FoilWheelCenterDist, m_FoilThickness * 0.5)), BuildFoil(), "Foil", world, false, 0, checkOverlaps);
+ for (int i = 0; i < SEASON_NS::FoilNbr; i++) {
+ new G4PVPlacement(
+ G4Transform3D(G4RotationMatrix(0, 0, 0),
+ G4ThreeVector(TMath::Sin((i + 1 * 0.5) * 2 * TMath::Pi() / SEASON_NS::FoilNbr + TMath::Pi()) *
+ SEASON_NS::FoilWheelCenterDist,
+ TMath::Cos((i + 1 * 0.5) * 2 * TMath::Pi() / SEASON_NS::FoilNbr + TMath::Pi()) *
+ SEASON_NS::FoilWheelCenterDist +
+ SEASON_NS::FoilWheelCenterDist,
+ m_FoilThickness * 0.5)),
+ BuildFoil(), "Foil", world, false, 0, checkOverlaps);
}
- }
- else if(m_EnableWheel){
- //
+ }
+ else if (m_EnableWheel) {
+ //
// Build Wheel
G4RotationMatrix zRot;
- zRot.rotateZ(TMath::Pi()*rad);
-
- new G4PVPlacement(G4Transform3D(zRot, G4ThreeVector(0, SEASON_NS::FoilWheelCenterDist, 0)), BuildWheel(), "Wheel", world, false, 0, checkOverlaps);
-
+ zRot.rotateZ(TMath::Pi() * rad);
+
+ new G4PVPlacement(G4Transform3D(zRot, G4ThreeVector(0, SEASON_NS::FoilWheelCenterDist, 0)), BuildWheel(), "Wheel",
+ world, false, 0, checkOverlaps);
+
// Build foils in the holes
- for ( int i = 0; i < SEASON_NS::FoilNbr; i++ )
- {
- new G4PVPlacement(G4Transform3D(G4RotationMatrix(0, 0, 0), G4ThreeVector(TMath::Sin(i*2*TMath::Pi()/SEASON_NS::FoilNbr + TMath::Pi()) * SEASON_NS::FoilWheelCenterDist, TMath::Cos(i*2*TMath::Pi()/SEASON_NS::FoilNbr + TMath::Pi()) * SEASON_NS::FoilWheelCenterDist + SEASON_NS::FoilWheelCenterDist, m_FoilThickness * 0.5)), BuildFoil(), "Foil", world, false, 0, checkOverlaps);
+ for (int i = 0; i < SEASON_NS::FoilNbr; i++) {
+ new G4PVPlacement(G4Transform3D(G4RotationMatrix(0, 0, 0),
+ G4ThreeVector(TMath::Sin(i * 2 * TMath::Pi() / SEASON_NS::FoilNbr + TMath::Pi()) *
+ SEASON_NS::FoilWheelCenterDist,
+ TMath::Cos(i * 2 * TMath::Pi() / SEASON_NS::FoilNbr + TMath::Pi()) *
+ SEASON_NS::FoilWheelCenterDist +
+ SEASON_NS::FoilWheelCenterDist,
+ m_FoilThickness * 0.5)),
+ BuildFoil(), "Foil", world, false, 0, checkOverlaps);
}
}
-
- if(m_EnableChamber){
- //
+
+ if (m_EnableChamber) {
+ //
// Build Wheel
- G4RotationMatrix zRot;
- new G4PVPlacement(G4Transform3D(zRot, G4ThreeVector(-70.9*mm,159.05*mm,m_DistDSSD1+SEASON_NS::DistDSSD1_AlWindow)), BuildChamber(), "Chamber", world, false, 0, checkOverlaps);
- }
+ G4RotationMatrix zRot;
+ new G4PVPlacement(
+ G4Transform3D(zRot, G4ThreeVector(-70.9 * mm, 159.05 * mm, m_DistDSSD1 + SEASON_NS::DistDSSD1_AlWindow)),
+ BuildChamber(), "Chamber", world, false, 0, checkOverlaps);
+ }
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method
-void SEASON::InitializeRootOutput(){
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- if(!pTree->FindBranch("SEASON")){
- pTree->Branch("SEASON", "TSEASONData", &m_Event) ;
+void SEASON::InitializeRootOutput() {
+ RootOutput* pAnalysis = RootOutput::getInstance();
+ TTree* pTree = pAnalysis->GetTree();
+ if (!pTree->FindBranch("SEASON")) {
+ pTree->Branch("SEASON", "TSEASONData", &m_Event);
}
- pTree->SetBranchAddress("SEASON", &m_Event) ;
+ pTree->SetBranchAddress("SEASON", &m_Event);
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
-void SEASON::ReadSensitive(const G4Event* event){
+void SEASON::ReadSensitive(const G4Event* event) {
m_Event->Clear();
-
+
///////////
- // Strip scorer
- DSSDScorers::PS_Images * SiScorer= (DSSDScorers::PS_Images*) m_SEASONScorer->GetPrimitive(0);
- map< unsigned int, pair<double,double> > mapFront;
- map< unsigned int, unsigned short int> mapRFront;
- map< unsigned int, pair<double,double> >::iterator it;
-
+ // Strip scorer
+ DSSDScorers::PS_Images* SiScorer = (DSSDScorers::PS_Images*)m_SEASONScorer->GetPrimitive(0);
+ map<unsigned int, pair<double, double>> mapFront;
+ map<unsigned int, unsigned short int> mapRFront;
+ map<unsigned int, pair<double, double>>::iterator it;
+
// X-side reading
- for(unsigned int i = 0 ; i < SiScorer->GetFrontMult() ; i++){
+ for (unsigned int i = 0; i < SiScorer->GetFrontMult(); i++) {
double EnergyX = SiScorer->GetEnergyFront(i);
unsigned int DetNbrX = SiScorer->GetDetectorFront(i);
double TimeX = SiScorer->GetTimeFront(i);
-
- unsigned int a,r,g,b;
- SiScorer->GetARGBFront(i,a,r,g,b);
-
- if(r==0)
- {
- mapFront[b+DetNbrX*1e6].first+=EnergyX;
- mapFront[b+DetNbrX*1e6].second+=TimeX;
- mapRFront[b+DetNbrX*1e6]+=r+i;
+
+ unsigned int a, r, g, b;
+ SiScorer->GetARGBFront(i, a, r, g, b);
+
+ if (r == 0) {
+ mapFront[b + DetNbrX * 1e6].first += EnergyX;
+ mapFront[b + DetNbrX * 1e6].second += TimeX;
+ mapRFront[b + DetNbrX * 1e6] += r + i;
}
- else{
+ else {
double rand = G4UniformRand();
- if (rand > 0.5){
+ if (rand > 0.5) {
double energy1 = EnergyX * rand;
- double energy2 = EnergyX * (rand-1);
- mapFront[b+DetNbrX*1e6].first+=energy1;
- mapFront[b+DetNbrX*1e6].second=TimeX;
- mapRFront[b+DetNbrX*1e6]+=r+i;
- mapFront[g+DetNbrX*1e6].first+=energy2;
- mapFront[g+DetNbrX*1e6].second=TimeX;
- mapRFront[g+DetNbrX*1e6]+=r+i;
+ double energy2 = EnergyX * (rand - 1);
+ mapFront[b + DetNbrX * 1e6].first += energy1;
+ mapFront[b + DetNbrX * 1e6].second = TimeX;
+ mapRFront[b + DetNbrX * 1e6] += r + i;
+ mapFront[g + DetNbrX * 1e6].first += energy2;
+ mapFront[g + DetNbrX * 1e6].second = TimeX;
+ mapRFront[g + DetNbrX * 1e6] += r + i;
}
else {
- double energy1 = - EnergyX * rand;
- double energy2 = EnergyX * (1-rand);
- mapFront[b+DetNbrX*1e6].first+=energy1;
- mapFront[b+DetNbrX*1e6].second=TimeX;
- mapRFront[b+DetNbrX*1e6]+=r+i;
- mapFront[g+DetNbrX*1e6].first+=energy2;
- mapFront[g+DetNbrX*1e6].second=TimeX;
- mapRFront[g+DetNbrX*1e6]+=r+i;
+ double energy1 = -EnergyX * rand;
+ double energy2 = EnergyX * (1 - rand);
+ mapFront[b + DetNbrX * 1e6].first += energy1;
+ mapFront[b + DetNbrX * 1e6].second = TimeX;
+ mapRFront[b + DetNbrX * 1e6] += r + i;
+ mapFront[g + DetNbrX * 1e6].first += energy2;
+ mapFront[g + DetNbrX * 1e6].second = TimeX;
+ mapRFront[g + DetNbrX * 1e6] += r + i;
}
}
}
-
- for(it=mapFront.begin();it!=mapFront.end();it++){
+
+ for (it = mapFront.begin(); it != mapFront.end(); it++) {
double EnergyX = it->second.first;
- if(abs(EnergyX) > m_EnergyThreshold ){
- if(EnergyX>2.) EnergyX = RandGauss::shoot(EnergyX,SEASON_NS::ResoAlpha);
- else EnergyX = RandGauss::shoot(EnergyX,SEASON_NS::ResoElectron);
+ if (abs(EnergyX) > m_EnergyThreshold) {
+ if (EnergyX > 2.)
+ EnergyX = RandGauss::shoot(EnergyX, SEASON_NS::ResoAlpha);
+ else
+ EnergyX = RandGauss::shoot(EnergyX, SEASON_NS::ResoElectron);
double TimeX = RandGauss::shoot(it->second.second, SEASON_NS::ResoTime);
- unsigned int strip = it->first-1000000*(it->first/1000000);
- unsigned int det = it->first/1000000;
+ unsigned int strip = it->first - 1000000 * (it->first / 1000000);
+ unsigned int det = it->first / 1000000;
m_Event->SetXEnergy(det, strip, EnergyX);
m_Event->SetXTime(det, strip, TimeX);
m_Event->SetXParticleID(mapRFront[it->first]);
}
}
-
- map< unsigned int, pair<double,double> > mapBack;
- map< unsigned int, unsigned short int> mapRBack;
+
+ map<unsigned int, pair<double, double>> mapBack;
+ map<unsigned int, unsigned short int> mapRBack;
// Y-side reading
- for(unsigned int i = 0 ; i < SiScorer->GetBackMult() ; i++){
+ for (unsigned int i = 0; i < SiScorer->GetBackMult(); i++) {
double EnergyY = SiScorer->GetEnergyBack(i);
- unsigned int DetNbrY = SiScorer->GetDetectorBack(i);
+ unsigned int DetNbrY = SiScorer->GetDetectorBack(i);
double TimeY = SiScorer->GetTimeBack(i);
-
- unsigned int a,r,g,b;
- SiScorer->GetARGBBack(i,a,r,g,b);
-
- if(r==0){
- mapBack[b+DetNbrY*1e6].first+=EnergyY;
- mapBack[b+DetNbrY*1e6].second+=TimeY;
- mapRBack[b+DetNbrY*1e6]+=r+i;
+
+ unsigned int a, r, g, b;
+ SiScorer->GetARGBBack(i, a, r, g, b);
+
+ if (r == 0) {
+ mapBack[b + DetNbrY * 1e6].first += EnergyY;
+ mapBack[b + DetNbrY * 1e6].second += TimeY;
+ mapRBack[b + DetNbrY * 1e6] += r + i;
}
- else{
+ else {
double rand = G4UniformRand();
double energy1 = rand * EnergyY;
- double energy2 = (1-rand) * EnergyY;
- mapBack[b+DetNbrY*1e6].first+=energy1;
- mapBack[b+DetNbrY*1e6].second+=TimeY;
- mapRBack[b+DetNbrY*1e6]+=r+i;
- mapBack[g+DetNbrY*1e6].first+=energy2;
- mapBack[g+DetNbrY*1e6].second+=TimeY;
- mapRBack[g+DetNbrY*1e6]+=r+i;
+ double energy2 = (1 - rand) * EnergyY;
+ mapBack[b + DetNbrY * 1e6].first += energy1;
+ mapBack[b + DetNbrY * 1e6].second += TimeY;
+ mapRBack[b + DetNbrY * 1e6] += r + i;
+ mapBack[g + DetNbrY * 1e6].first += energy2;
+ mapBack[g + DetNbrY * 1e6].second += TimeY;
+ mapRBack[g + DetNbrY * 1e6] += r + i;
}
}
-
- for(it=mapBack.begin();it!=mapBack.end();it++){
+
+ for (it = mapBack.begin(); it != mapBack.end(); it++) {
double EnergyY = it->second.first;
- if(EnergyY > m_EnergyThreshold ){
- if(EnergyY>2.) EnergyY = RandGauss::shoot(EnergyY,SEASON_NS::ResoAlpha);
- else EnergyY = RandGauss::shoot(EnergyY,SEASON_NS::ResoElectron);
+ if (EnergyY > m_EnergyThreshold) {
+ if (EnergyY > 2.)
+ EnergyY = RandGauss::shoot(EnergyY, SEASON_NS::ResoAlpha);
+ else
+ EnergyY = RandGauss::shoot(EnergyY, SEASON_NS::ResoElectron);
double TimeY = RandGauss::shoot(it->second.second, SEASON_NS::ResoTime);
- unsigned int strip = it->first-1000000*(it->first/1000000);
- unsigned int det = it->first/1000000;
+ unsigned int strip = it->first - 1000000 * (it->first / 1000000);
+ unsigned int det = it->first / 1000000;
m_Event->SetYEnergy(det, strip, EnergyY);
m_Event->SetYTime(det, strip, TimeY);
m_Event->SetYParticleID(mapRBack[it->first]);
@@ -634,28 +701,30 @@ void SEASON::ReadSensitive(const G4Event* event){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-////////////////////////////////////////////////////////////////
-void SEASON::InitializeScorers() {
+////////////////////////////////////////////////////////////////
+void SEASON::InitializeScorers() {
// This check is necessary in case the geometry is reloaded
- bool already_exist = false;
-
- if(already_exist)
- return ;
-
- m_SEASONScorer = CheckScorer("SEASONScorer",already_exist);
-
+ bool already_exist = false;
+
+ if (already_exist)
+ return;
+
+ m_SEASONScorer = CheckScorer("SEASONScorer", already_exist);
+
cout << "Size of scorers :\n";
cout << SEASON_NS::DetSize << " " << SEASON_NS::DetSize;
- cout << " " << SEASON_NS::NumberOfStripsX << " " << SEASON_NS::NumberOfStripsY << "\n";
-
-
+ cout << " " << SEASON_NS::NumberOfStripsX << " " << SEASON_NS::NumberOfStripsY << "\n";
+
string nptool = getenv("NPTOOL");
- G4VPrimitiveScorer* Scorer = new DSSDScorers::PS_Images("Scorer", nptool + "/NPLib/Detectors/SEASON/ressources/maskFront.png", nptool + "/NPLib/Detectors/SEASON/ressources/maskBack.png", 0.005, 0.005, 0, 0, 0xffff0000, 0);
-
- m_SEASONScorer->RegisterPrimitive(Scorer);
-
+ G4VPrimitiveScorer* Scorer = new DSSDScorers::PS_Images(
+ "Scorer", nptool + "/NPLib/Detectors/SEASON/ressources/maskFront.png",
+ nptool + "/NPLib/Detectors/SEASON/ressources/maskBack.png", 0.005, 0.005, 0, 0, 0xffff0000, 0);
+
+ m_SEASONScorer->RegisterPrimitive(Scorer);
+
// Interaction scorer
- m_SEASONScorer->RegisterPrimitive(new InteractionScorers::PS_Interactions("InteractionScorerSEASON", ms_InterCoord, 0));
+ m_SEASONScorer->RegisterPrimitive(
+ new InteractionScorers::PS_Interactions("InteractionScorerSEASON", ms_InterCoord, 0));
G4SDManager::GetSDMpointer()->AddNewDetector(m_SEASONScorer);
}
@@ -665,22 +734,20 @@ void SEASON::InitializeScorers() {
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
-NPS::VDetector* SEASON::Construct(){
- return (NPS::VDetector*) new SEASON();
-}
+NPS::VDetector* SEASON::Construct() { return (NPS::VDetector*)new SEASON(); }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
-extern"C" {
- class proxy_nps_SEASON{
- public:
- proxy_nps_SEASON(){
- NPS::DetectorFactory::getInstance()->AddToken("SEASON","SEASON");
- NPS::DetectorFactory::getInstance()->AddDetector("SEASON",SEASON::Construct);
- }
- };
-
- proxy_nps_SEASON p_nps_SEASON;
+extern "C" {
+class proxy_nps_SEASON {
+ public:
+ proxy_nps_SEASON() {
+ NPS::DetectorFactory::getInstance()->AddToken("SEASON", "SEASON");
+ NPS::DetectorFactory::getInstance()->AddDetector("SEASON", SEASON::Construct);
+ }
+};
+
+proxy_nps_SEASON p_nps_SEASON;
}
diff --git a/NPSimulation/Process/Decay.cc b/NPSimulation/Process/Decay.cc
index 30c82cdd4cfc5a142ba3456f28abc6608a4567b1..58be2c37c4e1bd5a71d6f6170b5ca4225f148993 100644
--- a/NPSimulation/Process/Decay.cc
+++ b/NPSimulation/Process/Decay.cc
@@ -21,67 +21,62 @@
* *
*****************************************************************************/
-#include <iostream>
-#include <string>
-#include <set>
-#include "NPFunction.h"
#include "Decay.hh"
-#include "NPOptionManager.h"
-#include "NPInputParser.h"
-#include "G4VPhysicalVolume.hh"
#include "G4Electron.hh"
#include "G4Gamma.hh"
-#include "G4SystemOfUnits.hh"
-#include "G4ParticleTable.hh"
#include "G4IonTable.hh"
-
+#include "G4ParticleTable.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4VPhysicalVolume.hh"
+#include "NPFunction.h"
+#include "NPInputParser.h"
+#include "NPOptionManager.h"
+#include <iostream>
+#include <set>
+#include <string>
using namespace NPS;
////////////////////////////////////////////////////////////////////////////////
-Decay::Decay(G4String modelName,G4Region* envelope) :
- G4VFastSimulationModel(modelName, envelope) {
- ReadConfiguration();
- m_PreviousEnergy=0 ;
- m_PreviousLength=0 ;
- }
-
+Decay::Decay(G4String modelName, G4Region* envelope) : G4VFastSimulationModel(modelName, envelope) {
+ ReadConfiguration();
+ m_PreviousEnergy = 0;
+ m_PreviousLength = 0;
+}
////////////////////////////////////////////////////////////////////////////////
-Decay::Decay(G4String modelName) :
- G4VFastSimulationModel(modelName) {
- }
+Decay::Decay(G4String modelName) : G4VFastSimulationModel(modelName) {}
////////////////////////////////////////////////////////////////////////////////
-Decay::~Decay() {
-}
+Decay::~Decay() {}
////////////////////////////////////////////////////////////////////////////////
-void Decay::ReadConfiguration(){
+void Decay::ReadConfiguration() {
NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
m_Decay.ReadConfiguration(input);
std::set<std::string> Mother = m_Decay.GetAllMotherName();
- std::set<std::string>::iterator it ;
- for(it = Mother.begin() ; it != Mother.end() ; it++){
+ std::set<std::string>::iterator it;
+ for (it = Mother.begin(); it != Mother.end(); it++) {
// ground state name, e.g. deuteron
m_MotherName.insert(NPL::ChangeNameToG4Standard(*it));
// excited state name e.g. H2
- m_MotherName.insert(NPL::ChangeNameToG4Standard(*it,true));
+ m_MotherName.insert(NPL::ChangeNameToG4Standard(*it, true));
}
}
////////////////////////////////////////////////////////////////////////////////
-G4bool Decay::IsApplicable( const G4ParticleDefinition& particleType) {
+G4bool Decay::IsApplicable(const G4ParticleDefinition& particleType) {
m_CurrentName = particleType.GetParticleName();
// Extract Ex from name
- if(m_CurrentName.find("[")!=std::string::npos)
- m_ExcitationEnergy = atof(m_CurrentName.substr(m_CurrentName.find("[")+1,m_CurrentName.find("]")-1).c_str())*keV;
+ if (m_CurrentName.find("[") != std::string::npos)
+ m_ExcitationEnergy =
+ atof(m_CurrentName.substr(m_CurrentName.find("[") + 1, m_CurrentName.find("]") - 1).c_str()) * keV;
else
- m_ExcitationEnergy=0;
+ m_ExcitationEnergy = 0;
// Strip name from excitation energy
- m_CurrentName = m_CurrentName.substr(0,m_CurrentName.find("["));
+ m_CurrentName = m_CurrentName.substr(0, m_CurrentName.find("["));
// If the decay exist
- if (m_MotherName.find(m_CurrentName)!=m_MotherName.end()) {
+ if (m_MotherName.find(m_CurrentName) != m_MotherName.end()) {
return true;
}
return false;
@@ -89,14 +84,14 @@ G4bool Decay::IsApplicable( const G4ParticleDefinition& particleType) {
////////////////////////////////////////////////////////////////////////////////
G4bool Decay::ModelTrigger(const G4FastTrack& fastTrack) {
- //FIXME: Solve the issue of long lived decay
- m_PreviousEnergy=fastTrack.GetPrimaryTrack()->GetKineticEnergy();
+ // FIXME: Solve the issue of long lived decay
+ m_PreviousEnergy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
// Check that a decay is possible:
- return m_Decay.AnyAboveThreshold(NPL::ChangeNameFromG4Standard(m_CurrentName),m_ExcitationEnergy);
+ return m_Decay.AnyAboveThreshold(NPL::ChangeNameFromG4Standard(m_CurrentName), m_ExcitationEnergy);
}
////////////////////////////////////////////////////////////////////////////////
-void Decay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
+void Decay::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) {
// Get the track info
const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
G4ThreeVector pdirection = PrimaryTrack->GetMomentum().unit();
@@ -111,12 +106,12 @@ void Decay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
// Randomize within the step
// Assume energy loss is linear within the step
// Assume no scattering
- double rand = G4RandFlat::shoot();
- double length = rand*(m_PreviousLength);
- energy += (1-rand)*(m_PreviousEnergy-energy);
+ double rand = G4RandFlat::shoot();
+ double length = rand * (m_PreviousLength);
+ energy += (1 - rand) * (m_PreviousEnergy - energy);
G4ThreeVector ldir = pdirection;
- ldir*=length;
- localPosition = localPosition - ldir;
+ ldir *= length;
+ localPosition = localPosition - ldir;
//////////////////////////////////////////////////
//////Define the kind of particle to shoot////////
//////////////////////////////////////////////////
@@ -127,54 +122,50 @@ void Decay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
std::vector<double> DPy;
std::vector<double> DPz;
- m_Decay.GenerateEvent(NPL::ChangeNameFromG4Standard(m_CurrentName),m_ExcitationEnergy,energy,
- pdirection.x(),pdirection.y(),pdirection.z(),
- Daughter, Ex,DEK,DPx,DPy,DPz);
+ m_Decay.GenerateEvent(NPL::ChangeNameFromG4Standard(m_CurrentName), m_ExcitationEnergy, energy, pdirection.x(),
+ pdirection.y(), pdirection.z(), Daughter, Ex, DEK, DPx, DPy, DPz);
-
- G4ParticleDefinition* DaughterDef;
+ G4ParticleDefinition* DaughterDef;
unsigned int size = Daughter.size();
- if(size == 0)
+ if (size == 0)
return;
- for(unsigned int i = 0 ; i < size ; i++){
+ for (unsigned int i = 0; i < size; i++) {
// Get the decaying particle
int DaughterZ = Daughter[i].GetZ();
int DaughterA = Daughter[i].GetA();
- DaughterDef=NULL;
+ DaughterDef = NULL;
// neutral particle
- if(DaughterZ==0){
- if(DaughterA==1)
- DaughterDef=G4ParticleTable::GetParticleTable()->FindParticle("neutron");
-
- else if(DaughterA==0){
- DaughterDef=G4ParticleTable::GetParticleTable()->FindParticle("gamma");
- }
+ if (DaughterZ == 0) {
+ if (DaughterA == 1)
+ DaughterDef = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
+ else if (DaughterA == 0) {
+ DaughterDef = G4ParticleTable::GetParticleTable()->FindParticle("gamma");
+ }
}
// proton
- else if (DaughterZ==1 && DaughterA==1 )
- DaughterDef=G4ParticleTable::GetParticleTable()->FindParticle("proton");
+ else if (DaughterZ == 1 && DaughterA == 1)
+ DaughterDef = G4ParticleTable::GetParticleTable()->FindParticle("proton");
// the rest
else
- DaughterDef=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(DaughterZ, DaughterA, Ex[i]);
+ DaughterDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(DaughterZ, DaughterA, Ex[i]);
// Set the momentum direction
- G4ThreeVector Momentum (DPx[i],DPy[i],DPz[i]);
- Momentum=Momentum.unit();
-
- G4DynamicParticle DynamicDaughter(DaughterDef,Momentum,DEK[i]);
+ G4ThreeVector Momentum(DPx[i], DPy[i], DPz[i]);
+ Momentum = Momentum.unit();
+
+ G4DynamicParticle DynamicDaughter(DaughterDef, Momentum, DEK[i]);
fastStep.CreateSecondaryTrack(DynamicDaughter, localPosition, time);
}
- if(size){
+ if (size) {
// Set the end of the step conditions
- fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0,pdirection);
- fastStep.SetPrimaryTrackFinalPosition(worldPosition);
+ fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0, pdirection);
+ fastStep.SetPrimaryTrackFinalPosition(worldPosition);
fastStep.SetTotalEnergyDeposited(0);
- fastStep.SetPrimaryTrackFinalTime (time);
+ fastStep.SetPrimaryTrackFinalTime(time);
fastStep.KillPrimaryTrack();
fastStep.SetPrimaryTrackPathLength(0.0);
- }
-
+ }
}
diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc
index d7a7195437d36312c842f2517a1e22ec2bd1a7a3..761651377df9315ba7914e4f655f3550d75a581b 100644
--- a/NPSimulation/Simulation.cc
+++ b/NPSimulation/Simulation.cc
@@ -1,15 +1,15 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// STL
-#include<cstdlib>
+#include <cstdlib>
-#include "G4RunManager.hh"
#include "G4PhysListFactory.hh"
+#include "G4RunManager.hh"
// UI
#include "G4UImanager.hh"
#include "QBBC.hh"
-#include "G4UIterminal.hh"
#include "G4UItcsh.hh"
+#include "G4UIterminal.hh"
#include "G4VisManager.hh"
#include "G4VisExecutive.hh"
@@ -22,175 +22,171 @@
#include "PrimaryGeneratorAction.hh"
// G4 General Source
-#include "SteppingVerbose.hh"
#include "Randomize.hh"
+#include "SteppingVerbose.hh"
// Root
#include "TRandom.h"
// NPS headers
-#include "SteppingAction.hh"
#include "EventAction.hh"
-#include "RunAction.hh"
#include "NPSimulationVersion.hh"
-//NPL headers
+#include "RunAction.hh"
+#include "SteppingAction.hh"
+// NPL headers
#include "NPOptionManager.h"
#include "RootOutput.h"
-int main(int argc, char** argv){
- // Initialize NPOptionManager object
- NPOptionManager* OptionManager = NPOptionManager::getInstance(argc, argv);
- OptionManager->SetIsSimulation();
- if(OptionManager->GetVerboseLevel() > 0){
- string line;
- line.resize(80,'*');
- cout << endl << line << endl;
- cout << "******************************** NPSimulation ********************************"<< endl;
- cout << line << endl;
- cout << "NPSimulation version: npsimulation-"<< NPS::version_major <<"-" << NPS::version_minor << "-" << NPS::version_dets <<endl;
- cout << " Copyright: NPTool Collaboration "<<endl;
- cout << " Gitlab: https://gitlab.in2p3.fr/np/nptool "<<endl; ;
- cout << line << endl;
- }
-
- // Test if input files are found. If not, exit
- if (OptionManager->IsDefault("EventGenerator"))
- OptionManager->SendErrorAndExit("EventGenerator");
- if (OptionManager->IsDefault("DetectorConfiguration"))
- OptionManager->SendErrorAndExit("DetectorConfiguration");
- // case when input files are here
- G4String EventGeneratorFileName = OptionManager->GetReactionFile();
- G4String DetectorFileName = OptionManager->GetDetectorFile();
-
-
- // initialize the state of the root and geant4 random generator
- if(OptionManager->GetRandomSeed()>0){
- std::cout << " Seeds for random generators set to: " << OptionManager->GetRandomSeed() << std::endl;
- gRandom->SetSeed(OptionManager->GetRandomSeed());
- CLHEP::HepRandom::setTheSeed(OptionManager->GetRandomSeed(),3);
- }
-
-
-
- // my Verbose output class
- G4VSteppingVerbose::SetInstance(new SteppingVerbose);
-
- ///////////////////////////////////////////////////////////////
- ///////////////// Initializing the Root Output ////////////////
- ///////////////////////////////////////////////////////////////
- RootOutput::getInstance(OptionManager->GetOutputFile());
-
- // Construct the default run manager
- G4RunManager* runManager = new G4RunManager;
-
- // set mandatory initialization classes
- DetectorConstruction* detector = new DetectorConstruction();
- runManager->SetUserInitialization(detector);
-
- PhysicsList* physicsList = new PhysicsList();
- physicsList->SetVerboseLevel(0);
- runManager->SetUserInitialization(physicsList);
- PrimaryGeneratorAction* primary = new PrimaryGeneratorAction(detector);
-
- // Initialize Geant4 kernel
- runManager->Initialize();
-
- ///////////////////////////////////////////////////////////////
- /////////// Define UI terminal for interactive mode ///////////
- ///////////////////////////////////////////////////////////////
- // interactive mode : define UI session
- // Get the pointer to the User Interface manager
- G4cout << "//////////// Starting UI ////////////"<< endl;
- G4UImanager* UImanager = G4UImanager::GetUIpointer();
- G4UIExecutive* ui = new G4UIExecutive(argc, argv);
-
-
- ///////////////////////////////////////////////////////////////
- ////////////////////// Reading Reaction ///////////////////////
- ///////////////////////////////////////////////////////////////
- primary->ReadEventGeneratorFile(EventGeneratorFileName);
- runManager->SetUserAction(primary);
-
- ///////////////////////////////////////////////////////////////
- ///////////////// Starting the Stepping Action ////////////////
- ///////////////////////////////////////////////////////////////
- SteppingAction* stepping_action = new SteppingAction() ;
- runManager->SetUserAction(stepping_action) ;
-
- ///////////////////////////////////////////////////////////////
- ////////////////// Starting the Event Action //////////////////
- ///////////////////////////////////////////////////////////////
- EventAction* event_action = new EventAction() ;
- event_action->SetDetector(detector) ;
- runManager->SetUserAction(event_action) ;
-
- ///////////////////////////////////////////////////////////////
- /////////////////// Starting the Run Action ///////////////////
- ///////////////////////////////////////////////////////////////
- RunAction* run_action = new RunAction() ;
- runManager->SetUserAction(run_action);
-
-
- G4VisManager* visManager=NULL;
-
- if(!OptionManager->GetG4BatchMode()){
- string Path_Macro = getenv("NPTOOL");
- Path_Macro+="/NPSimulation/ressources/macro/";
- UImanager->ApplyCommand("/control/execute " +Path_Macro+"verbose.mac");
-
- UImanager->ApplyCommand("/control/execute " +Path_Macro+"aliases.mac");
- visManager = new G4VisExecutive("Quiet");
- visManager->Initialize();
- UImanager->ApplyCommand("/control/execute " +Path_Macro+"vis.mac");
- if (ui->IsGUI()){
- UImanager->ApplyCommand("/control/execute " +Path_Macro+"gui.mac");
- }
-
- #ifdef __APPLE__
- string command= "osascript ";
- command+= getenv("NPTOOL");
- command+="/NPSimulation/ressources/scripts/bringtofront.osa & ";
- int res =system(command.c_str());
- res =0;
- #endif
- }
- else{// if batch mode do not accumulate any track
- UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate 0");
- }
- // Execute user macro
- if(!OptionManager->IsDefault("G4MacroPath")){
- UImanager->ApplyCommand("/control/execute "+ OptionManager->GetG4MacroPath());
+int main(int argc, char** argv) {
+ // Initialize NPOptionManager object
+ NPOptionManager* OptionManager = NPOptionManager::getInstance(argc, argv);
+ OptionManager->SetIsSimulation();
+ if (OptionManager->GetVerboseLevel() > 0) {
+ string line;
+ line.resize(80, '*');
+ cout << endl << line << endl;
+ cout << "******************************** NPSimulation ********************************" << endl;
+ cout << line << endl;
+ cout << "NPSimulation version: npsimulation-" << NPS::version_major << "-" << NPS::version_minor << "-"
+ << NPS::version_dets << endl;
+ cout << " Copyright: NPTool Collaboration " << endl;
+ cout << " Gitlab: https://gitlab.in2p3.fr/np/nptool " << endl;
+ ;
+ cout << line << endl;
+ }
+
+ // Test if input files are found. If not, exit
+ if (OptionManager->IsDefault("EventGenerator"))
+ OptionManager->SendErrorAndExit("EventGenerator");
+ if (OptionManager->IsDefault("DetectorConfiguration"))
+ OptionManager->SendErrorAndExit("DetectorConfiguration");
+ // case when input files are here
+ G4String EventGeneratorFileName = OptionManager->GetReactionFile();
+ G4String DetectorFileName = OptionManager->GetDetectorFile();
+
+ // initialize the state of the root and geant4 random generator
+ if (OptionManager->GetRandomSeed() > 0) {
+ std::cout << " Seeds for random generators set to: " << OptionManager->GetRandomSeed() << std::endl;
+ gRandom->SetSeed(OptionManager->GetRandomSeed());
+ CLHEP::HepRandom::setTheSeed(OptionManager->GetRandomSeed(), 3);
+ }
+
+ // my Verbose output class
+ G4VSteppingVerbose::SetInstance(new SteppingVerbose);
+
+ ///////////////////////////////////////////////////////////////
+ ///////////////// Initializing the Root Output ////////////////
+ ///////////////////////////////////////////////////////////////
+ RootOutput::getInstance(OptionManager->GetOutputFile());
+
+ // Construct the default run manager
+ G4RunManager* runManager = new G4RunManager;
+
+ // set mandatory initialization classes
+ DetectorConstruction* detector = new DetectorConstruction();
+ runManager->SetUserInitialization(detector);
+
+ PhysicsList* physicsList = new PhysicsList();
+ physicsList->SetVerboseLevel(0);
+ runManager->SetUserInitialization(physicsList);
+ PrimaryGeneratorAction* primary = new PrimaryGeneratorAction(detector);
+
+ // Initialize Geant4 kernel
+ runManager->Initialize();
+
+ ///////////////////////////////////////////////////////////////
+ /////////// Define UI terminal for interactive mode ///////////
+ ///////////////////////////////////////////////////////////////
+ // interactive mode : define UI session
+ // Get the pointer to the User Interface manager
+ G4cout << "//////////// Starting UI ////////////" << endl;
+ G4UImanager* UImanager = G4UImanager::GetUIpointer();
+ G4UIExecutive* ui = new G4UIExecutive(argc, argv);
+
+ ///////////////////////////////////////////////////////////////
+ ////////////////////// Reading Reaction ///////////////////////
+ ///////////////////////////////////////////////////////////////
+ primary->ReadEventGeneratorFile(EventGeneratorFileName);
+ runManager->SetUserAction(primary);
+
+ ///////////////////////////////////////////////////////////////
+ ///////////////// Starting the Stepping Action ////////////////
+ ///////////////////////////////////////////////////////////////
+ SteppingAction* stepping_action = new SteppingAction();
+ runManager->SetUserAction(stepping_action);
+
+ ///////////////////////////////////////////////////////////////
+ ////////////////// Starting the Event Action //////////////////
+ ///////////////////////////////////////////////////////////////
+ EventAction* event_action = new EventAction();
+ event_action->SetDetector(detector);
+ runManager->SetUserAction(event_action);
+
+ ///////////////////////////////////////////////////////////////
+ /////////////////// Starting the Run Action ///////////////////
+ ///////////////////////////////////////////////////////////////
+ RunAction* run_action = new RunAction();
+ runManager->SetUserAction(run_action);
+
+ G4VisManager* visManager = NULL;
+
+ if (!OptionManager->GetG4BatchMode()) {
+ string Path_Macro = getenv("NPTOOL");
+ Path_Macro += "/NPSimulation/ressources/macro/";
+ UImanager->ApplyCommand("/control/execute " + Path_Macro + "verbose.mac");
+
+ UImanager->ApplyCommand("/control/execute " + Path_Macro + "aliases.mac");
+ visManager = new G4VisExecutive("Quiet");
+ visManager->Initialize();
+ UImanager->ApplyCommand("/control/execute " + Path_Macro + "vis.mac");
+ if (ui->IsGUI()) {
+ UImanager->ApplyCommand("/control/execute " + Path_Macro + "gui.mac");
}
-
- // Start the session
- if(!OptionManager->GetG4BatchMode())
- ui->SessionStart();
-
-
- delete ui;
-
- if(visManager)
- delete visManager;
-
- ///////////////////////////////////////////////////////////////
- ////////////////////// Job termination ////////////////////////
- ///////////////////////////////////////////////////////////////
- // Save the Geant4 random generator internal generator state in a TASCII
- // file store with the root output
- std::ofstream file(".geant4_random_state");
- CLHEP::HepRandom::saveFullState(file);
- file.close();
- TAsciiFile* aFile = new TAsciiFile();
- aFile->SetNameTitle("G4RandomFinalState",".geant4_random_state");
- aFile->Append(".geant4_random_state");
- RootOutput::getInstance()->GetFile()->cd();
- aFile->Write(0,TAsciiFile::kOverwrite);
- int dummy = system("rm .geant4_random_state");
- dummy*=2;
- // delete primary; delete detector;
-
- delete runManager;
-
- RootOutput::getInstance()->Destroy();
- return 0;
+
+#ifdef __APPLE__
+ string command = "osascript ";
+ command += getenv("NPTOOL");
+ command += "/NPSimulation/ressources/scripts/bringtofront.osa & ";
+ int res = system(command.c_str());
+ res = 0;
+#endif
+ }
+ else { // if batch mode do not accumulate any track
+ UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate 0");
+ }
+ // Execute user macro
+ if (!OptionManager->IsDefault("G4MacroPath")) {
+ UImanager->ApplyCommand("/control/execute " + OptionManager->GetG4MacroPath());
+ }
+
+ // Start the session
+ if (!OptionManager->GetG4BatchMode())
+ ui->SessionStart();
+
+ delete ui;
+
+ if (visManager)
+ delete visManager;
+
+ ///////////////////////////////////////////////////////////////
+ ////////////////////// Job termination ////////////////////////
+ ///////////////////////////////////////////////////////////////
+ // Save the Geant4 random generator internal generator state in a TASCII
+ // file store with the root output
+ std::ofstream file(".geant4_random_state");
+ CLHEP::HepRandom::saveFullState(file);
+ file.close();
+ TAsciiFile* aFile = new TAsciiFile();
+ aFile->SetNameTitle("G4RandomFinalState", ".geant4_random_state");
+ aFile->Append(".geant4_random_state");
+ RootOutput::getInstance()->GetFile()->cd();
+ aFile->Write(0, TAsciiFile::kOverwrite);
+ int dummy = system("rm .geant4_random_state");
+ dummy *= 2;
+ // delete primary; delete detector;
+
+ delete runManager;
+
+ RootOutput::getInstance()->Destroy();
+ return 0;
}
diff --git a/Projects/E805/Analysis.cxx b/Projects/E805/Analysis.cxx
index db34beabf4293a75347e9eae6b2b1c801722c178..1dde6158af8f09f8054044c857b9781f2196139b 100755
--- a/Projects/E805/Analysis.cxx
+++ b/Projects/E805/Analysis.cxx
@@ -53,17 +53,21 @@ void Analysis::Init(){
//
ProtonSi = NPL::EnergyLoss(Path+ "proton_Si.G4table", "G4Table", 100);
+
+ Cal = CalibrationManager::getInstance();
}
///////////////////////////// Initialize some important parameters //////////////////////////////////
bool Analysis::UnallocateBeforeBuild(){
- //std::cout << "test unallocate" << std::endl;
- GATCONFMASTER = **GATCONFMASTER_;
- return (GATCONFMASTER > 0);
- //DATATRIG_CATS = **DATATRIG_CATS_;
- //return (DATATRIG_CATS > 0);
- //return true;
+ // std::cout << "test unallocate" << std::endl;
+ //GATCONFMASTER = **GATCONFMASTER_;
+ //return (GATCONFMASTER > 0);
+ return true;
+}
+
+bool Analysis::UnallocateBeforeTreat(){
+ return true;
}
bool Analysis::FillOutputCondition(){
@@ -73,8 +77,86 @@ bool Analysis::FillOutputCondition(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
-
+
+ // if(M2->CsI_E.size() > 0)
+ // std::cout << "Analysis test " << M2->CsI_E[0] << " " << M2->CsI_E.size() << " " << "\n \n";
+
ReInit();
+ //////////////////// MUST2 Part ////////////////////
+ int M2_size = M2->Si_E.size();
+ for(unsigned int countMust2 = 0 ; countMust2 < M2_size ; countMust2++){
+ M2_TelescopeM++;
+ // MUST2
+ int TelescopeNumber = M2->TelescopeNumber[countMust2];
+
+ // Part 1 : Impact Angle
+ ThetaM2Surface = 0;
+ ThetaNormalTarget = 0;
+
+ //BeamImpact = TVector3(Xf,Yf,0);
+ BeamImpact = TVector3(0,0,0);
+
+ //BeamDirection = TVector3(CATS2_X - CATS1_X,CATS2_Y - CATS1_Y,497);
+ BeamDirection = TVector3(0,0,1);
+
+ TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
+ M2_ThetaLab.push_back(HitDirection.Angle( BeamDirection ));
+
+ //std::cout << M2 -> GetPositionOfInteraction(countMust2).X() << " " << M2 -> GetPositionOfInteraction(countMust2).Y() << " "<< M2 -> GetPositionOfInteraction(countMust2).Z() << std::endl;
+ M2_X.push_back(M2 -> GetPositionOfInteraction(countMust2).X());
+ M2_Y.push_back(M2 -> GetPositionOfInteraction(countMust2).Y());
+ M2_Z.push_back(M2 -> GetPositionOfInteraction(countMust2).Z());
+
+ ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+
+ // Part 2 : Impact Energy
+ Si_E_M2 = 0;
+ CsI_E_M2 = 0;
+ int CristalNb = M2->CsI_N[countMust2];
+
+ Si_E_M2 = M2->Si_E[countMust2];
+ CsI_E_M2= M2->CsI_E[countMust2];
+ if(CsI_E_M2 > 0)
+ std::cout << "Analysis" << CsI_E_M2 << " " << M2->CsI_E.size() << " " << CristalNb << "\n \n";
+ for(unsigned int i = 0; i < ParticleType.size(); i++){
+ Energy[ParticleType[i]] = 0;
+
+ if(CsI_E_M2>8192 ){
+ // The energy in CsI is calculate form dE/dx Table because
+ Energy[ParticleType[i]] = Cal->ApplyCalibration(ParticleType[i]+"MUST2/T"+NPL::itoa(TelescopeNumber)+"_CsI"+NPL::itoa(CristalNb)+"_E",CsI_E_M2);
+ std::cout << Energy[ParticleType[i]] << "\n";
+ //Energy = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
+ //Energy+=Si_E_M2;
+ }
+
+// else
+ Energy[ParticleType[i]] += Si_E_M2;
+
+
+ }
+ // Evaluate energy using the thickness
+ //Energy = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
+ // Target Correction
+ //M2_ELab.push_back(LightTarget.EvaluateInitialEnergy(Energy ,TargetThickness*0.5, ThetaNormalTarget));
+ M2_ELab.push_back(Energy["proton"]);
+ M2_dE.push_back(Si_E_M2);
+
+ // Part 3 : Excitation Energy Calculation
+ M2_Ex_p.push_back(reaction->ReconstructRelativistic( Energy["proton"] , M2_ThetaLab[countMust2] ));
+ M2_Ex_d.push_back(reaction->ReconstructRelativistic( Energy["deuteron"] , M2_ThetaLab[countMust2] ));
+ M2_Ex_t.push_back(reaction->ReconstructRelativistic( Energy["triton"] , M2_ThetaLab[countMust2] ));
+ M2_Ex_a.push_back(reaction->ReconstructRelativistic( Energy["alpha"] , M2_ThetaLab[countMust2] ));
+ M2_ThetaLab[countMust2]=M2_ThetaLab[countMust2]/deg;
+
+ // Part 4 : Theta CM Calculation
+ M2_ThetaCM.push_back(reaction->EnergyLabToThetaCM( M2_ELab[countMust2] , M2_ThetaLab[countMust2])/deg);
+
+ //if(proton_cut[TelescopeNumber]){
+ // std::cout << "Test :" << proton_cut[TelescopeNumber] << " \n";
+ //
+ }//end loop MUST2
+
//for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
//Si_E_M2 = M2->Si_E[countMust2];
@@ -101,9 +183,24 @@ void Analysis::TreatEvent(){
void Analysis::InitOutputBranch() {
+ RootOutput::getInstance()->GetTree()->Branch("M2_TelescopeM",&M2_TelescopeM,"M2_TelescopeM/s");
+ RootOutput::getInstance()->GetTree()->Branch("M2_Ex_p",&M2_Ex_p);
+ RootOutput::getInstance()->GetTree()->Branch("M2_Ex_d",&M2_Ex_d);
+ RootOutput::getInstance()->GetTree()->Branch("M2_Ex_t",&M2_Ex_t);
+ RootOutput::getInstance()->GetTree()->Branch("M2_Ex_a",&M2_Ex_a);
+ //RootOutput::getInstance()->GetTree()->Branch("ExNoBeam",&ExNoBeam,"ExNoBeam/D");
+ //RootOutput::getInstance()->GetTree()->Branch("ExNoProton",&ExNoProton,"ExNoProton/D");
+ RootOutput::getInstance()->GetTree()->Branch("M2_ELab",&M2_ELab);
+ RootOutput::getInstance()->GetTree()->Branch("M2_ThetaLab",&M2_ThetaLab);
+ RootOutput::getInstance()->GetTree()->Branch("M2_ThetaCM",&M2_ThetaCM);
+ RootOutput::getInstance()->GetTree()->Branch("M2_X",&M2_X);
+ RootOutput::getInstance()->GetTree()->Branch("M2_Y",&M2_Y);
+ RootOutput::getInstance()->GetTree()->Branch("M2_Z",&M2_Z);
+ RootOutput::getInstance()->GetTree()->Branch("M2_dE",&M2_dE);
+ RootOutput::getInstance()->GetTree()->Branch("CsI_E_M2",&CsI_E_M2);
RootOutput::getInstance()->GetTree()->Branch("M2_ECsI_from_deltaE",&M2_ECsI_from_deltaE);
- RootOutput::getInstance()->GetTree()->Branch("Beta_from_deltaE",&Beta_from_deltaE);
- RootOutput::getInstance()->GetTree()->Branch("Beth_from_deltaE",&Beth_from_deltaE);
+ RootOutput::getInstance()->GetTree()->Branch("GATCONF",&GATCONFMASTER);
+
}
void Analysis::UnallocateVariables(){
@@ -111,15 +208,32 @@ void Analysis::UnallocateVariables(){
void Analysis::InitInputBranch(){
TTreeReader* inputTreeReader = RootInput::getInstance()->GetTreeReader();
- GATCONFMASTER_ = new TTreeReaderValue<unsigned short>(*inputTreeReader,"GATCONFMASTER");
+ // GATCONFMASTER_ = new TTreeReaderValue<unsigned short>(*inputTreeReader,"GATCONFMASTER");
//DATATRIG_CATS_ = new TTreeReaderValue<unsigned short>(*inputTreeReader,"DATATRIG_CATS");
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInit(){
- //M2_ECsI_from_deltaE.clear();
- //Beta_from_deltaE.clear();
- //Beth_from_deltaE.clear();
+
+ Theta_seg = -1000;
+ Phi_seg = -1000;
+ M2_TelescopeM= 0;
+ M2_Ex_p.clear();
+ M2_Ex_d.clear();
+ M2_Ex_t.clear();
+ M2_Ex_a.clear();
+ M2_ECsI_from_deltaE.clear();
+ //ExNoBeam=ExNoProto.clear();
+ //EDC.clear();
+ M2_ELab.clear();
+ //BeamEnergy .clear();
+ M2_ThetaLab.clear();
+ M2_ThetaCM.clear();
+ M2_X.clear();
+ M2_Y.clear();
+ M2_Z.clear();
+ M2_dE.clear();
+
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/E805/Analysis.h b/Projects/E805/Analysis.h
index 4ad3a6d6289db1ee3849f935fc18eed262fc4a32..f12284ee3d45e6e51182f1daac83ca61b23f92d5 100755
--- a/Projects/E805/Analysis.h
+++ b/Projects/E805/Analysis.h
@@ -50,7 +50,7 @@ class Analysis: public NPL::VAnalysis{
void Init();
bool UnallocateBeforeBuild();
bool FillOutputCondition();
- bool UnallocateBeforeTreat(){return true;};
+ bool UnallocateBeforeTreat();
void TreatEvent();
void End();
void ReInit();
@@ -87,7 +87,10 @@ class Analysis: public NPL::VAnalysis{
// MUST2 info
unsigned short M2_TelescopeM;
- std::vector<double> M2_Ex;
+ std::vector<double> M2_Ex_p;
+ std::vector<double> M2_Ex_d;
+ std::vector<double> M2_Ex_t;
+ std::vector<double> M2_Ex_a;
std::vector<double> M2_ExNoBeam;
std::vector<double> M2_ExNoProton;
std::vector<double> M2_EDC;
@@ -162,7 +165,7 @@ class Analysis: public NPL::VAnalysis{
float EnergyAddBackDoppler;
float EnergyAddBack;
int ExogamDetNb[3];
- int CristalNb[3];
+ // int CristalNb[3];
int SegmentNb[3];
std::vector<int> event1;
@@ -241,7 +244,8 @@ class Analysis: public NPL::VAnalysis{
double ThetaMGSurface ;
double Si_E_M2 ;
double CsI_E_M2 ;
- double Energy ;
+ std::vector<string> ParticleType{"proton","deuteron","triton","alpha"};
+ std::map<TString, double> Energy ;
double BeamEnergy;
double ThetaGDSurface ;
@@ -285,7 +289,7 @@ class Analysis: public NPL::VAnalysis{
double ZDD_R = 1027*mm;
private:
- TMust2Data* M2_Raw;
+ // TMust2Data* M2_Raw;
TMust2Physics* M2;
TCATSPhysics* CATS;
TZDDPhysics* ZDD;
@@ -354,6 +358,7 @@ class Analysis: public NPL::VAnalysis{
default_random_engine generator;
normal_distribution<double> distribution = normal_distribution<double>(0.0,1.0);
+ CalibrationManager* Cal;
};
#endif
\ No newline at end of file
diff --git a/Projects/E805/CalibrationAphaTest.txt b/Projects/E805/CalibrationAphaTest.txt
index e152e0f0fce629547b7fb4c0e482119aa0e820a8..bd1bdebf7b86abfc232b67c19966d43f5b6d4442 100644
--- a/Projects/E805/CalibrationAphaTest.txt
+++ b/Projects/E805/CalibrationAphaTest.txt
@@ -8,5 +8,29 @@ CalibrationFilePath
./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM2.cal
./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM3.cal
./data/NPRoot/Calibration/CalAlpha_r0388_06+/Cal_Str_Y_E_MM4.cal
-
-
\ No newline at end of file
+
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM1_r0183.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM1_r0184.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM2_r0183.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM2_r0184.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM3_r0189.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM3_r0190.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_X_T_MM4_r0187.cal
+ ./data/NPRoot/Calibration/Time/Coeff/Cal_Str_Y_T_MM4_r0188_manual.cal
+
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_proton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_deuteron.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_triton.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM1_alpha.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM2_alpha.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM3_alpha.txt
+ ./data/NPRoot/Calibration/Run_CSICalib/CsI_MM4_alpha.txt
diff --git a/Projects/E805/RunToTreat.txt b/Projects/E805/RunToTreat.txt
index 48d7cb36b564021d576a138ddb7c57bc0c093022..ef2289f4dcda12c6627dd33451c491f502e5eb28 100644
--- a/Projects/E805/RunToTreat.txt
+++ b/Projects/E805/RunToTreat.txt
@@ -1,4 +1,5 @@
TTreeName
- RD
+ SimulatedTree
RootFileName
- ./RootR/r0316_000r.root
+ ./data/NPRoot/Simulation/test.root
+
diff --git a/Projects/E805/project.config b/Projects/E805/project.config
index e7d7ae93f84611521fd60633b734e74b87245043..f4d49abf59994412aa11998b16b6959b5abd9a3e 100644
--- a/Projects/E805/project.config
+++ b/Projects/E805/project.config
@@ -3,5 +3,6 @@ Project E805
SimulationOutput= ./data/NPRoot/Simulation/
EnergyLoss= ../../Inputs/EnergyLoss/
CalibrationOutput= ./data/NPRoot/Calibration/
+ Cuts= ./data/NPRoot/Cuts/
diff --git a/Projects/Gaspard/132Sndp.reaction b/Projects/Gaspard/132Sndp.reaction
index dc79bccafb3d933f5c61c24fac20950b03d30ac3..a8a4e8403f1c744e8f63c14272b0fb6576de654f 100644
--- a/Projects/Gaspard/132Sndp.reaction
+++ b/Projects/Gaspard/132Sndp.reaction
@@ -6,11 +6,11 @@ Beam
Particle= 132Sn
Energy= 1320
SigmaEnergy= 0 MeV
- SigmaThetaX= 0.42 MeV
+ SigmaThetaX= 0.42 deg
SigmaPhiY= 0.42 deg
+ MeanThetaX= 0.3 deg
SigmaX= 1.7 mm
SigmaY= 1.7 mm
- MeanThetaX= 0.3 deg
MeanPhiY= 0 deg
MeanX= 0.0 mm
MeanY= 0.0 mm
diff --git a/Projects/Strasse/geometry/strasse_CAD.detector b/Projects/Strasse/geometry/strasse_CAD.detector
index 8e917df8dce733c15ca9602739d556a3e86b85e2..56acc97eda84301e75ea19b4ea78b2dd6cfd2549 100644
--- a/Projects/Strasse/geometry/strasse_CAD.detector
+++ b/Projects/Strasse/geometry/strasse_CAD.detector
@@ -88,10 +88,10 @@ Strasse Outer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Strasse InactiveMaterial
- Chamber= ./geometry/STRASSE_Chamber.stl
- Stars= ./geometry/STRASSE_StarSupports.stl
- Base= ./geometry/STRASSE_Base.stl
- Blades= ./geometry/STRASSE_Blades.stl
+ %Chamber= ./geometry/STRASSE_Chamber.stl
+ %Stars= ./geometry/STRASSE_StarSupports.stl
+ %Base= ./geometry/STRASSE_Base.stl
+ %Blades= ./geometry/STRASSE_Blades.stl
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
Catana CSV
Path= geometry/Catana.csv
diff --git a/Projects/Strasse/reaction/Ca54.source b/Projects/Strasse/reaction/Ca54.source
index d5bfa4982e64b2902a23f00caed141bce71117df..cd3283d5049d4bfaaace39c7c25243a1e090c136 100755
--- a/Projects/Strasse/reaction/Ca54.source
+++ b/Projects/Strasse/reaction/Ca54.source
@@ -11,7 +11,7 @@ Beam
MeanPhiY= 0 deg
MeanX= 0 mm
MeanY= 0 mm
- ExcitationEnergy= 3680 keV
+ ExcitationEnergy= 2466 keV
ZEmmission= 0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LevelData 54Ca
diff --git a/Projects/Strasse/reaction/Ca54p2p.reaction b/Projects/Strasse/reaction/Ca54p2p.reaction
index 586513a6bd731fa3aa19ee10306887d4e7c5a89d..53851ab68c2f11ff451dcc6c163a135cea702b8f 100755
--- a/Projects/Strasse/reaction/Ca54p2p.reaction
+++ b/Projects/Strasse/reaction/Ca54p2p.reaction
@@ -23,7 +23,7 @@ QFSReaction
ExcitationEnergyHeavy= 3.6800 MeV
MomentumSigma= 50.0
ShootHeavy= 1
- ShootLight= 1
+ ShootLight= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LevelData 54Ca