gwm17/Mask
1
0
Fork 0
mirror of https://github.com/gwm17/Mask.git synced 2024-11-22 10:18:50 -05:00
Code Issues Projects Releases Wiki Activity

Implemented multithreading in detectors, overhaul application api for detectors, implement detector input file. Need robust testing to verify performance gain in detectors. Verified that code compiles and executes.

Browse Source
This commit is contained in:
Gordon McCann 2022-09-07 11:29:57 -04:00
parent 0456f230f5
commit 2f7b0969dc
19 changed files with 476 additions and 538 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
detector_input.txt Normal file
View 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

45
kinematics_input.txt Normal file
View File

@ -0,0 +1,45 @@
----------Data Information----------
OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
NumberOfThreads: 6
----------Reaction Information----------
NumberOfSamples: 1000000
begin_chain
begin_step
Type: Reaction
begin_nuclei
5 10
2 3
2 4
end_nuclei
BeamEnergyMean(MeV): 24.0
BeamEnergySigma(MeV): 0.0
ThetaType: Lab
ThetaMin(deg): 15.0
ThetaMax(deg): 15.0
PhiMin(deg): 0.0
PhMax(deg): 0.0
ResidualExcitationMean(MeV): 16.8
ResidualExcitationSigma(MeV): 0.023
end_step
begin_step
Type: Decay
begin_nuclei
5 9
2 4
end_nuclei
PhiMin(deg): 0.0
PhMax(deg): 360.0
ResidualExcitationMean(MeV): 0.0
ResidualExcitationSigma(MeV): 0.0
AngularDistributionFile: ./etc/isotropic_dist.txt
end_step
end_chain
----------Target Information----------
NumberOfLayers: 1
begin_layer
Thickness(ug/cm^2): 74
begin_elements (Z, A, Stoich.)
element 5 10 1
end_elements
end_layer
--------------------------------------

240
src/Detectors/AnasenArray.cpp
View File

@ -230,7 +230,7 @@ double AnasenArray::RunConsistencyCheck()
}
DetectorResult AnasenArray::IsRing1(Mask::Nucleus& nucleus)
DetectorResult AnasenArray::IsRing1(const Mask::Nucleus& nucleus)
{
DetectorResult observation;
double thetaIncident;
@ -255,7 +255,7 @@ DetectorResult AnasenArray::IsRing1(Mask::Nucleus& nucleus)
return observation;
}
DetectorResult AnasenArray::IsRing2(Mask::Nucleus& nucleus)
DetectorResult AnasenArray::IsRing2(const Mask::Nucleus& nucleus)
{
DetectorResult observation;
double thetaIncident;
@ -280,7 +280,7 @@ DetectorResult AnasenArray::IsRing2(Mask::Nucleus& nucleus)
return observation;
}
DetectorResult AnasenArray::IsQQQ(Mask::Nucleus& nucleus)
DetectorResult AnasenArray::IsQQQ(const Mask::Nucleus& nucleus)
{
DetectorResult observation;
double thetaIncident;
@ -324,7 +324,7 @@ DetectorResult AnasenArray::IsQQQ(Mask::Nucleus& nucleus)
return observation;
}
DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus)
DetectorResult AnasenArray::IsDetected(const Mask::Nucleus& nucleus)
{
DetectorResult result;
if(nucleus.GetKE() <= s_energyThreshold)
@ -337,236 +337,4 @@ DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus)
if(!result.detectFlag)
result = IsQQQ(nucleus);
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;
}

16
src/Detectors/AnasenArray.h
View File

@ -6,8 +6,8 @@
#include "DetectorArray.h"
#include "SX3Detector.h"
#include "QQQDetector.h"
#include "Target.h"
#include "Nucleus.h"
#include "Mask/Target.h"
#include "Mask/Nucleus.h"
#include "AnasenDeadChannelMap.h"
class AnasenArray : public DetectorArray
@ -15,17 +15,15 @@ class AnasenArray : public DetectorArray
public:
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 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:
DetectorResult IsRing1(Mask::Nucleus& nucleus);
DetectorResult IsRing2(Mask::Nucleus& nucleus);
DetectorResult IsQQQ(Mask::Nucleus& nucleus);
DetectorResult IsAnasen(Mask::Nucleus& nucleus);
void CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts);
DetectorResult IsRing1(const Mask::Nucleus& nucleus);
DetectorResult IsRing2(const Mask::Nucleus& nucleus);
DetectorResult IsQQQ(const Mask::Nucleus& nucleus);
std::vector<SX3Detector> m_Ring1;
std::vector<SX3Detector> m_Ring2;

3
src/Detectors/CMakeLists.txt
View File

@ -18,6 +18,9 @@ target_sources(Detectors PUBLIC
SabreArray.h
SX3Detector.cpp
SX3Detector.h
DetectorApp.h
DetectorApp.cpp
DetectorArray.cpp
)
target_link_libraries(Detectors

144
src/Detectors/DetectorApp.cpp Normal file
View 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
View 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
View 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;
}

17
src/Detectors/DetectorArray.h
View File

@ -5,6 +5,7 @@
#include <cmath>
#include "Math/Point3D.h"
#include "Mask/Nucleus.h"
struct DetectorResult
{
@ -14,15 +15,23 @@ struct DetectorResult
std::string det_name = "";
};
enum class ArrayType
{
None,
Anasen,
Sabre
};
class DetectorArray
{
public:
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 double RunConsistencyCheck() = 0;
virtual void SetDeadChannelMap(const std::string& filename) = 0;
protected:
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;
};
DetectorArray* CreateDetectorArray(ArrayType type);
std::string ArrayTypeToString(ArrayType type);
ArrayType StringToArrayType(const std::string& value);
#endif

237
src/Detectors/SabreArray.cpp
View File

@ -31,161 +31,6 @@ 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)
{
std::ofstream output(filename);
@ -259,7 +104,7 @@ double SabreArray::RunConsistencyCheck()
}
/*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;
if(nucleus.GetKE() <= s_energyThreshold)
@ -297,82 +142,4 @@ DetectorResult SabreArray::IsSabre(Mask::Nucleus& nucleus)
observation.detectFlag = false;
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]++;
}
}
}
}

8
src/Detectors/SabreArray.h
View File

@ -3,7 +3,7 @@
#include "DetectorArray.h"
#include "SabreDetector.h"
#include "Target.h"
#include "Mask/Target.h"
#include "SabreDeadChannelMap.h"
#include "Mask/Nucleus.h"
@ -12,14 +12,12 @@ class SabreArray : public DetectorArray
public:
SabreArray();
~SabreArray();
void SetDeadChannelMap(const std::string& filename) { m_deadMap.LoadMapfile(filename); };
void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
virtual void SetDeadChannelMap(const std::string& filename) override { m_deadMap.LoadMapfile(filename); };
virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus) override;
void DrawDetectorSystem(const std::string& filename) override;
double RunConsistencyCheck() override;
private:
DetectorResult IsSabre(Mask::Nucleus& nucleus);
void CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts);
std::vector<SabreDetector> m_detectors;

40
src/Detectors/main.cpp
View File

@ -1,5 +1,4 @@
#include "SabreArray.h"
#include "AnasenArray.h"
#include "DetectorApp.h"
#include "KinematicsExceptions.h"
#include <iostream>
#include <string>
@ -7,7 +6,7 @@
int main(int argc, char** argv)
{
if(argc != 4)
if(argc != 2)
{
std::cerr<<"Incorrect number of commandline arguments! Returning."<<std::endl;
return 1;
@ -19,26 +18,21 @@ int main(int argc, char** argv)
return 1;
}
std::string inputname = argv[1];
std::string outputname = argv[2];
std::string statsname = argv[3];
try
{
DetectorApp app;
if(!app.LoadConfig(argv[1]))
{
std::cerr << "Unable to load config file " << argv[1] << ". Shutting down." << std::endl;
return 1;
}
SabreArray sabre;
std::string mapfile = "./etc/sabreDeadChannels_May2022.txt";
sabre.SetDeadChannelMap(mapfile);
sabre.CalculateEfficiency(inputname, outputname, statsname);
//std::cout<<"Running consistency check(0=success): "<<sabre.RunConsistencyCheck()<<std::endl;
//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");
*/
app.Run();
}
catch(std::exception& e)
{
std::cerr << e.what() << std::endl;
return 1;
}
return 0;
}

2
src/Mask/CMakeLists.txt
View File

@ -53,6 +53,8 @@ target_sources(Mask PRIVATE
ThreadPool.h
FileWriter.h
FileWriter.cpp
FileReader.h
FileReader.cpp
)
set(THREADS_PREFER_PTHREAD_FLAG On)

66
src/Mask/FileReader.cpp Normal file
View 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
View 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

6
src/Mask/FileWriter.h
View File

@ -19,15 +19,15 @@ namespace Mask {
FileWriter(const std::string& filename, const std::string& treename);
~FileWriter();
bool IsOpen() const { return m_file->IsOpen(); }
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 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);
void Open(const std::string& filename, const std::string& treename); //Not thread safe!
void Close(); //Not thread safe!
private:

42
src/Mask/MaskApp.cpp
View File

@ -39,9 +39,16 @@ namespace Mask {
input>>junk>>m_nsamples;
std::cout<<"Allocating resources... Asking for " << m_nthreads << " threads...";
m_resources = std::make_unique<ThreadPool>(m_nthreads);
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");
@ -174,12 +181,27 @@ namespace Mask {
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(auto system : m_systemList)
m_resources->PushJob({std::bind(&MaskApp::RunChunk, std::ref(*this), std::placeholders::_1), system});
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;
@ -243,18 +265,4 @@ namespace Mask {
std::cout<<"---------------------------------------------"<<std::endl;
}
void MaskApp::RunChunk(ReactionSystem* system)
{
if(system == nullptr)
return;
uint64_t samples = m_nsamples / m_nthreads;
for(uint64_t i=0; i<samples; i++)
{
system->RunSystem();
m_fileWriter.PushData(*(system->GetNuclei()));
}
}
}

3
src/Mask/MaskApp.h
View File

@ -36,8 +36,9 @@ namespace Mask {
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> m_resources;
std::unique_ptr<ThreadPool<ReactionSystem*, uint64_t>> m_resources;
};

19
src/Mask/ThreadPool.h
View File

@ -20,16 +20,19 @@
namespace Mask {
using JobFunction = std::function<void(ReactionSystem*)>;
struct Job
{
JobFunction func;
ReactionSystem* argument;
};
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)
{
@ -94,7 +97,7 @@ namespace Mask {
m_numberRunning++;
m_queueSize--;
job.func(job.argument);
std::apply(job.func, job.argument);
m_numberRunning--;
}