diff --git a/NPLib/CMakeLists.txt b/NPLib/CMakeLists.txt
index dceaa08a0feade1d72f9625edf8ab45af486b964..dae59391ea14715418ae5877ced3730b23deaf82 100644
--- a/NPLib/CMakeLists.txt
+++ b/NPLib/CMakeLists.txt
@@ -75,10 +75,11 @@ include_directories("Core/")
include_directories("Physics/")
include_directories("TrackReconstruction/")
include_directories("Online/")
+include_directories("Neutron/")
# Call the macro
subdirlist(SUB_DIRECTORY ${CMAKE_BINARY_DIR})
-set(SUB_DIRECTORY ${SUB_DIRECTORY} Core Physics TrackReconstruction Calibration Online Utility)
+set(SUB_DIRECTORY ${SUB_DIRECTORY} Core Physics TrackReconstruction Calibration Online Utility Neutron)
# Add each sub folder to the project
set(TARGET_LIST "")
foreach(subdir ${SUB_DIRECTORY})
diff --git a/NPLib/Detectors/AnnularS1/TAnnularS1Physics.cxx b/NPLib/Detectors/AnnularS1/TAnnularS1Physics.cxx
index 6440f49750e4ca17083915ccf8fbb8d4a74999f4..048a9e73641c9b123c152076812821deebbbafd2 100644
--- a/NPLib/Detectors/AnnularS1/TAnnularS1Physics.cxx
+++ b/NPLib/Detectors/AnnularS1/TAnnularS1Physics.cxx
@@ -40,6 +40,10 @@ using namespace NPUNITS;
#include "TChain.h"
///////////////////////////////////////////////////////////////////////////
+// Global
+TRandom *Rand = new TRandom3();
+
+
ClassImp(TAnnularS1Physics)
///////////////////////////////////////////////////////////////////////////
TAnnularS1Physics::TAnnularS1Physics(){
@@ -150,7 +154,7 @@ void TAnnularS1Physics::PreTreat(){
}
// Sector T
- unsigned int sizeSectorT = m_EventData->GetS1ThetaTMult();
+ unsigned int sizeSectorT = m_EventData->GetS1PhiTMult();
for(unsigned int i = 0 ; i < sizeSectorT ; ++i){
m_PreTreatedData->SetS1PhiTDetectorNbr( m_EventData->GetS1PhiTDetectorNbr(i) );
m_PreTreatedData->SetS1PhiTStripNbr( m_EventData->GetS1PhiTStripNbr(i) );
@@ -239,7 +243,7 @@ void TAnnularS1Physics::ReadAnalysisConfig(){
getline(AnalysisConfigFile, LineBuffer);
// search for "header"
- if (LineBuffer.compare(0, 11, "ConfigAnnularS1") == 0) ReadingStatus = true;
+ if (LineBuffer.compare(0, 15, "ConfigAnnularS1") == 0) ReadingStatus = true;
// loop on tokens and data
while (ReadingStatus ) {
@@ -255,7 +259,7 @@ void TAnnularS1Physics::ReadAnalysisConfig(){
else if (whatToDo=="MAX_STRIP_MULTIPLICITY") {
AnalysisConfigFile >> DataBuffer;
m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str() );
- cout << "MAXIMUN STRIP MULTIPLICITY " << m_MaximumStripMultiplicityAllowed << endl;
+ cout << "MAXIMUM STRIP MULTIPLICITY " << m_MaximumStripMultiplicityAllowed << endl;
}
else if (whatToDo=="STRIP_ENERGY_MATCHING_SIGMA") {
@@ -535,6 +539,32 @@ TVector3 TAnnularS1Physics::GetPositionOfInteraction(const int i) const{
return(Position) ;
}
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+TVector3 TAnnularS1Physics::GetRandomisedPositionOfInteraction(const int i) const{
+
+ double R_Min = 24;
+ double R_Max = 48;
+
+ double Phi_Min = 0 ;
+ double Phi_Max = 360;
+
+ int NumberofRing = 16 ; //Per Quadrant
+ int NumberofSector = 16 ; //Per detector, ( 4 in each Quad)
+
+ double rho_half_range = 0.5*(R_Max-R_Min)/NumberofRing ; // ring strip spacing in mm
+ double phi_half_range = 0.5*(Phi_Max-Phi_Min)/NumberofSector ; //radial strip spacing in rad
+ TVector3 OriginalPosition = GetPositionOfInteraction(i);
+ double rho = OriginalPosition.Perp();
+ double phi = OriginalPosition.Phi();
+ double z = OriginalPosition.Z();
+ // randomises within a given detector ring and sector
+ double rho_min2 = (rho-rho_half_range)*(rho-rho_half_range) ; // ^2 to reproduce a randomization in the arc
+ double rho_max2 = (rho+rho_half_range)*(rho+rho_half_range) ;
+ double rho_rand = sqrt(Rand->Uniform(rho_min2,rho_max2));// sqrt is necessary for realistic randomise!
+ double phi_rand = phi + Rand->Uniform(-phi_half_range, +phi_half_range)*M_PI/180;
+
+ return( TVector3(rho_rand*cos(phi_rand),rho_rand*sin(phi_rand),z) ) ;
+}
/////////////////////////////
void TAnnularS1Physics::InitializeStandardParameter(){
// Enable all channel
diff --git a/NPLib/Detectors/AnnularS1/TAnnularS1Physics.h b/NPLib/Detectors/AnnularS1/TAnnularS1Physics.h
index e0a96bb071c93ea8e9f2829653e0e9931388ed39..ebb49c5dbde4b1827639aa726ab5656ed95d9eac 100644
--- a/NPLib/Detectors/AnnularS1/TAnnularS1Physics.h
+++ b/NPLib/Detectors/AnnularS1/TAnnularS1Physics.h
@@ -34,6 +34,7 @@
#include "TVector2.h"
#include "TVector3.h"
#include "TObject.h"
+#include "TRandom3.h"
using namespace std ;
@@ -146,6 +147,7 @@ class TAnnularS1Physics : public TObject, public NPL::VDetector{
int GetEventMultiplicity() const { return EventMultiplicity; };
TVector3 GetPositionOfInteraction(const int i) const;
+ TVector3 GetRandomisedPositionOfInteraction(const int i) const;
TVector3 GetDetectorNormal(const int i) const;
private: // Parameter used in the analysis
diff --git a/NPLib/Detectors/MUST2/TMust2Physics.cxx b/NPLib/Detectors/MUST2/TMust2Physics.cxx
index 71cb70555d951407e9f41a3b054edb98ad7ae91e..10435fc52d3bcdb69a5f1132c4da95e69b557860 100644
--- a/NPLib/Detectors/MUST2/TMust2Physics.cxx
+++ b/NPLib/Detectors/MUST2/TMust2Physics.cxx
@@ -60,10 +60,10 @@ ClassImp(TMust2Physics)
m_SiLi_E_RAW_Threshold = 8200;
m_CsI_E_RAW_Threshold = 8200;
// Calibrated Threshold
- m_Si_X_E_Threshold = 0;
- m_Si_Y_E_Threshold = 0;
- m_SiLi_E_Threshold = 0;
- m_CsI_E_Threshold = 0;
+ m_Si_X_E_Threshold = 0*keV;
+ m_Si_Y_E_Threshold = 0*keV;
+ m_SiLi_E_Threshold = 0*keV;
+ m_CsI_E_Threshold = 0*keV;
m_Ignore_not_matching_SiLi = false;
m_Ignore_not_matching_CsI = false;
@@ -80,7 +80,7 @@ ClassImp(TMust2Physics)
m_SiLi_MatchingY.resize(16, 0);
m_SiLi_MatchingX[0] = 112;
- m_SiLi_MatchingY[1] = 112;
+ m_SiLi_MatchingY[0] = 112;
m_SiLi_MatchingX[1] = 112;
m_SiLi_MatchingY[1] = 80;
@@ -674,7 +674,7 @@ void TMust2Physics::ReadAnalysisConfig() {
else if (whatToDo == "MAX_STRIP_MULTIPLICITY") {
AnalysisConfigFile >> DataBuffer;
m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str());
- cout << "MAXIMUN STRIP MULTIPLICITY "
+ cout << "MAXIMUM STRIP MULTIPLICITY "
<< m_MaximumStripMultiplicityAllowed << endl;
}
@@ -803,31 +803,30 @@ void TMust2Physics::ReadAnalysisConfig() {
else if (whatToDo == "SI_X_E_THRESHOLD") {
AnalysisConfigFile >> DataBuffer;
- m_Si_X_E_Threshold = atof(DataBuffer.c_str());
- cout << whatToDo << " " << m_Si_X_E_Threshold << endl;
+ m_Si_X_E_Threshold = atof(DataBuffer.c_str()) *keV ;
+ cout << whatToDo << " " << m_Si_X_E_Threshold << " MeV" << endl;
}
else if (whatToDo == "SI_Y_E_THRESHOLD") {
AnalysisConfigFile >> DataBuffer;
- m_Si_Y_E_Threshold = atof(DataBuffer.c_str());
- cout << whatToDo << " " << m_Si_Y_E_Threshold << endl;
+ m_Si_Y_E_Threshold = atof(DataBuffer.c_str()) *keV;
+ cout << whatToDo << " " << m_Si_Y_E_Threshold << " MeV" << endl;
}
else if (whatToDo == "SILI_E_THRESHOLD") {
AnalysisConfigFile >> DataBuffer;
- m_SiLi_E_Threshold = atof(DataBuffer.c_str());
- cout << whatToDo << " " << m_SiLi_E_Threshold << endl;
+ m_SiLi_E_Threshold = atof(DataBuffer.c_str()) *keV;
+ cout << whatToDo << " " << m_SiLi_E_Threshold << " MeV" << endl;
}
else if (whatToDo == "CSI_E_THRESHOLD") {
AnalysisConfigFile >> DataBuffer;
- m_CsI_E_Threshold = atof(DataBuffer.c_str());
- cout << whatToDo << " " << m_CsI_E_Threshold << endl;
+ m_CsI_E_Threshold = atof(DataBuffer.c_str()) *keV;
+ cout << whatToDo << " " << m_CsI_E_Threshold << " MeV" << endl;
}
- else {
- ReadingStatus = false;
- }
+ else if (AnalysisConfigFile.eof()) ReadingStatus = false;
+
}
}
}
@@ -1398,7 +1397,6 @@ double fSiLi_E(const TMust2Data* m_EventData, const int& i) {
name += "_SiLi";
name += NPL::itoa(m_EventData->GetMMSiLiEPadNbr(i));
name += "_E";
-
return CalibrationManager::getInstance()->ApplyCalibration(
name, m_EventData->GetMMSiLiEEnergy(i),1);
}
diff --git a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h
index 790f89e0d5579a97c65ab28215e70c935d0ce4f8..7ee6615c139dcc921f306026bd0dc30ba095cfd7 100644
--- a/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h
+++ b/NPLib/Detectors/Samurai/TSamuraiFDC0Physics.h
@@ -171,6 +171,7 @@ class TSamuraiFDC0Physics : public TObject, public NPL::VDetector{
double GetPosX(){return PosX;}
double GetPosY(){return PosY;}
double GetThetaX(){return ThetaX;}
+ double GetPhiY(){return PhiY;}
double GetDevX(){return devX;}
double GetDevY(){return devY;}
int GetPileUp(){return PileUp;}
diff --git a/NPLib/Neutron/CMakeLists.txt b/NPLib/Neutron/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..cf61d53dd863e07c1e0c240d1911ab08ff59dc49
--- /dev/null
+++ b/NPLib/Neutron/CMakeLists.txt
@@ -0,0 +1,7 @@
+# add_custom_command(OUTPUT CrossTalkDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh CrossTalk.h CrossTalkDict.cxx CrossTalk.rootmap libNPNeutron.so NPNeutronLinkDef.h DEPENDS CrossTalk.h)
+
+add_library(NPNeutron SHARED NPCrossTalk.cxx)
+
+target_link_libraries(NPNeutron ${ROOT_LIBRARIES} NPCore)
+
+install(FILES NPCrossTalk.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Neutron/NPCrossTalk.cxx b/NPLib/Neutron/NPCrossTalk.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f0aedf6baf59e4ce9c29cb6572b7ac19594245ed
--- /dev/null
+++ b/NPLib/Neutron/NPCrossTalk.cxx
@@ -0,0 +1,170 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2020 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 : Cyril Lenain contact address: lenain@lpcaen.in2p3.fr *
+ * *
+ * Creation Date : November 2020 *
+ * Last update : November 2020 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class find multi-neutrons in an neutron detection array by *
+ * rejecting "crosstalk" *
+ *****************************************************************************/
+
+#include"NPCrossTalk.h"
+#include "Math/Factory.h"
+#include "TError.h"
+#include "TGraph.h"
+
+#include <iostream>
+using namespace std;
+using namespace NPL;
+
+////////////////////////////////////////////////////////////////////////////////
+CrossTalk::CrossTalk(){
+ // Radius of a cluster = Coef * dXdYdZ
+ coef = 3;
+ // this avoid error
+ gErrorIgnoreLevel = kError;
+ /* FirstHit = -1; */
+
+}
+////////////////////////////////////////////////////////////////////////////////
+CrossTalk::~CrossTalk(){
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+void CrossTalk::AddHitVector(const vector<double>& X, const vector<double>& Y, const vector<double>& Z,const vector<double>& dX, const vector<double>& dY, const vector<double>& dZ, const vector<double> &T){
+
+ HitX = &X;
+ HitY = &Y;
+ HitZ = &Z;
+
+ HitdX = &dX;
+ HitdY = &dY;
+ HitdZ = &dZ;
+
+ HitT = &T;
+ sizeHit = X.size();
+}
+
+vector<int> CrossTalk::ComputeCrossTalk(){
+
+ FirstHit = -1;
+
+ static double x1,y1,z1,dx1,dy1,dz1,t1;
+ static double x2,y2,z2,dx2,dy2,dz2,t2;
+ static double Dist, dR1, dR2;
+
+ // A different Cluster number (starting at 1) is assigned to each hit
+ static vector<int> ID_ClustHit;
+ ID_ClustHit.clear();
+ for(int i =0 ; i < sizeHit; i++){
+ ID_ClustHit.push_back(i+1);
+ }
+
+ //Test each Dist(n-n) to find clusters
+ //When 2 neutrons are part of a the same cluster, change their ID_ClustHit for the lowest
+ for(int j = 0; j < sizeHit-1; j++){
+ x1 = (*HitX)[j], y1 = (*HitY)[j], z1 = (*HitZ)[j], dx1 = (*HitdX)[j], dy1 = (*HitdY)[j], dz1 = (*HitdZ)[j], t1 = (*HitT)[j];
+ dR1 = sqrt( dx1*dx1 + dy1*dy1 + dz1*dz1);
+ for(int jj = j+1; jj < sizeHit ; jj++){
+ x2 = (*HitX)[jj], y2 = (*HitY)[jj], z2 = (*HitZ)[jj];
+ Dist = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
+ if(Dist < coef*dR1){
+ if(ID_ClustHit[jj] < ID_ClustHit[j]){
+ ID_ClustHit[j] = ID_ClustHit[jj];
+ }
+ else if(ID_ClustHit[jj] > ID_ClustHit[j]){
+ ID_ClustHit[jj] = ID_ClustHit[j];
+ }
+ }
+ }
+ }
+ //Create a map to hold the Hit_ID for each cluster
+ static map<unsigned int, vector<unsigned int>> mapOfClust;
+ mapOfClust.clear();
+ for(unsigned int i = 0; i < sizeHit; i++){
+ mapOfClust[ID_ClustHit[i]].push_back(i);
+ }
+
+ //Put first hit in time of each cluster in a new map mapOfFirstN
+ //And also find the first hit of all : FirstN
+ static map< unsigned int, unsigned int > mapOfFirstN;
+ static unsigned int NbrOfClust, FirstN;
+ static int FirstClust;
+ NbrOfClust = 0;
+ NbrOfClust = mapOfClust.size();
+ static double GeneralTmin;
+ GeneralTmin = 1000000;
+ for(int i = 1; i <= NbrOfClust; i++){
+ static unsigned int clustSize;
+ static double TminClust;
+ TminClust = 10000000;
+ clustSize = mapOfClust[i].size();
+ for(int j = 0; j < clustSize; j++){
+ if((*HitT)[mapOfClust[i][j]] < TminClust){
+ TminClust = (*HitT)[mapOfClust[i][j]];
+ mapOfFirstN[i] = mapOfClust[i][j];
+ }
+ if((*HitT)[mapOfClust[i][j]] < GeneralTmin){
+ GeneralTmin = (*HitT)[mapOfClust[i][j]];
+ FirstN = mapOfClust[i][j];
+ FirstClust = i;
+ }
+ }
+ }
+
+ FirstHit = FirstN;
+ static double v_n, Dmax;
+ static vector<double> CrossTalk;
+ CrossTalk.clear();
+ //Test now the "Friend" clusters
+ for(int i = 1; i < NbrOfClust-1; i++){
+ x1 = (*HitX)[mapOfFirstN[i]], y1 = (*HitY)[mapOfFirstN[i]], z1 = (*HitZ)[mapOfFirstN[i]], t1 = (*HitT)[mapOfFirstN[i]];
+ v_n = sqrt(x1*x1+y1*y1+z1*z1)/t1;
+ for(int j = i+1; i < NbrOfClust; i++){
+ x2 = (*HitX)[mapOfFirstN[i]], y2 = (*HitY)[mapOfFirstN[i]], z2 = (*HitZ)[mapOfFirstN[i]], t2 = (*HitT)[mapOfFirstN[i]];
+ Dist = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
+ Dmax = (t2-t1)*v_n;
+ cout << Dmax << endl;
+ if(Dist < Dmax){
+ if(t1 < t2){
+ CrossTalk.push_back(j);
+ }
+ else{
+ CrossTalk.push_back(i);
+ }
+ }
+ }
+ }
+
+ //elimate potential CrossTalk in mapOfFirstN
+ static unsigned int sizeCT;
+ sizeCT = CrossTalk.size();
+ for(unsigned int i = 0; i < sizeCT; i++ ){
+ mapOfFirstN.erase(i);
+ }
+
+ //return vector of real neutrons IDs
+ m_Neutrons.clear();
+ for(auto itr = mapOfFirstN.begin(); itr != mapOfFirstN.end(); itr++){
+ m_Neutrons.push_back(itr->second);
+ }
+
+ return m_Neutrons;
+
+
+}
+
+int CrossTalk::GetFirstN(){
+ return FirstHit;
+}
+
diff --git a/NPLib/Neutron/NPCrossTalk.h b/NPLib/Neutron/NPCrossTalk.h
new file mode 100644
index 0000000000000000000000000000000000000000..30323548b93a0b91b1300de0df13ad4dbccaabfb
--- /dev/null
+++ b/NPLib/Neutron/NPCrossTalk.h
@@ -0,0 +1,62 @@
+#ifndef NPCROSSTALK_H
+#define NPCROSSTALK_H
+/*****************************************************************************
+ * Copyright (C) 2009-2020 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 : Cyril LENAIN contact address: lenain@lpcaen.in2p3.fr *
+ * *
+ * Creation Date : November 2020 *
+ * Last update : November 2020 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class find multi-neutrons in an neutron detection array by *
+ * rejecting "crosstalk" *
+ * Notation: *
+ *****************************************************************************/
+#include <vector>
+#include "TVector3.h"
+#include "Math/Functor.h"
+#include <map>
+
+class TGraph;
+namespace NPL{
+
+ class CrossTalk{
+ public:
+ CrossTalk();
+ ~CrossTalk();
+
+ public:
+
+ void AddHitVector(const std::vector<double>& X, const std::vector<double>& Y,const std::vector<double>& Z, const std::vector<double>& dX, const std::vector<double>& dY, const std::vector<double>& dZ, const std::vector<double>& T);
+
+ std::vector<int> ComputeCrossTalk();
+
+ int GetFirstN();
+
+ private: // private member used by
+ ROOT::Math::Functor m_func;
+ const std::vector<double>* HitX;
+ const std::vector<double>* HitY;
+ const std::vector<double>* HitZ;
+ const std::vector<double>* HitdX;
+ const std::vector<double>* HitdY;
+ const std::vector<double>* HitdZ;
+ const std::vector<double>* HitT;
+
+ std::vector<int> ClustHit;
+ std::vector<int> m_Neutrons;
+ double coef;
+ unsigned int sizeHit ;
+ int FirstHit;
+
+ };
+}
+
+#endif
diff --git a/NPLib/Online/NPOnlineGUI.cxx b/NPLib/Online/NPOnlineGUI.cxx
index 77ea3c58f2fd9c8f24128a7b173fcbbfaf73140a..9cb9304aebbfa61c3e8aea802110fe39f9f642ed 100644
--- a/NPLib/Online/NPOnlineGUI.cxx
+++ b/NPLib/Online/NPOnlineGUI.cxx
@@ -136,26 +136,28 @@ void NPL::OnlineGUI::ResetAll(){
}
////////////////////////////////////////////////////////////////////////////////
void NPL::OnlineGUI::ResetCurrent(){
- TCanvas* c = m_EmbeddedCanvas->GetCanvas();
+ TCanvas* c = m_EmbeddedCanvas->GetCanvas();
if (!c)
return;
+ // reset log scale attribute
gPad->SetLogx(false);
gPad->SetLogy(false);
gPad->SetLogz(false);
-
- TList* list = gPad->GetListOfPrimitives();
- int Hsize = list->GetSize();
- for(int h = 0 ; h < Hsize ; h++){
- TObject* obj = list->At(h);
- if(obj->InheritsFrom(TH1::Class())){
+ // loop on histograms, reset axis and content
+ TList* list = gPad->GetListOfPrimitives();
+ int Hsize = list->GetSize();
+ for(int h = 0 ; h < Hsize ; h++){
+ TObject* obj = list->At(h);
+ if(obj->InheritsFrom(TH1::Class())){
TH1* h = (TH1*) obj;
h->GetXaxis()->UnZoom();
h->GetYaxis()->UnZoom();
h->GetZaxis()->UnZoom();
- }
- }
+ h->Reset("ICESM");
+ }
+ }
c->Update();
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/NPSimulation/Detectors/MUST2/MUST2Array.cc b/NPSimulation/Detectors/MUST2/MUST2Array.cc
index 117b02d53a3179190a62479d9bdee6f78bd3643a..1545d2bdce5bf84f798016164862679157e26a07 100644
--- a/NPSimulation/Detectors/MUST2/MUST2Array.cc
+++ b/NPSimulation/Detectors/MUST2/MUST2Array.cc
@@ -695,8 +695,110 @@ void MUST2Array::ReadConfiguration(NPL::InputParser parser) {
if (blocks[i]->GetString("VIS") == "all")
m_non_sensitive_part_visiualisation = true;
}
+
+ ////////////////////
+ //Read the thresholds from the analysis config
+ ////////////////////
+
+ bool ReadingStatus = false;
+
+ // path to file
+ string FileName = "./configs/ConfigMust2.dat";
+
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigMust2.dat found: Default parameters loaded for "
+ "Analysis "
+ << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameters for Analysis from ConfigMust2.dat " << endl;
+
+ // read analysis config file
+ string LineBuffer, DataBuffer, whatToDo;
+
+ while (!AnalysisConfigFile.eof()) {
+ // Pick-up next line
+ getline(AnalysisConfigFile, LineBuffer);
+ // search for "header"
+ if (LineBuffer.compare(0, 11, "ConfigMust2") == 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');
+ }
+ //Resolutions
+ else if (whatToDo == "SI_E_RESOLUTION") {
+ AnalysisConfigFile >> DataBuffer;
+ ResoStrip = atof(DataBuffer.c_str());
+ ResoStrip = ResoStrip*keV/2.35;
+ cout << whatToDo << " " << ResoStrip << " MeV/2.35 "<< endl;
+ }
+ else if (whatToDo == "SILI_E_RESOLUTION") {
+ AnalysisConfigFile >> DataBuffer;
+ ResoSiLi = atof(DataBuffer.c_str());
+ ResoSiLi = ResoSiLi*keV/2.35;
+ cout << whatToDo << " " << ResoSiLi << " MeV/2.35 "<< endl;
+ }
+ else if (whatToDo == "CSI_E_RESOLUTION") {
+ AnalysisConfigFile >> DataBuffer;
+ ResoCsI = atof(DataBuffer.c_str());
+ ResoCsI = ResoCsI*keV/2.35;
+ cout << whatToDo << " " << ResoCsI << " MeV/2.35 "<< endl;
+ }
+ //Time
+ else if (whatToDo == "MUST_T_RESOLUTION") {
+ AnalysisConfigFile >> DataBuffer;
+ ResoTimeMust = atof(DataBuffer.c_str());
+ ResoTimeMust = ResoTimeMust*ns/2.35;
+ cout << whatToDo << " " << ResoTimeMust << " ns/2.35 "<< endl;
+ }
+ else if (whatToDo == "SI_T_OFFSET") {
+ AnalysisConfigFile >> DataBuffer;
+ TimeOffset = atof(DataBuffer.c_str());
+ TimeOffset = TimeOffset*ns;
+ cout << whatToDo << " " << TimeOffset << " ns "<< endl;
+ }
+ //Thresholds
+ else if (whatToDo == "SI_X_E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ ThresholdSiX = atof(DataBuffer.c_str());
+ ThresholdSiX = ThresholdSiX*keV;
+ cout << whatToDo << " " << ThresholdSiX << " MeV "<< endl;
+ }
+ else if (whatToDo == "SI_Y_E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ ThresholdSiY = atof(DataBuffer.c_str());
+ ThresholdSiY = ThresholdSiY*keV;
+ cout << whatToDo << " " << ThresholdSiY << " MeV "<< endl;
+ }
+ else if (whatToDo == "SILI_E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ ThresholdSiLi = atof(DataBuffer.c_str());
+ ThresholdSiLi = ThresholdSiLi*keV;
+ cout << whatToDo << " " << ThresholdSiLi << " MeV "<< endl;
+ }
+ else if (whatToDo == "CSI_E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ ThresholdCsI = atof(DataBuffer.c_str());
+ ThresholdCsI = ThresholdCsI*keV;
+ cout << whatToDo << " " << ThresholdCsI << " MeV "<< endl;
+ }
+ else if (AnalysisConfigFile.eof()) ReadingStatus = false;
+ }
+ }
}
+
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
void MUST2Array::ConstructDetector(G4LogicalVolume* world) {
@@ -860,7 +962,7 @@ void MUST2Array::ReadSensitive(const G4Event*) {
double timeX = TimeOffset - RandGauss::shoot(it->second.second, ResoTimeMust);
unsigned int strip = it->first-1000000*(it->first/1000000);
unsigned int det = it->first/1000000;
- if (energyX > ThresholdSi) {
+ if (energyX > ThresholdSiX) {
trig.insert(det);
SiScoredHit = true;
m_Event->SetStripXE(det, strip ,
@@ -906,7 +1008,7 @@ void MUST2Array::ReadSensitive(const G4Event*) {
double timeY = TimeOffset - RandGauss::shoot(it->second.second, ResoTimeMust);
unsigned int strip = it->first-1000000*(it->first/1000000);
unsigned int det = it->first/1000000;
- if (energyY > ThresholdSi) {
+ if (energyY > ThresholdSiY) {
trig.insert(det);
SiScoredHit = true;
m_Event->SetStripYE(det, strip ,
diff --git a/NPSimulation/Detectors/MUST2/MUST2Array.hh b/NPSimulation/Detectors/MUST2/MUST2Array.hh
index 36f1efc84ee7e170f5db7a38cd219fdd25aecfda..ad1c00a80c41706d888faff7f5a134fbaa814bcc 100644
--- a/NPSimulation/Detectors/MUST2/MUST2Array.hh
+++ b/NPSimulation/Detectors/MUST2/MUST2Array.hh
@@ -34,16 +34,17 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace MUST2
{
- // Resolution
- const G4double ResoTimeMust = 0.212765957 ;// = 500ps // Unit is ns/2.35
- const G4double ResoSiLi = 0.028 ;// = 130 keV of resolution // Unit is MeV/2.35
- const G4double ResoCsI = 0.08 ;// = 188 kev of resolution // Unit is MeV/2.35
- const G4double ResoStrip = 0.0149 ;// 0.0223 = 52keV of Resolution // Unit is MeV/2.35 14.861996
- const G4double TimeOffset = 500 ;// 500 ns stop
- // Threshold
- const G4double ThresholdSi = 500 * keV;
- const G4double ThresholdSiLi = 500 * keV;
- const G4double ThresholdCsI = 500 * keV;
+ // Default Resolution
+ G4double ResoTimeMust = 0.2127 ;// = 500ps // Unit is ns/2.35
+ G4double ResoStrip = 0.022 ;// 0.0223 = 52keV of Resolution // Unit is MeV/2.35 14.861996
+ G4double ResoSiLi = 0.055 ;// = 130 keV of resolution // Unit is MeV/2.35
+ G4double ResoCsI = 0.080 ;// = 188 kev of resolution // Unit is MeV/2.35
+ G4double TimeOffset = 500 ;// 500 ns stop
+ // Default Threshold
+ G4double ThresholdSiX = 500 * keV;
+ G4double ThresholdSiY = 500 * keV;
+ G4double ThresholdSiLi = 200 * keV; // typically shielded by the DSSD
+ G4double ThresholdCsI = 200 * keV;
// Geometry
const G4double FaceFront = 11.*cm ;
const G4double FaceBack = 16.5*cm ;
diff --git a/Projects/ComptonTelescope/online/CanvasList.txt b/Projects/ComptonTelescope/online/CanvasList.txt
index aa772776412cce4c6f50343d190486cede5f7d7f..e8da3d622a0c2c31299a091bb357c01c3fd81f09 100644
--- a/Projects/ComptonTelescope/online/CanvasList.txt
+++ b/Projects/ComptonTelescope/online/CanvasList.txt
@@ -14,7 +14,7 @@ Canvas
Canvas
Path= Calorimeter
Divide= 1 2
- Histo= CT1_CALOR_SPECTRUM CT1_CALOR_RAW_TRIGGER CT1_CALOR_POS
+ Histo= CT1_CALOR_SPECTRUM CT1_CALOR_POS
Canvas
Path= ComptonTelescope
diff --git a/Projects/eAGanil/resolution.cxx b/Projects/eAGanil/resolution.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..eeb0d67b703a8ca46a3fd9c7d3ec33fd5e8b2140
--- /dev/null
+++ b/Projects/eAGanil/resolution.cxx
@@ -0,0 +1,63 @@
+void AddNuclei(string,double);
+void resolution(){
+
+ auto f = new TFile("../../Outputs/Analysis/PhysicsTree.root");
+ auto tree= (TTree*) f->FindObjectAny("PhysicsTree");
+
+ vector<double> Resolution={5e-2,1e-2,5e-3,1e-3,5e-4,1e-4,5e-5,1e-5,5e-6,1e-6};
+ unsigned int sizeR = Resolution.size();
+ auto c = new TCanvas();
+ c->Divide(3,2);
+ vector<double> R ;
+ for(unsigned int i = 0 ; i < sizeR ; i++){
+ c->cd(i+1);
+ if(i<3)
+ tree->Draw(Form("Ex>>h%d(1000,-10,10)",i),Form("Resolution==%f",Resolution[i]),"");
+ else
+ tree->Draw(Form("Ex>>h%d(1000,-0.5,0.5)",i),Form("Resolution==%f",Resolution[i]),"");
+
+ auto h = (TH1*) gDirectory->FindObjectAny(Form("h%d",i));
+ h->Fit("gaus");
+ R.push_back(2*sqrt(2*log(2))*h->GetFunction("gaus")->GetParameter(2)*1000.);
+ cout << R[i] << endl;
+ }
+
+ auto g = new TGraph(R.size(),&Resolution[0],&R[0]);
+ new TCanvas();
+ g->Draw("apc");
+ gPad->SetLogx();
+
+ g->GetYaxis()->SetRangeUser(-100,1400);
+ g->GetXaxis()->SetTitle("E/#deltaE");
+ g->GetYaxis()->SetTitle("FWHM (keV)");
+
+
+ AddNuclei("134Sn",725.6);
+ AddNuclei("132Sn",64.4);
+ AddNuclei("133Sn",853.7);
+ AddNuclei("70Ni",183.11);
+ AddNuclei("71Ni",252.2);
+ AddNuclei("72Ni",455);
+ //AddNuclei("73Ni",239.2);
+ AddNuclei("74Ni",739);
+}
+
+void AddNuclei(string name,double keV){
+ static bool alt = true;
+ static double off =1e-6;
+ auto L = new TLine(1e-6,keV,5e-2,keV);
+ L->SetLineStyle(2);
+ L->SetLineColor(kAzure+7);
+ L->Draw("same");
+ TLatex* t;
+ if(alt){
+ t = new TLatex(3e-6-off,keV-30,name.c_str());
+ alt=false;
+ }
+ else{
+ t = new TLatex(3e-6+off,keV+10,name.c_str());
+ alt=true;
+ }
+ t->SetTextSize(0.02);
+ t->Draw();
+}