mirror of
https://github.com/gwm17/Mask.git
synced 2024-11-22 18:28:51 -05:00
Compare commits
5 Commits
c12e649ea3
...
25fee0f39b
| Author | SHA1 | Date | |
|---|---|---|---|
|
25fee0f39b | ||
 Gordon McCann
|
8b81dda70e | ||
|
|
3db9e510f5 | ||
|
Gordon McCann
|
2f7b0969dc | ||
|
Gordon McCann
|
0456f230f5 |
|
|
@ -34,11 +34,13 @@ Detector geometry is encoded using ROOT math libraries in the `src/Detectors` fo
|
||||||
|
|
||||||
To choose which detector scheme is run, modify the main function in `src/Detectors/main.cpp`. The included geometries also have options to do an internal geometry consistency check and print out coordinates for drawing the detector arrays, which can be useful for testing.
|
To choose which detector scheme is run, modify the main function in `src/Detectors/main.cpp`. The included geometries also have options to do an internal geometry consistency check and print out coordinates for drawing the detector arrays, which can be useful for testing.
|
||||||
|
|
||||||
|
To run the geometry code, one needs to provide an input file containing the following: the path of a Mask kinematics data file, the path to which data should be written, the path to a file containing a list of dead channels (optional, if not used, write None for the path), the number of threads to be used by the thread pool, and a keyword for the array type (current options are Sabre or Anasen)
|
||||||
|
|
||||||
To run Detectors use the format
|
To run Detectors use the format
|
||||||
|
|
||||||
`./bin/Detectors <kinematics_datafile> <new_detection_datafile> <new_detection_statsfile>`
|
`./bin/Detectors <input_file>`
|
||||||
|
|
||||||
where the detection datafile contains all of the kinematics data as well as information about which particles are detected and the statsfile is a text file containing efficiency statistics.
|
An example input file is provided with the repository.
|
||||||
|
|
||||||
## Data visualization
|
## Data visualization
|
||||||
All data is saved as ROOT trees of std::vectors of Mask::Nucleus classes. To enable this, a ROOT dictionary is generated and linked into a shared library found in the `lib` directory of the repository. This allows the user to link to the shared library for accessing and analyzing the data generated by Mask.
|
All data is saved as ROOT trees of std::vectors of Mask::Nucleus classes. To enable this, a ROOT dictionary is generated and linked into a shared library found in the `lib` directory of the repository. This allows the user to link to the shared library for accessing and analyzing the data generated by Mask.
|
||||||
|
|
|
||||||
5
detector_input.txt
Normal file
5
detector_input.txt
Normal file
|
|
@ -0,0 +1,5 @@
|
||||||
|
InputDataFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
|
||||||
|
OutputDataFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B_Sabre.root
|
||||||
|
DeadChannelFile: etc/sabreDeadChannels_May2022.txt
|
||||||
|
NumberOfThreads: 5
|
||||||
|
ArrayType: Sabre
|
||||||
|
|
@ -1,7 +1,8 @@
|
||||||
----------Data Information----------
|
----------Data Information----------
|
||||||
OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
|
OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
|
||||||
|
NumberOfThreads: 6
|
||||||
----------Reaction Information----------
|
----------Reaction Information----------
|
||||||
NumberOfSamples: 100000
|
NumberOfSamples: 1000000
|
||||||
begin_chain
|
begin_chain
|
||||||
begin_step
|
begin_step
|
||||||
Type: Reaction
|
Type: Reaction
|
||||||
|
|
@ -230,7 +230,7 @@ double AnasenArray::RunConsistencyCheck()
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
DetectorResult AnasenArray::IsRing1(Mask::Nucleus& nucleus)
|
DetectorResult AnasenArray::IsRing1(const Mask::Nucleus& nucleus)
|
||||||
{
|
{
|
||||||
DetectorResult observation;
|
DetectorResult observation;
|
||||||
double thetaIncident;
|
double thetaIncident;
|
||||||
|
|
@ -255,7 +255,7 @@ DetectorResult AnasenArray::IsRing1(Mask::Nucleus& nucleus)
|
||||||
return observation;
|
return observation;
|
||||||
}
|
}
|
||||||
|
|
||||||
DetectorResult AnasenArray::IsRing2(Mask::Nucleus& nucleus)
|
DetectorResult AnasenArray::IsRing2(const Mask::Nucleus& nucleus)
|
||||||
{
|
{
|
||||||
DetectorResult observation;
|
DetectorResult observation;
|
||||||
double thetaIncident;
|
double thetaIncident;
|
||||||
|
|
@ -280,7 +280,7 @@ DetectorResult AnasenArray::IsRing2(Mask::Nucleus& nucleus)
|
||||||
return observation;
|
return observation;
|
||||||
}
|
}
|
||||||
|
|
||||||
DetectorResult AnasenArray::IsQQQ(Mask::Nucleus& nucleus)
|
DetectorResult AnasenArray::IsQQQ(const Mask::Nucleus& nucleus)
|
||||||
{
|
{
|
||||||
DetectorResult observation;
|
DetectorResult observation;
|
||||||
double thetaIncident;
|
double thetaIncident;
|
||||||
|
|
@ -324,7 +324,7 @@ DetectorResult AnasenArray::IsQQQ(Mask::Nucleus& nucleus)
|
||||||
return observation;
|
return observation;
|
||||||
}
|
}
|
||||||
|
|
||||||
DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus)
|
DetectorResult AnasenArray::IsDetected(const Mask::Nucleus& nucleus)
|
||||||
{
|
{
|
||||||
DetectorResult result;
|
DetectorResult result;
|
||||||
if(nucleus.GetKE() <= s_energyThreshold)
|
if(nucleus.GetKE() <= s_energyThreshold)
|
||||||
|
|
@ -337,236 +337,4 @@ DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus)
|
||||||
if(!result.detectFlag)
|
if(!result.detectFlag)
|
||||||
result = IsQQQ(nucleus);
|
result = IsQQQ(nucleus);
|
||||||
return result;
|
return result;
|
||||||
}
|
|
||||||
|
|
||||||
void AnasenArray::CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts)
|
|
||||||
{
|
|
||||||
if (data.size() == 3 && data[1].isDetected && data[2].isDetected)
|
|
||||||
{
|
|
||||||
counts[0]++;
|
|
||||||
}
|
|
||||||
else if (data.size() == 4 && data[2].isDetected && data[3].isDetected)
|
|
||||||
{
|
|
||||||
counts[0]++;
|
|
||||||
}
|
|
||||||
else if(data.size() == 6)
|
|
||||||
{
|
|
||||||
if(data[2].isDetected && data[4].isDetected)
|
|
||||||
{
|
|
||||||
counts[0]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[5].isDetected)
|
|
||||||
{
|
|
||||||
counts[1]++;
|
|
||||||
}
|
|
||||||
if(data[4].isDetected && data[5].isDetected)
|
|
||||||
{
|
|
||||||
counts[2]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[4].isDetected && data[5].isDetected)
|
|
||||||
{
|
|
||||||
counts[3]++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
else if(data.size() == 8)
|
|
||||||
{
|
|
||||||
if(data[2].isDetected && data[4].isDetected)
|
|
||||||
{
|
|
||||||
counts[0]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[6].isDetected)
|
|
||||||
{
|
|
||||||
counts[1]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[2]++;
|
|
||||||
}
|
|
||||||
if(data[4].isDetected && data[6].isDetected)
|
|
||||||
{
|
|
||||||
counts[3]++;
|
|
||||||
}
|
|
||||||
if(data[4].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[4]++;
|
|
||||||
}
|
|
||||||
if(data[6].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[5]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[4].isDetected && data[6].isDetected)
|
|
||||||
{
|
|
||||||
counts[6]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[4].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[7]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[6].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[8]++;
|
|
||||||
}
|
|
||||||
if(data[4].isDetected && data[6].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[9]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[4].isDetected && data[6].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[10]++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void AnasenArray::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) {
|
|
||||||
|
|
||||||
if(!dmap.IsValid())
|
|
||||||
{
|
|
||||||
std::cerr<<"Invalid Dead Channel Map at AnasenArray::CalculateEfficiency()! If you want to run with all possible";
|
|
||||||
std::cerr<<"channels active, simply load an empty file."<<std::endl;
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
std::cout<<"----------ANASEN Efficiency Calculation----------"<<std::endl;
|
|
||||||
std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
|
|
||||||
std::cout<<"Running efficiency calculation..."<<std::endl;
|
|
||||||
|
|
||||||
TFile* input = TFile::Open(inputname.c_str(), "READ");
|
|
||||||
TFile* output = TFile::Open(outputname.c_str(), "RECREATE");
|
|
||||||
std::ofstream stats(statsname);
|
|
||||||
stats<<std::setprecision(5);
|
|
||||||
|
|
||||||
TTree* intree = (TTree*) input->Get("SimTree");
|
|
||||||
std::vector<Mask::Nucleus>* dataHandle = new std::vector<Mask::Nucleus>();
|
|
||||||
intree->SetBranchAddress("nuclei", &dataHandle);
|
|
||||||
|
|
||||||
output->cd();
|
|
||||||
TTree* outtree = new TTree("SimTree", "SimTree");
|
|
||||||
outtree->Branch("nuclei", dataHandle);
|
|
||||||
|
|
||||||
input->cd();
|
|
||||||
|
|
||||||
stats<<"Efficiency statistics for data from "<<inputname<<" using the ANASEN geometry"<<std::endl;
|
|
||||||
stats<<"Given in order of target=0, projectile=1, ejectile=2, residual=3, .... etc."<<std::endl;
|
|
||||||
|
|
||||||
intree->GetEntry(1);
|
|
||||||
std::vector<int> counts;
|
|
||||||
std::vector<int> coinc_counts;
|
|
||||||
counts.resize(dataHandle->size());
|
|
||||||
switch(counts.size())
|
|
||||||
{
|
|
||||||
case 3: coinc_counts.resize(1, 0); break;
|
|
||||||
case 4: coinc_counts.resize(1, 0); break;
|
|
||||||
case 6: coinc_counts.resize(4, 0); break;
|
|
||||||
case 8: coinc_counts.resize(11, 0); break;
|
|
||||||
default:
|
|
||||||
{
|
|
||||||
std::cerr<<"Bad reaction type at AnasenArray::CalculateEfficiency (given value: "<<counts.size()<<"). Quiting..."<<std::endl;
|
|
||||||
input->Close();
|
|
||||||
output->Close();
|
|
||||||
stats.close();
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
uint64_t nentries = intree->GetEntries();
|
|
||||||
uint64_t percent5 = nentries*0.05;
|
|
||||||
uint64_t count = 0;
|
|
||||||
uint64_t npercent = 0;
|
|
||||||
|
|
||||||
Mask::Nucleus nucleus;
|
|
||||||
for(uint64_t i=0; i<nentries; i++)
|
|
||||||
{
|
|
||||||
intree->GetEntry(i);
|
|
||||||
if(++count == percent5)
|
|
||||||
{//Show progress every 5%
|
|
||||||
npercent++;
|
|
||||||
count = 0;
|
|
||||||
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
|
|
||||||
}
|
|
||||||
|
|
||||||
for(std::size_t j=0; j<dataHandle->size(); j++)
|
|
||||||
{
|
|
||||||
Mask::Nucleus& nucleus = (*dataHandle)[j];
|
|
||||||
DetectorResult result = IsAnasen(nucleus);
|
|
||||||
if(result.detectFlag)
|
|
||||||
{
|
|
||||||
nucleus.isDetected = true;
|
|
||||||
nucleus.detectedKE = result.energy_deposited;
|
|
||||||
nucleus.detectedTheta = result.direction.Theta();
|
|
||||||
nucleus.detectedPhi = result.direction.Phi();
|
|
||||||
counts[j]++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
CountCoincidences(*dataHandle, coinc_counts);
|
|
||||||
|
|
||||||
outtree->Fill();
|
|
||||||
}
|
|
||||||
input->Close();
|
|
||||||
output->cd();
|
|
||||||
outtree->Write(outtree->GetName(), TObject::kOverwrite);
|
|
||||||
output->Close();
|
|
||||||
|
|
||||||
delete dataHandle;
|
|
||||||
|
|
||||||
stats<<std::setw(10)<<"Index"<<"|"<<std::setw(10)<<"Efficiency"<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
for(unsigned int i=0; i<counts.size(); i++)
|
|
||||||
{
|
|
||||||
stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
stats<<"Coincidence Efficiency"<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
if(dataHandle->size() == 3)
|
|
||||||
{
|
|
||||||
stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
else if(dataHandle->size() == 4)
|
|
||||||
{
|
|
||||||
stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
else if(dataHandle->size() == 6)
|
|
||||||
{
|
|
||||||
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
else if(dataHandle->size() == 8)
|
|
||||||
{
|
|
||||||
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
stats.close();
|
|
||||||
|
|
||||||
std::cout<<std::endl;
|
|
||||||
std::cout<<"Complete."<<std::endl;
|
|
||||||
std::cout<<"---------------------------------------------"<<std::endl;
|
|
||||||
}
|
}
|
||||||
|
|
@ -6,8 +6,8 @@
|
||||||
#include "DetectorArray.h"
|
#include "DetectorArray.h"
|
||||||
#include "SX3Detector.h"
|
#include "SX3Detector.h"
|
||||||
#include "QQQDetector.h"
|
#include "QQQDetector.h"
|
||||||
#include "Target.h"
|
#include "Mask/Target.h"
|
||||||
#include "Nucleus.h"
|
#include "Mask/Nucleus.h"
|
||||||
#include "AnasenDeadChannelMap.h"
|
#include "AnasenDeadChannelMap.h"
|
||||||
|
|
||||||
class AnasenArray : public DetectorArray
|
class AnasenArray : public DetectorArray
|
||||||
|
|
@ -15,17 +15,15 @@ class AnasenArray : public DetectorArray
|
||||||
public:
|
public:
|
||||||
AnasenArray();
|
AnasenArray();
|
||||||
~AnasenArray();
|
~AnasenArray();
|
||||||
virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
|
virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus);
|
||||||
virtual void DrawDetectorSystem(const std::string& filename) override;
|
virtual void DrawDetectorSystem(const std::string& filename) override;
|
||||||
virtual double RunConsistencyCheck() override;
|
virtual double RunConsistencyCheck() override;
|
||||||
inline void SetDeadChannelMap(const std::string& filename) { dmap.LoadMapfile(filename); }
|
virtual void SetDeadChannelMap(const std::string& filename) override { dmap.LoadMapfile(filename); }
|
||||||
|
|
||||||
private:
|
private:
|
||||||
DetectorResult IsRing1(Mask::Nucleus& nucleus);
|
DetectorResult IsRing1(const Mask::Nucleus& nucleus);
|
||||||
DetectorResult IsRing2(Mask::Nucleus& nucleus);
|
DetectorResult IsRing2(const Mask::Nucleus& nucleus);
|
||||||
DetectorResult IsQQQ(Mask::Nucleus& nucleus);
|
DetectorResult IsQQQ(const Mask::Nucleus& nucleus);
|
||||||
DetectorResult IsAnasen(Mask::Nucleus& nucleus);
|
|
||||||
void CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts);
|
|
||||||
|
|
||||||
std::vector<SX3Detector> m_Ring1;
|
std::vector<SX3Detector> m_Ring1;
|
||||||
std::vector<SX3Detector> m_Ring2;
|
std::vector<SX3Detector> m_Ring2;
|
||||||
|
|
|
||||||
|
|
@ -18,6 +18,9 @@ target_sources(Detectors PUBLIC
|
||||||
SabreArray.h
|
SabreArray.h
|
||||||
SX3Detector.cpp
|
SX3Detector.cpp
|
||||||
SX3Detector.h
|
SX3Detector.h
|
||||||
|
DetectorApp.h
|
||||||
|
DetectorApp.cpp
|
||||||
|
DetectorArray.cpp
|
||||||
)
|
)
|
||||||
|
|
||||||
target_link_libraries(Detectors
|
target_link_libraries(Detectors
|
||||||
|
|
|
||||||
144
src/Detectors/DetectorApp.cpp
Normal file
144
src/Detectors/DetectorApp.cpp
Normal file
|
|
@ -0,0 +1,144 @@
|
||||||
|
#include "DetectorApp.h"
|
||||||
|
#include <fstream>
|
||||||
|
#include <iostream>
|
||||||
|
|
||||||
|
DetectorApp::DetectorApp() :
|
||||||
|
m_resources(nullptr)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
DetectorApp::~DetectorApp()
|
||||||
|
{
|
||||||
|
for(std::size_t i=0; i<m_detectorList.size(); i++)
|
||||||
|
delete m_detectorList[i];
|
||||||
|
}
|
||||||
|
|
||||||
|
bool DetectorApp::LoadConfig(const std::string& filename)
|
||||||
|
{
|
||||||
|
std::cout<<"----------Detector Efficiency Calculation----------"<<std::endl;
|
||||||
|
std::ifstream input(filename);
|
||||||
|
if(!input.is_open())
|
||||||
|
{
|
||||||
|
std::cerr<<"Unable to open input config file "<<filename<<std::endl;
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::string junk;
|
||||||
|
std::string typeName;
|
||||||
|
|
||||||
|
input >> junk >> m_inputFileName;
|
||||||
|
input >> junk >> m_outputFileName;
|
||||||
|
input >> junk >> m_deadChannelFileName;
|
||||||
|
|
||||||
|
input >> junk >> m_nthreads;
|
||||||
|
input >> junk >> typeName;
|
||||||
|
|
||||||
|
std::cout << "Creating " << m_nthreads << " detector arrays..." << std::endl;
|
||||||
|
ArrayType type = StringToArrayType(typeName);
|
||||||
|
for(uint64_t i=0; i<m_nthreads; i++)
|
||||||
|
{
|
||||||
|
m_detectorList.push_back(CreateDetectorArray(type));
|
||||||
|
if(m_deadChannelFileName != "None")
|
||||||
|
m_detectorList.back()->SetDeadChannelMap(m_deadChannelFileName);
|
||||||
|
}
|
||||||
|
std::cout << "Done" << std::endl;
|
||||||
|
|
||||||
|
std::cout << "Allocating " << m_nthreads << " threads..." << std::endl;
|
||||||
|
m_resources = std::make_unique<Mask::ThreadPool<DetectorArray*, uint64_t>>(m_nthreads);
|
||||||
|
std::cout << "Done" << std::endl;
|
||||||
|
|
||||||
|
std::cout << "Opening input data file " << m_inputFileName << "..." << std::endl;
|
||||||
|
m_fileReader.Open(m_inputFileName, "SimTree");
|
||||||
|
if(!m_fileReader.IsOpen() || !m_fileReader.IsTree())
|
||||||
|
{
|
||||||
|
std::cerr << "Unable to open input data file " << m_inputFileName << std::endl;
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
m_nentries = m_fileReader.GetSize();
|
||||||
|
std::cout << "Done. Detected " << m_nentries << " events in the file." << std::endl;
|
||||||
|
|
||||||
|
//Little bit of integer division mangling to make sure we read every event in file
|
||||||
|
uint64_t quotient = m_nentries / m_nthreads;
|
||||||
|
uint64_t remainder = m_nentries % m_nthreads;
|
||||||
|
m_chunkSamples.push_back(quotient + remainder);
|
||||||
|
for(uint64_t i=1; i<m_nthreads; i++)
|
||||||
|
m_chunkSamples.push_back(quotient);
|
||||||
|
|
||||||
|
std::cout << "Opening output data file " << m_outputFileName << std::endl;
|
||||||
|
m_fileWriter.Open(m_outputFileName, "SimTree");
|
||||||
|
if(!m_fileWriter.IsOpen() || !m_fileWriter.IsTree())
|
||||||
|
{
|
||||||
|
std::cerr << "Unable to open output data file " << m_outputFileName << std::endl;
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
std::cout << "Done" << std::endl;
|
||||||
|
|
||||||
|
std::cout << "Ready to launch!" << std::endl;
|
||||||
|
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
|
void DetectorApp::Run()
|
||||||
|
{
|
||||||
|
std::cout<<"Running efficiency calculation..."<<std::endl;
|
||||||
|
|
||||||
|
if(m_detectorList.size() != m_nthreads)
|
||||||
|
{
|
||||||
|
std::cerr << "Detector list not equal to number of threads" << std::endl;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
for(uint64_t i=0; i<m_detectorList.size(); i++)
|
||||||
|
{
|
||||||
|
//Create a job for the thread pool, using a lambda and providing a tuple of the arguments
|
||||||
|
m_resources->PushJob({[this](DetectorArray* array, uint64_t chunkSamples)
|
||||||
|
{
|
||||||
|
if(system == nullptr)
|
||||||
|
return;
|
||||||
|
|
||||||
|
std::vector<Mask::Nucleus> data;
|
||||||
|
DetectorResult result;
|
||||||
|
for(uint64_t i=0; i<chunkSamples; i++)
|
||||||
|
{
|
||||||
|
m_fileReader.Read(data);
|
||||||
|
for(auto& nucleus : data)
|
||||||
|
{
|
||||||
|
result = array->IsDetected(nucleus);
|
||||||
|
if(result.detectFlag)
|
||||||
|
{
|
||||||
|
nucleus.isDetected = true;
|
||||||
|
nucleus.detectedKE = result.energy_deposited;
|
||||||
|
nucleus.detectedTheta = result.direction.Theta();
|
||||||
|
nucleus.detectedPhi = result.direction.Phi();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
m_fileWriter.PushData(data);
|
||||||
|
}
|
||||||
|
},
|
||||||
|
{m_detectorList[i], m_chunkSamples[i]} //arguments to function, in order
|
||||||
|
}
|
||||||
|
);
|
||||||
|
}
|
||||||
|
|
||||||
|
uint64_t size = m_fileReader.GetSize();
|
||||||
|
uint64_t count = 0;
|
||||||
|
double percent = 0.05;
|
||||||
|
uint64_t flushVal = size*percent;
|
||||||
|
uint64_t flushCount = 0;
|
||||||
|
|
||||||
|
while(true)
|
||||||
|
{
|
||||||
|
if(count == flushVal)
|
||||||
|
{
|
||||||
|
count = 0;
|
||||||
|
++flushCount;
|
||||||
|
std::cout<<"\rPercent of data written to disk: "<<percent*flushCount*100<<"%"<<std::flush;
|
||||||
|
}
|
||||||
|
|
||||||
|
if(m_resources->IsFinished() && m_fileWriter.GetQueueSize() == 0)
|
||||||
|
break;
|
||||||
|
else if(m_fileWriter.Write())
|
||||||
|
++count;
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
40
src/Detectors/DetectorApp.h
Normal file
40
src/Detectors/DetectorApp.h
Normal file
|
|
@ -0,0 +1,40 @@
|
||||||
|
#ifndef DETECTOR_APP_H
|
||||||
|
#define DETECTOR_APP_H
|
||||||
|
|
||||||
|
#include "DetectorArray.h"
|
||||||
|
#include "Mask/FileWriter.h"
|
||||||
|
#include "Mask/FileReader.h"
|
||||||
|
#include "Mask/ThreadPool.h"
|
||||||
|
|
||||||
|
#include <string>
|
||||||
|
#include <vector>
|
||||||
|
#include <memory>
|
||||||
|
|
||||||
|
class DetectorApp
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
DetectorApp();
|
||||||
|
~DetectorApp();
|
||||||
|
|
||||||
|
bool LoadConfig(const std::string& filename);
|
||||||
|
|
||||||
|
void Run();
|
||||||
|
|
||||||
|
private:
|
||||||
|
|
||||||
|
std::vector<DetectorArray*> m_detectorList; //One array per thread
|
||||||
|
std::vector<uint64_t> m_chunkSamples;
|
||||||
|
Mask::FileWriter m_fileWriter;
|
||||||
|
Mask::FileReader m_fileReader;
|
||||||
|
|
||||||
|
std::string m_inputFileName;
|
||||||
|
std::string m_outputFileName;
|
||||||
|
std::string m_deadChannelFileName;
|
||||||
|
|
||||||
|
uint64_t m_nthreads;
|
||||||
|
uint64_t m_nentries;
|
||||||
|
|
||||||
|
std::unique_ptr<Mask::ThreadPool<DetectorArray*, uint64_t>> m_resources;
|
||||||
|
};
|
||||||
|
|
||||||
|
#endif
|
||||||
35
src/Detectors/DetectorArray.cpp
Normal file
35
src/Detectors/DetectorArray.cpp
Normal file
|
|
@ -0,0 +1,35 @@
|
||||||
|
#include "DetectorArray.h"
|
||||||
|
#include "AnasenArray.h"
|
||||||
|
#include "SabreArray.h"
|
||||||
|
|
||||||
|
DetectorArray* CreateDetectorArray(ArrayType type)
|
||||||
|
{
|
||||||
|
switch(type)
|
||||||
|
{
|
||||||
|
case ArrayType::None: return nullptr;
|
||||||
|
case ArrayType::Anasen: return new AnasenArray();
|
||||||
|
case ArrayType::Sabre: return new SabreArray();
|
||||||
|
}
|
||||||
|
return nullptr;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::string ArrayTypeToString(ArrayType type)
|
||||||
|
{
|
||||||
|
switch(type)
|
||||||
|
{
|
||||||
|
case ArrayType::None: return "None";
|
||||||
|
case ArrayType::Anasen: return "Anasen";
|
||||||
|
case ArrayType::Sabre: return "Sabre";
|
||||||
|
}
|
||||||
|
return "None";
|
||||||
|
}
|
||||||
|
|
||||||
|
ArrayType StringToArrayType(const std::string& value)
|
||||||
|
{
|
||||||
|
if (value == "Anasen")
|
||||||
|
return ArrayType::Anasen;
|
||||||
|
else if (value == "Sabre")
|
||||||
|
return ArrayType::Sabre;
|
||||||
|
else
|
||||||
|
return ArrayType::None;
|
||||||
|
}
|
||||||
|
|
@ -5,6 +5,7 @@
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
|
|
||||||
#include "Math/Point3D.h"
|
#include "Math/Point3D.h"
|
||||||
|
#include "Mask/Nucleus.h"
|
||||||
|
|
||||||
struct DetectorResult
|
struct DetectorResult
|
||||||
{
|
{
|
||||||
|
|
@ -14,15 +15,23 @@ struct DetectorResult
|
||||||
std::string det_name = "";
|
std::string det_name = "";
|
||||||
};
|
};
|
||||||
|
|
||||||
|
enum class ArrayType
|
||||||
|
{
|
||||||
|
None,
|
||||||
|
Anasen,
|
||||||
|
Sabre
|
||||||
|
};
|
||||||
|
|
||||||
class DetectorArray
|
class DetectorArray
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
DetectorArray() {};
|
DetectorArray() {};
|
||||||
virtual ~DetectorArray() {};
|
virtual ~DetectorArray() {};
|
||||||
|
|
||||||
virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) = 0;
|
virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus) = 0;
|
||||||
virtual void DrawDetectorSystem(const std::string& filename) = 0;
|
virtual void DrawDetectorSystem(const std::string& filename) = 0;
|
||||||
virtual double RunConsistencyCheck() = 0;
|
virtual double RunConsistencyCheck() = 0;
|
||||||
|
virtual void SetDeadChannelMap(const std::string& filename) = 0;
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < s_epsilon ? true : false; };
|
bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < s_epsilon ? true : false; };
|
||||||
|
|
@ -30,4 +39,10 @@ protected:
|
||||||
static constexpr double s_epsilon = 1.0e-6;
|
static constexpr double s_epsilon = 1.0e-6;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
DetectorArray* CreateDetectorArray(ArrayType type);
|
||||||
|
|
||||||
|
std::string ArrayTypeToString(ArrayType type);
|
||||||
|
|
||||||
|
ArrayType StringToArrayType(const std::string& value);
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
|
@ -31,161 +31,6 @@ SabreArray::SabreArray() :
|
||||||
|
|
||||||
SabreArray::~SabreArray() {}
|
SabreArray::~SabreArray() {}
|
||||||
|
|
||||||
void SabreArray::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname)
|
|
||||||
{
|
|
||||||
std::cout<<"----------SABRE Efficiency Calculation----------"<<std::endl;
|
|
||||||
std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
|
|
||||||
std::cout<<"Running efficiency calculation..."<<std::endl;
|
|
||||||
|
|
||||||
if(!m_deadMap.IsValid())
|
|
||||||
{
|
|
||||||
std::cerr<<"Unable to run SABRE Efficiency without a dead channel map."<<std::endl;
|
|
||||||
std::cerr<<"If you have no dead channels, simply make a file that's empty"<<std::endl;
|
|
||||||
std::cerr<<"Exiting."<<std::endl;
|
|
||||||
std::cout<<"---------------------------------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
TFile* input = TFile::Open(inputname.c_str(), "READ");
|
|
||||||
TFile* output = TFile::Open(outputname.c_str(), "RECREATE");
|
|
||||||
std::ofstream stats(statsname);
|
|
||||||
stats<<std::setprecision(5);
|
|
||||||
|
|
||||||
TTree* intree = (TTree*) input->Get("SimTree");
|
|
||||||
std::vector<Mask::Nucleus>* dataHandle = new std::vector<Mask::Nucleus>();
|
|
||||||
intree->SetBranchAddress("nuclei", &dataHandle);
|
|
||||||
|
|
||||||
output->cd();
|
|
||||||
TTree* outtree = new TTree("SimTree", "SimTree");
|
|
||||||
outtree->Branch("nuclei", dataHandle);
|
|
||||||
|
|
||||||
input->cd();
|
|
||||||
|
|
||||||
stats<<"Efficiency statistics for data from "<<inputname<<" using the ANASEN geometry"<<std::endl;
|
|
||||||
stats<<"Given in order of target=0, projectile=1, ejectile=2, residual=3, .... etc."<<std::endl;
|
|
||||||
|
|
||||||
intree->GetEntry(1);
|
|
||||||
std::vector<int> counts;
|
|
||||||
std::vector<int> coinc_counts;
|
|
||||||
counts.resize(dataHandle->size());
|
|
||||||
switch(counts.size())
|
|
||||||
{
|
|
||||||
case 3: coinc_counts.resize(1, 0); break;
|
|
||||||
case 4: coinc_counts.resize(1, 0); break;
|
|
||||||
case 6: coinc_counts.resize(4, 0); break;
|
|
||||||
case 8: coinc_counts.resize(11, 0); break;
|
|
||||||
default:
|
|
||||||
{
|
|
||||||
std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<<counts.size()<<"). Quiting..."<<std::endl;
|
|
||||||
input->Close();
|
|
||||||
output->Close();
|
|
||||||
stats.close();
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
uint64_t nentries = intree->GetEntries();
|
|
||||||
uint64_t percent5 = nentries*0.05;
|
|
||||||
uint64_t count = 0;
|
|
||||||
uint64_t npercent = 0;
|
|
||||||
|
|
||||||
Mask::Nucleus nucleus;
|
|
||||||
for(uint64_t i=0; i<nentries; i++)
|
|
||||||
{
|
|
||||||
intree->GetEntry(i);
|
|
||||||
if(++count == percent5)
|
|
||||||
{//Show progress every 5%
|
|
||||||
npercent++;
|
|
||||||
count = 0;
|
|
||||||
std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
|
|
||||||
}
|
|
||||||
|
|
||||||
for(std::size_t j=0; j<dataHandle->size(); j++)
|
|
||||||
{
|
|
||||||
Mask::Nucleus& nucleus = (*dataHandle)[j];
|
|
||||||
DetectorResult result = IsSabre(nucleus);
|
|
||||||
if(result.detectFlag)
|
|
||||||
{
|
|
||||||
nucleus.isDetected = true;
|
|
||||||
nucleus.detectedKE = result.energy_deposited;
|
|
||||||
nucleus.detectedTheta = result.direction.Theta();
|
|
||||||
nucleus.detectedPhi = result.direction.Phi();
|
|
||||||
counts[j]++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
CountCoincidences(*dataHandle, coinc_counts);
|
|
||||||
|
|
||||||
outtree->Fill();
|
|
||||||
}
|
|
||||||
input->Close();
|
|
||||||
output->cd();
|
|
||||||
outtree->Write(outtree->GetName(), TObject::kOverwrite);
|
|
||||||
output->Close();
|
|
||||||
|
|
||||||
delete dataHandle;
|
|
||||||
|
|
||||||
stats<<std::setw(10)<<"Index"<<"|"<<std::setw(10)<<"Efficiency"<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
for(unsigned int i=0; i<counts.size(); i++) {
|
|
||||||
stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
stats<<"Coincidence Efficiency"<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
if(counts.size() == 3)
|
|
||||||
{
|
|
||||||
stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
else if(counts.size() == 4)
|
|
||||||
{
|
|
||||||
stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
else if(counts.size() == 6)
|
|
||||||
{
|
|
||||||
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
else if(counts.size() == 8)
|
|
||||||
{
|
|
||||||
stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/nentries)<<std::endl;
|
|
||||||
stats<<"---------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
stats.close();
|
|
||||||
|
|
||||||
std::cout<<std::endl;
|
|
||||||
std::cout<<"Complete."<<std::endl;
|
|
||||||
std::cout<<"---------------------------------------------"<<std::endl;
|
|
||||||
}
|
|
||||||
|
|
||||||
void SabreArray::DrawDetectorSystem(const std::string& filename)
|
void SabreArray::DrawDetectorSystem(const std::string& filename)
|
||||||
{
|
{
|
||||||
std::ofstream output(filename);
|
std::ofstream output(filename);
|
||||||
|
|
@ -259,7 +104,7 @@ double SabreArray::RunConsistencyCheck()
|
||||||
}
|
}
|
||||||
|
|
||||||
/*Returns if detected, as well as total energy deposited in SABRE*/
|
/*Returns if detected, as well as total energy deposited in SABRE*/
|
||||||
DetectorResult SabreArray::IsSabre(Mask::Nucleus& nucleus)
|
DetectorResult SabreArray::IsDetected(const Mask::Nucleus& nucleus)
|
||||||
{
|
{
|
||||||
DetectorResult observation;
|
DetectorResult observation;
|
||||||
if(nucleus.GetKE() <= s_energyThreshold)
|
if(nucleus.GetKE() <= s_energyThreshold)
|
||||||
|
|
@ -297,82 +142,4 @@ DetectorResult SabreArray::IsSabre(Mask::Nucleus& nucleus)
|
||||||
|
|
||||||
observation.detectFlag = false;
|
observation.detectFlag = false;
|
||||||
return observation;
|
return observation;
|
||||||
}
|
}
|
||||||
|
|
||||||
void SabreArray::CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts)
|
|
||||||
{
|
|
||||||
if (data.size() == 3 && data[1].isDetected && data[2].isDetected)
|
|
||||||
{
|
|
||||||
counts[0]++;
|
|
||||||
}
|
|
||||||
else if (data.size() == 4 && data[2].isDetected && data[3].isDetected)
|
|
||||||
{
|
|
||||||
counts[0]++;
|
|
||||||
}
|
|
||||||
else if(data.size() == 6)
|
|
||||||
{
|
|
||||||
if(data[2].isDetected && data[4].isDetected)
|
|
||||||
{
|
|
||||||
counts[0]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[5].isDetected)
|
|
||||||
{
|
|
||||||
counts[1]++;
|
|
||||||
}
|
|
||||||
if(data[4].isDetected && data[5].isDetected)
|
|
||||||
{
|
|
||||||
counts[2]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[4].isDetected && data[5].isDetected)
|
|
||||||
{
|
|
||||||
counts[3]++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
else if(data.size() == 8)
|
|
||||||
{
|
|
||||||
if(data[2].isDetected && data[4].isDetected)
|
|
||||||
{
|
|
||||||
counts[0]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[6].isDetected)
|
|
||||||
{
|
|
||||||
counts[1]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[2]++;
|
|
||||||
}
|
|
||||||
if(data[4].isDetected && data[6].isDetected)
|
|
||||||
{
|
|
||||||
counts[3]++;
|
|
||||||
}
|
|
||||||
if(data[4].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[4]++;
|
|
||||||
}
|
|
||||||
if(data[6].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[5]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[4].isDetected && data[6].isDetected)
|
|
||||||
{
|
|
||||||
counts[6]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[4].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[7]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[6].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[8]++;
|
|
||||||
}
|
|
||||||
if(data[4].isDetected && data[6].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[9]++;
|
|
||||||
}
|
|
||||||
if(data[2].isDetected && data[4].isDetected && data[6].isDetected && data[7].isDetected)
|
|
||||||
{
|
|
||||||
counts[10]++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
@ -3,7 +3,7 @@
|
||||||
|
|
||||||
#include "DetectorArray.h"
|
#include "DetectorArray.h"
|
||||||
#include "SabreDetector.h"
|
#include "SabreDetector.h"
|
||||||
#include "Target.h"
|
#include "Mask/Target.h"
|
||||||
#include "SabreDeadChannelMap.h"
|
#include "SabreDeadChannelMap.h"
|
||||||
#include "Mask/Nucleus.h"
|
#include "Mask/Nucleus.h"
|
||||||
|
|
||||||
|
|
@ -12,14 +12,12 @@ class SabreArray : public DetectorArray
|
||||||
public:
|
public:
|
||||||
SabreArray();
|
SabreArray();
|
||||||
~SabreArray();
|
~SabreArray();
|
||||||
void SetDeadChannelMap(const std::string& filename) { m_deadMap.LoadMapfile(filename); };
|
virtual void SetDeadChannelMap(const std::string& filename) override { m_deadMap.LoadMapfile(filename); };
|
||||||
void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
|
virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus) override;
|
||||||
void DrawDetectorSystem(const std::string& filename) override;
|
void DrawDetectorSystem(const std::string& filename) override;
|
||||||
double RunConsistencyCheck() override;
|
double RunConsistencyCheck() override;
|
||||||
|
|
||||||
private:
|
private:
|
||||||
DetectorResult IsSabre(Mask::Nucleus& nucleus);
|
|
||||||
void CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts);
|
|
||||||
|
|
||||||
std::vector<SabreDetector> m_detectors;
|
std::vector<SabreDetector> m_detectors;
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -1,5 +1,4 @@
|
||||||
#include "SabreArray.h"
|
#include "DetectorApp.h"
|
||||||
#include "AnasenArray.h"
|
|
||||||
#include "KinematicsExceptions.h"
|
#include "KinematicsExceptions.h"
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
#include <string>
|
#include <string>
|
||||||
|
|
@ -7,7 +6,7 @@
|
||||||
int main(int argc, char** argv)
|
int main(int argc, char** argv)
|
||||||
{
|
{
|
||||||
|
|
||||||
if(argc != 4)
|
if(argc != 2)
|
||||||
{
|
{
|
||||||
std::cerr<<"Incorrect number of commandline arguments! Returning."<<std::endl;
|
std::cerr<<"Incorrect number of commandline arguments! Returning."<<std::endl;
|
||||||
return 1;
|
return 1;
|
||||||
|
|
@ -19,26 +18,21 @@ int main(int argc, char** argv)
|
||||||
return 1;
|
return 1;
|
||||||
}
|
}
|
||||||
|
|
||||||
std::string inputname = argv[1];
|
try
|
||||||
std::string outputname = argv[2];
|
{
|
||||||
std::string statsname = argv[3];
|
DetectorApp app;
|
||||||
|
if(!app.LoadConfig(argv[1]))
|
||||||
|
{
|
||||||
|
std::cerr << "Unable to load config file " << argv[1] << ". Shutting down." << std::endl;
|
||||||
|
return 1;
|
||||||
|
}
|
||||||
|
|
||||||
|
app.Run();
|
||||||
SabreArray sabre;
|
}
|
||||||
std::string mapfile = "./etc/sabreDeadChannels_May2022.txt";
|
catch(std::exception& e)
|
||||||
sabre.SetDeadChannelMap(mapfile);
|
{
|
||||||
sabre.CalculateEfficiency(inputname, outputname, statsname);
|
std::cerr << e.what() << std::endl;
|
||||||
//std::cout<<"Running consistency check(0=success): "<<sabre.RunConsistencyCheck()<<std::endl;
|
return 1;
|
||||||
//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.txt");
|
}
|
||||||
|
|
||||||
|
|
||||||
/*
|
|
||||||
AnasenArray anasen;
|
|
||||||
std::string mapfile = "./etc/AnasenDeadChannels.txt";
|
|
||||||
anasen.SetDeadChannelMap(mapfile);
|
|
||||||
anasen.CalculateEfficiency(inputname, outputname, statsname);
|
|
||||||
//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
|
|
||||||
//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_fixedZ.txt");
|
|
||||||
*/
|
|
||||||
return 0;
|
return 0;
|
||||||
}
|
}
|
||||||
|
|
|
||||||
|
|
@ -24,6 +24,7 @@ int main(int argc, char** argv)
|
||||||
return 1;
|
return 1;
|
||||||
}
|
}
|
||||||
calculator.Run();
|
calculator.Run();
|
||||||
|
//calculator.RunSingleThread();
|
||||||
}
|
}
|
||||||
catch(const std::exception& e)
|
catch(const std::exception& e)
|
||||||
{
|
{
|
||||||
|
|
|
||||||
|
|
@ -50,7 +50,15 @@ target_sources(Mask PRIVATE
|
||||||
ThreeStepSystem.h
|
ThreeStepSystem.h
|
||||||
TwoStepSystem.cpp
|
TwoStepSystem.cpp
|
||||||
TwoStepSystem.h
|
TwoStepSystem.h
|
||||||
|
ThreadPool.h
|
||||||
|
FileWriter.h
|
||||||
|
FileWriter.cpp
|
||||||
|
FileReader.h
|
||||||
|
FileReader.cpp
|
||||||
)
|
)
|
||||||
|
|
||||||
target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES})
|
set(THREADS_PREFER_PTHREAD_FLAG On)
|
||||||
|
find_package(Threads REQUIRED)
|
||||||
|
|
||||||
|
target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES} Threads::Threads)
|
||||||
set_target_properties(Mask PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})
|
set_target_properties(Mask PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})
|
||||||
66
src/Mask/FileReader.cpp
Normal file
66
src/Mask/FileReader.cpp
Normal file
|
|
@ -0,0 +1,66 @@
|
||||||
|
#include "FileReader.h"
|
||||||
|
|
||||||
|
namespace Mask {
|
||||||
|
|
||||||
|
FileReader::FileReader() :
|
||||||
|
m_file(nullptr), m_tree(nullptr), m_branchHandle(nullptr), m_currentEntry(0), m_size(0)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
FileReader::FileReader(const std::string& filename, const std::string& treename) :
|
||||||
|
m_file(nullptr), m_tree(nullptr), m_branchHandle(nullptr), m_currentEntry(0), m_size(0)
|
||||||
|
{
|
||||||
|
Open(filename, treename);
|
||||||
|
}
|
||||||
|
|
||||||
|
FileReader::~FileReader()
|
||||||
|
{
|
||||||
|
Close();
|
||||||
|
}
|
||||||
|
|
||||||
|
void FileReader::Open(const std::string& filename, const std::string& treename)
|
||||||
|
{
|
||||||
|
if(m_file != nullptr || m_tree != nullptr)
|
||||||
|
Close();
|
||||||
|
|
||||||
|
m_file = TFile::Open(filename.c_str(), "READ");
|
||||||
|
if(m_file != nullptr && m_file->IsOpen())
|
||||||
|
{
|
||||||
|
m_tree = (TTree*) m_file->Get(treename.c_str());
|
||||||
|
if(m_tree == nullptr)
|
||||||
|
{
|
||||||
|
m_file->Close();
|
||||||
|
delete m_file;
|
||||||
|
m_file = nullptr;
|
||||||
|
}
|
||||||
|
m_branchHandle = new std::vector<Nucleus>();
|
||||||
|
m_tree->SetBranchAddress("nuclei", &m_branchHandle);
|
||||||
|
m_size = m_tree->GetEntries();
|
||||||
|
m_currentEntry = 0; //Reset file position
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void FileReader::Close()
|
||||||
|
{
|
||||||
|
if(m_file != nullptr && m_file->IsOpen())
|
||||||
|
{
|
||||||
|
m_file->Close();
|
||||||
|
delete m_file;
|
||||||
|
m_file = nullptr;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
bool FileReader::Read(std::vector<Nucleus>& dataHandle)
|
||||||
|
{
|
||||||
|
std::scoped_lock<std::mutex> guard(m_fileMutex);
|
||||||
|
int bytes = m_tree->GetEntry(m_currentEntry);
|
||||||
|
if(bytes != 0)
|
||||||
|
{
|
||||||
|
dataHandle = *m_branchHandle;
|
||||||
|
m_currentEntry++;
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
}
|
||||||
46
src/Mask/FileReader.h
Normal file
46
src/Mask/FileReader.h
Normal file
|
|
@ -0,0 +1,46 @@
|
||||||
|
#ifndef FILE_READER_H
|
||||||
|
#define FILE_READER_H
|
||||||
|
|
||||||
|
#include "Nucleus.h"
|
||||||
|
|
||||||
|
#include "TFile.h"
|
||||||
|
#include "TTree.h"
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
#include <string>
|
||||||
|
#include <mutex>
|
||||||
|
#include <atomic>
|
||||||
|
|
||||||
|
namespace Mask {
|
||||||
|
|
||||||
|
class FileReader
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
FileReader();
|
||||||
|
FileReader(const std::string& filename, const std::string& treename);
|
||||||
|
~FileReader();
|
||||||
|
|
||||||
|
void Open(const std::string& filename, const std::string& treename); //Not thread safe
|
||||||
|
void Close(); //Not thread safe
|
||||||
|
|
||||||
|
/*
|
||||||
|
Read: fills entry to given dataHandle. Returns true if data was successfully filled, otherwise returns false
|
||||||
|
*/
|
||||||
|
bool Read(std::vector<Nucleus>& dataHandle); //Thread safe
|
||||||
|
uint64_t GetSize() { return m_size; }//In entries (implicitly thread safe)
|
||||||
|
bool IsOpen() { return m_file == nullptr ? false : m_file->IsOpen(); } //Should be safe?
|
||||||
|
bool IsTree() { return m_tree != nullptr; } //Should be safe?
|
||||||
|
|
||||||
|
private:
|
||||||
|
TFile* m_file;
|
||||||
|
TTree* m_tree;
|
||||||
|
|
||||||
|
std::vector<Nucleus>* m_branchHandle;
|
||||||
|
|
||||||
|
std::mutex m_fileMutex;
|
||||||
|
std::atomic<uint64_t> m_currentEntry;
|
||||||
|
std::atomic<uint64_t> m_size; //in entries
|
||||||
|
};
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
75
src/Mask/FileWriter.cpp
Normal file
75
src/Mask/FileWriter.cpp
Normal file
|
|
@ -0,0 +1,75 @@
|
||||||
|
#include "FileWriter.h"
|
||||||
|
|
||||||
|
namespace Mask {
|
||||||
|
|
||||||
|
FileWriter::FileWriter() :
|
||||||
|
m_file(nullptr), m_tree(nullptr), m_queueSize(0)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
FileWriter::FileWriter(const std::string& filename, const std::string& treename) :
|
||||||
|
m_file(nullptr), m_tree(nullptr)
|
||||||
|
{
|
||||||
|
m_file = TFile::Open(filename.c_str(), "RECREATE");
|
||||||
|
if(m_file != nullptr && m_file->IsOpen())
|
||||||
|
{
|
||||||
|
m_tree = new TTree(treename.c_str(), treename.c_str());
|
||||||
|
m_tree->Branch("nuclei", &m_dataHandle);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
FileWriter::~FileWriter()
|
||||||
|
{
|
||||||
|
Close();
|
||||||
|
}
|
||||||
|
|
||||||
|
void FileWriter::Open(const std::string& filename, const std::string& treename)
|
||||||
|
{
|
||||||
|
if(m_file != nullptr || m_tree != nullptr)
|
||||||
|
Close();
|
||||||
|
|
||||||
|
m_file = TFile::Open(filename.c_str(), "RECREATE");
|
||||||
|
if(m_file != nullptr && m_file->IsOpen())
|
||||||
|
{
|
||||||
|
m_tree = new TTree(treename.c_str(), treename.c_str());
|
||||||
|
m_tree->Branch("nuclei", &m_dataHandle);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void FileWriter::Close()
|
||||||
|
{
|
||||||
|
if(m_file->IsOpen())
|
||||||
|
{
|
||||||
|
if(m_tree != nullptr)
|
||||||
|
m_tree->Write(m_tree->GetName(), TObject::kOverwrite);
|
||||||
|
|
||||||
|
m_file->Close();
|
||||||
|
delete m_file;
|
||||||
|
m_file = nullptr;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void FileWriter::PushData(const std::vector<Nucleus>& data)
|
||||||
|
{
|
||||||
|
std::scoped_lock<std::mutex> guard(m_queueMutex);
|
||||||
|
m_queue.push(data);
|
||||||
|
++m_queueSize;
|
||||||
|
}
|
||||||
|
|
||||||
|
bool FileWriter::Write()
|
||||||
|
{
|
||||||
|
if(m_queueSize == 0)
|
||||||
|
return false;
|
||||||
|
|
||||||
|
//Aquire lock for as short a time as possible
|
||||||
|
{
|
||||||
|
std::scoped_lock<std::mutex> guard(m_queueMutex);
|
||||||
|
m_dataHandle = m_queue.front();
|
||||||
|
m_queue.pop();
|
||||||
|
}
|
||||||
|
--m_queueSize;
|
||||||
|
|
||||||
|
m_tree->Fill();
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
}
|
||||||
45
src/Mask/FileWriter.h
Normal file
45
src/Mask/FileWriter.h
Normal file
|
|
@ -0,0 +1,45 @@
|
||||||
|
#ifndef FILE_WRITER_H
|
||||||
|
#define FILE_WRITER_H
|
||||||
|
|
||||||
|
#include "Nucleus.h"
|
||||||
|
|
||||||
|
#include "TFile.h"
|
||||||
|
#include "TTree.h"
|
||||||
|
|
||||||
|
#include <queue>
|
||||||
|
#include <mutex>
|
||||||
|
#include <atomic>
|
||||||
|
|
||||||
|
namespace Mask {
|
||||||
|
|
||||||
|
class FileWriter
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
FileWriter();
|
||||||
|
FileWriter(const std::string& filename, const std::string& treename);
|
||||||
|
~FileWriter();
|
||||||
|
|
||||||
|
bool IsOpen() const { return m_file == nullptr ? false : m_file->IsOpen(); }
|
||||||
|
bool IsTree() const { return m_tree == nullptr ? false : true; }
|
||||||
|
|
||||||
|
std::size_t GetQueueSize() const { return m_queueSize; } //Implicitly thread-safe
|
||||||
|
|
||||||
|
void PushData(const std::vector<Nucleus>& data); //Thread-safe
|
||||||
|
bool Write(); //Not completely thread-safe, should only be used by main application loop (acess to queue is safe, but all else not)
|
||||||
|
|
||||||
|
void Open(const std::string& filename, const std::string& treename); //Not thread safe!
|
||||||
|
void Close(); //Not thread safe!
|
||||||
|
|
||||||
|
private:
|
||||||
|
TFile* m_file;
|
||||||
|
TTree* m_tree;
|
||||||
|
|
||||||
|
std::vector<Nucleus> m_dataHandle;
|
||||||
|
|
||||||
|
std::mutex m_queueMutex;
|
||||||
|
std::atomic<std::size_t> m_queueSize;
|
||||||
|
std::queue<std::vector<Nucleus>> m_queue;
|
||||||
|
};
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
|
|
@ -8,7 +8,7 @@
|
||||||
namespace Mask {
|
namespace Mask {
|
||||||
|
|
||||||
MaskApp::MaskApp() :
|
MaskApp::MaskApp() :
|
||||||
m_system(nullptr)
|
m_system(nullptr), m_resources(nullptr)
|
||||||
{
|
{
|
||||||
std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
|
std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
|
||||||
}
|
}
|
||||||
|
|
@ -16,6 +16,8 @@ namespace Mask {
|
||||||
MaskApp::~MaskApp()
|
MaskApp::~MaskApp()
|
||||||
{
|
{
|
||||||
delete m_system;
|
delete m_system;
|
||||||
|
for(std::size_t i=0; i<m_systemList.size(); i++)
|
||||||
|
delete m_systemList[i];
|
||||||
}
|
}
|
||||||
|
|
||||||
bool MaskApp::LoadConfig(const std::string& filename)
|
bool MaskApp::LoadConfig(const std::string& filename)
|
||||||
|
|
@ -31,10 +33,25 @@ namespace Mask {
|
||||||
std::string junk;
|
std::string junk;
|
||||||
std::getline(input, junk);
|
std::getline(input, junk);
|
||||||
input>>junk>>m_outputName;
|
input>>junk>>m_outputName;
|
||||||
|
input>>junk>>m_nthreads;
|
||||||
std::getline(input, junk);
|
std::getline(input, junk);
|
||||||
std::getline(input, junk);
|
std::getline(input, junk);
|
||||||
input>>junk>>m_nsamples;
|
input>>junk>>m_nsamples;
|
||||||
|
|
||||||
|
std::cout<<"Allocating resources... Asking for " << m_nthreads << " threads...";
|
||||||
|
m_resources = std::make_unique<ThreadPool<ReactionSystem*, uint64_t>>(m_nthreads);
|
||||||
|
std::cout<<" Complete."<<std::endl;
|
||||||
|
|
||||||
|
//Little bit of integer division mangling to make sure we do the total number of samples
|
||||||
|
uint64_t quotient = m_nsamples / m_nthreads;
|
||||||
|
uint64_t remainder = m_nsamples % m_nthreads;
|
||||||
|
m_chunkSamples.push_back(quotient + remainder);
|
||||||
|
for(uint64_t i=1; i<m_nthreads; i++)
|
||||||
|
m_chunkSamples.push_back(quotient);
|
||||||
|
|
||||||
|
std::cout<<"Outputing data to file: " << m_outputName <<std::endl;
|
||||||
|
m_fileWriter.Open(m_outputName, "SimTree");
|
||||||
|
|
||||||
std::vector<StepParameters> params;
|
std::vector<StepParameters> params;
|
||||||
int z, a;
|
int z, a;
|
||||||
while(input>>junk)
|
while(input>>junk)
|
||||||
|
|
@ -98,10 +115,20 @@ namespace Mask {
|
||||||
m_system = CreateSystem(params);
|
m_system = CreateSystem(params);
|
||||||
if(m_system == nullptr || !m_system->IsValid())
|
if(m_system == nullptr || !m_system->IsValid())
|
||||||
{
|
{
|
||||||
std::cerr<<"Param size: "<<params.size()<<std::endl;
|
|
||||||
std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
|
std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
for(uint32_t i=0; i<m_nthreads; i++)
|
||||||
|
{
|
||||||
|
m_systemList.push_back(CreateSystem(params));
|
||||||
|
if(m_systemList.back() == nullptr || !m_systemList.back()->IsValid())
|
||||||
|
{
|
||||||
|
std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
std::getline(input, junk);
|
std::getline(input, junk);
|
||||||
std::getline(input, junk);
|
std::getline(input, junk);
|
||||||
|
|
||||||
|
|
@ -132,9 +159,13 @@ namespace Mask {
|
||||||
}
|
}
|
||||||
m_system->SetLayeredTarget(target);
|
m_system->SetLayeredTarget(target);
|
||||||
|
|
||||||
std::cout<<"Outputing data to file: "<<m_outputName<<std::endl;
|
for(auto system : m_systemList)
|
||||||
std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
|
{
|
||||||
std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
|
system->SetLayeredTarget(target);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::cout<<"Reaction equation: "<<m_system->GetSystemEquation()<<std::endl;
|
||||||
|
std::cout<<"Number of samples: "<<m_nsamples<<std::endl;
|
||||||
|
|
||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
|
|
@ -144,8 +175,59 @@ namespace Mask {
|
||||||
{
|
{
|
||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
|
|
||||||
void MaskApp::Run()
|
void MaskApp::Run()
|
||||||
|
{
|
||||||
|
std::cout<<"Running simulation..."<<std::endl;
|
||||||
|
if(m_systemList.size() != m_nthreads)
|
||||||
|
{
|
||||||
|
std::cerr << "System list not equal to number of threads" << std::endl;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
//Give our thread pool some tasks
|
||||||
|
for(std::size_t i=0; i<m_systemList.size(); i++)
|
||||||
|
{
|
||||||
|
//bind a lambda to the job, taking in a ReactionSystem, and then provide a reaction system as the tuple arguments.
|
||||||
|
m_resources->PushJob({[this](ReactionSystem* system, uint64_t chunkSamples)
|
||||||
|
{
|
||||||
|
if(system == nullptr)
|
||||||
|
return;
|
||||||
|
|
||||||
|
for(uint64_t i=0; i<chunkSamples; i++)
|
||||||
|
{
|
||||||
|
system->RunSystem();
|
||||||
|
m_fileWriter.PushData(*(system->GetNuclei()));
|
||||||
|
}
|
||||||
|
},
|
||||||
|
{m_systemList[i], m_chunkSamples[i]}});
|
||||||
|
}
|
||||||
|
|
||||||
|
uint64_t count = 0;
|
||||||
|
double percent = 0.05;
|
||||||
|
uint64_t flushVal = m_nsamples*percent;
|
||||||
|
uint64_t flushCount = 0;
|
||||||
|
while(true)
|
||||||
|
{
|
||||||
|
if(count == flushVal)
|
||||||
|
{
|
||||||
|
count = 0;
|
||||||
|
++flushCount;
|
||||||
|
std::cout<<"\rPercent of data written to disk: "<<percent*flushCount*100<<"%"<<std::flush;
|
||||||
|
}
|
||||||
|
|
||||||
|
if(m_resources->IsFinished() && m_fileWriter.GetQueueSize() == 0)
|
||||||
|
break;
|
||||||
|
else if(m_fileWriter.Write())
|
||||||
|
++count;
|
||||||
|
}
|
||||||
|
|
||||||
|
std::cout<<std::endl;
|
||||||
|
std::cout<<"Complete."<<std::endl;
|
||||||
|
std::cout<<"---------------------------------------------"<<std::endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
void MaskApp::RunSingleThread()
|
||||||
{
|
{
|
||||||
std::cout<<"Running simulation..."<<std::endl;
|
std::cout<<"Running simulation..."<<std::endl;
|
||||||
if(m_system == nullptr)
|
if(m_system == nullptr)
|
||||||
|
|
@ -158,11 +240,11 @@ namespace Mask {
|
||||||
tree->Branch("nuclei", m_system->GetNuclei());
|
tree->Branch("nuclei", m_system->GetNuclei());
|
||||||
|
|
||||||
//For progress tracking
|
//For progress tracking
|
||||||
uint32_t percent5 = 0.05*m_nsamples;
|
uint64_t percent5 = 0.05*m_nsamples;
|
||||||
uint32_t count = 0;
|
uint64_t count = 0;
|
||||||
uint32_t npercent = 0;
|
uint64_t npercent = 0;
|
||||||
|
|
||||||
for(uint32_t i=0; i<m_nsamples; i++)
|
for(uint64_t i=0; i<m_nsamples; i++)
|
||||||
{
|
{
|
||||||
if(++count == percent5)
|
if(++count == percent5)
|
||||||
{
|
{
|
||||||
|
|
|
||||||
|
|
@ -7,6 +7,10 @@
|
||||||
#include "TwoStepSystem.h"
|
#include "TwoStepSystem.h"
|
||||||
#include "ThreeStepSystem.h"
|
#include "ThreeStepSystem.h"
|
||||||
#include "RxnType.h"
|
#include "RxnType.h"
|
||||||
|
#include "ThreadPool.h"
|
||||||
|
#include "FileWriter.h"
|
||||||
|
|
||||||
|
#include <memory>
|
||||||
|
|
||||||
namespace Mask {
|
namespace Mask {
|
||||||
|
|
||||||
|
|
@ -17,18 +21,24 @@ namespace Mask {
|
||||||
~MaskApp();
|
~MaskApp();
|
||||||
bool LoadConfig(const std::string& filename);
|
bool LoadConfig(const std::string& filename);
|
||||||
bool SaveConfig(const std::string& filename);
|
bool SaveConfig(const std::string& filename);
|
||||||
int GetNumberOfSamples() const { return m_nsamples; }
|
|
||||||
const std::string GetSystemName() const { return m_system == nullptr ? "Invalid System" : m_system->GetSystemEquation(); }
|
|
||||||
const std::string GetOutputName() const { return m_outputName; }
|
|
||||||
const RxnType GetReactionType() const { return m_rxnType; }
|
|
||||||
|
|
||||||
void Run();
|
void Run();
|
||||||
|
void RunSingleThread();
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
void RunChunk(ReactionSystem* system);
|
||||||
|
|
||||||
ReactionSystem* m_system;
|
ReactionSystem* m_system;
|
||||||
|
|
||||||
std::string m_outputName;
|
std::string m_outputName;
|
||||||
RxnType m_rxnType;
|
RxnType m_rxnType;
|
||||||
int m_nsamples;
|
uint64_t m_nsamples;
|
||||||
|
uint32_t m_nthreads;
|
||||||
|
|
||||||
|
std::vector<ReactionSystem*> m_systemList; //One system for each thread
|
||||||
|
std::vector<uint64_t> m_chunkSamples;
|
||||||
|
FileWriter m_fileWriter;
|
||||||
|
std::unique_ptr<ThreadPool<ReactionSystem*, uint64_t>> m_resources;
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -1,8 +1,7 @@
|
||||||
#include "RandomGenerator.h"
|
#include "RandomGenerator.h"
|
||||||
|
|
||||||
namespace Mask {
|
namespace Mask {
|
||||||
RandomGenerator* RandomGenerator::s_instance = new RandomGenerator();
|
|
||||||
|
|
||||||
RandomGenerator::RandomGenerator()
|
RandomGenerator::RandomGenerator()
|
||||||
{
|
{
|
||||||
std::random_device rd;
|
std::random_device rd;
|
||||||
|
|
|
||||||
|
|
@ -10,12 +10,15 @@ namespace Mask {
|
||||||
public:
|
public:
|
||||||
~RandomGenerator();
|
~RandomGenerator();
|
||||||
std::mt19937& GetGenerator() { return rng; }
|
std::mt19937& GetGenerator() { return rng; }
|
||||||
static RandomGenerator& GetInstance() { return *s_instance; }
|
static RandomGenerator& GetInstance()
|
||||||
|
{
|
||||||
|
static thread_local RandomGenerator s_instance;
|
||||||
|
return s_instance;
|
||||||
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
RandomGenerator();
|
RandomGenerator();
|
||||||
|
|
||||||
static RandomGenerator* s_instance;
|
|
||||||
std::mt19937 rng;
|
std::mt19937 rng;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
||||||
118
src/Mask/ThreadPool.h
Normal file
118
src/Mask/ThreadPool.h
Normal file
|
|
@ -0,0 +1,118 @@
|
||||||
|
/*
|
||||||
|
ThreadPool.h
|
||||||
|
Simple thread pool implementation. Jobs are q'd, and as threads are made available, jobs are dispatched to these threads.
|
||||||
|
|
||||||
|
GWM -- Aug 2022
|
||||||
|
*/
|
||||||
|
#ifndef THREAD_POOL_H
|
||||||
|
#define THREAD_POOL_H
|
||||||
|
|
||||||
|
#include "ReactionSystem.h"
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
#include <atomic>
|
||||||
|
#include <thread>
|
||||||
|
#include <mutex>
|
||||||
|
#include <functional>
|
||||||
|
#include <queue>
|
||||||
|
#include <condition_variable>
|
||||||
|
|
||||||
|
|
||||||
|
namespace Mask {
|
||||||
|
|
||||||
|
template<typename... Types>
|
||||||
|
class ThreadPool
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
|
||||||
|
using JobFunction = std::function<void(Types...)>;
|
||||||
|
struct Job
|
||||||
|
{
|
||||||
|
JobFunction func;
|
||||||
|
std::tuple<Types...> argument;
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
|
ThreadPool(uint32_t nthreads) :
|
||||||
|
m_isStopped(false), m_numberRunning(0), m_queueSize(0), m_initShutdown(false)
|
||||||
|
{
|
||||||
|
for(uint32_t i=0; i<nthreads; i++)
|
||||||
|
{
|
||||||
|
m_pool.push_back(std::thread(std::bind(&ThreadPool::ExecuteJob, std::ref(*this))));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
~ThreadPool()
|
||||||
|
{
|
||||||
|
if(!m_isStopped)
|
||||||
|
Shutdown();
|
||||||
|
}
|
||||||
|
|
||||||
|
void PushJob(Job job)
|
||||||
|
{
|
||||||
|
{
|
||||||
|
std::unique_lock<std::mutex> guard(m_poolMutex);
|
||||||
|
m_queue.push(job);
|
||||||
|
m_queueSize++;
|
||||||
|
}
|
||||||
|
m_wakeCondition.notify_one();
|
||||||
|
}
|
||||||
|
|
||||||
|
bool IsFinished()
|
||||||
|
{
|
||||||
|
return m_numberRunning == 0 && m_queueSize == 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
void Shutdown()
|
||||||
|
{
|
||||||
|
m_initShutdown = true;
|
||||||
|
m_wakeCondition.notify_all();
|
||||||
|
for(auto& thread : m_pool)
|
||||||
|
{
|
||||||
|
thread.join();
|
||||||
|
}
|
||||||
|
|
||||||
|
m_pool.clear();
|
||||||
|
m_isStopped = true;
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
void ExecuteJob()
|
||||||
|
{
|
||||||
|
while(true)
|
||||||
|
{
|
||||||
|
Job job;
|
||||||
|
{
|
||||||
|
std::unique_lock<std::mutex> guard(m_poolMutex);
|
||||||
|
m_wakeCondition.wait(guard, [this](){
|
||||||
|
return (!m_queue.empty() || m_initShutdown);
|
||||||
|
});
|
||||||
|
if(m_initShutdown)
|
||||||
|
return;
|
||||||
|
|
||||||
|
job = m_queue.front();
|
||||||
|
m_queue.pop();
|
||||||
|
}
|
||||||
|
|
||||||
|
//Change number running first so that no crash
|
||||||
|
m_numberRunning++;
|
||||||
|
m_queueSize--;
|
||||||
|
|
||||||
|
std::apply(job.func, job.argument);
|
||||||
|
m_numberRunning--;
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
std::queue<Job, std::deque<Job>> m_queue;
|
||||||
|
std::vector<std::thread> m_pool;
|
||||||
|
std::mutex m_poolMutex;
|
||||||
|
std::condition_variable m_wakeCondition;
|
||||||
|
|
||||||
|
bool m_isStopped;
|
||||||
|
std::atomic<int> m_numberRunning;
|
||||||
|
std::atomic<int> m_queueSize;
|
||||||
|
std::atomic<bool> m_initShutdown;
|
||||||
|
};
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
Loading…
Reference in New Issue
Block a user