mirror of
https://github.com/gwm17/Mask.git
synced 2024-11-22 10:18:50 -05:00
Gordon McCann
GitHub
commit
25fee0f39b
|
|
@ -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 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
|
||||
|
||||
`./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
|
||||
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----------
|
||||
OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
|
||||
NumberOfThreads: 6
|
||||
----------Reaction Information----------
|
||||
NumberOfSamples: 100000
|
||||
NumberOfSamples: 1000000
|
||||
begin_chain
|
||||
begin_step
|
||||
Type: Reaction
|
||||
|
|
@ -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)
|
||||
|
|
@ -338,235 +338,3 @@ DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus)
|
|||
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;
|
||||
}
|
||||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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
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 "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
|
||||
|
|
@ -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)
|
||||
|
|
@ -298,81 +143,3 @@ 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]++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -24,6 +24,7 @@ int main(int argc, char** argv)
|
|||
return 1;
|
||||
}
|
||||
calculator.Run();
|
||||
//calculator.RunSingleThread();
|
||||
}
|
||||
catch(const std::exception& e)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -50,7 +50,15 @@ target_sources(Mask PRIVATE
|
|||
ThreeStepSystem.h
|
||||
TwoStepSystem.cpp
|
||||
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})
|
||||
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 {
|
||||
|
||||
MaskApp::MaskApp() :
|
||||
m_system(nullptr)
|
||||
m_system(nullptr), m_resources(nullptr)
|
||||
{
|
||||
std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
|
||||
}
|
||||
|
|
@ -16,6 +16,8 @@ namespace Mask {
|
|||
MaskApp::~MaskApp()
|
||||
{
|
||||
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)
|
||||
|
|
@ -31,10 +33,25 @@ namespace Mask {
|
|||
std::string junk;
|
||||
std::getline(input, junk);
|
||||
input>>junk>>m_outputName;
|
||||
input>>junk>>m_nthreads;
|
||||
std::getline(input, junk);
|
||||
std::getline(input, junk);
|
||||
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;
|
||||
int z, a;
|
||||
while(input>>junk)
|
||||
|
|
@ -98,10 +115,20 @@ namespace Mask {
|
|||
m_system = CreateSystem(params);
|
||||
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;
|
||||
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);
|
||||
|
||||
|
|
@ -132,9 +159,13 @@ namespace Mask {
|
|||
}
|
||||
m_system->SetLayeredTarget(target);
|
||||
|
||||
std::cout<<"Outputing data to file: "<<m_outputName<<std::endl;
|
||||
std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
|
||||
std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
|
||||
for(auto system : m_systemList)
|
||||
{
|
||||
system->SetLayeredTarget(target);
|
||||
}
|
||||
|
||||
std::cout<<"Reaction equation: "<<m_system->GetSystemEquation()<<std::endl;
|
||||
std::cout<<"Number of samples: "<<m_nsamples<<std::endl;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
|
@ -146,6 +177,57 @@ namespace Mask {
|
|||
}
|
||||
|
||||
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;
|
||||
if(m_system == nullptr)
|
||||
|
|
@ -158,11 +240,11 @@ namespace Mask {
|
|||
tree->Branch("nuclei", m_system->GetNuclei());
|
||||
|
||||
//For progress tracking
|
||||
uint32_t percent5 = 0.05*m_nsamples;
|
||||
uint32_t count = 0;
|
||||
uint32_t npercent = 0;
|
||||
uint64_t percent5 = 0.05*m_nsamples;
|
||||
uint64_t count = 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)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -7,6 +7,10 @@
|
|||
#include "TwoStepSystem.h"
|
||||
#include "ThreeStepSystem.h"
|
||||
#include "RxnType.h"
|
||||
#include "ThreadPool.h"
|
||||
#include "FileWriter.h"
|
||||
|
||||
#include <memory>
|
||||
|
||||
namespace Mask {
|
||||
|
||||
|
|
@ -17,18 +21,24 @@ namespace Mask {
|
|||
~MaskApp();
|
||||
bool LoadConfig(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 RunSingleThread();
|
||||
|
||||
private:
|
||||
void RunChunk(ReactionSystem* system);
|
||||
|
||||
ReactionSystem* m_system;
|
||||
|
||||
std::string m_outputName;
|
||||
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,7 +1,6 @@
|
|||
#include "RandomGenerator.h"
|
||||
|
||||
namespace Mask {
|
||||
RandomGenerator* RandomGenerator::s_instance = new RandomGenerator();
|
||||
|
||||
RandomGenerator::RandomGenerator()
|
||||
{
|
||||
|
|
|
|||
|
|
@ -10,12 +10,15 @@ namespace Mask {
|
|||
public:
|
||||
~RandomGenerator();
|
||||
std::mt19937& GetGenerator() { return rng; }
|
||||
static RandomGenerator& GetInstance() { return *s_instance; }
|
||||
static RandomGenerator& GetInstance()
|
||||
{
|
||||
static thread_local RandomGenerator s_instance;
|
||||
return s_instance;
|
||||
}
|
||||
|
||||
private:
|
||||
RandomGenerator();
|
||||
|
||||
static RandomGenerator* s_instance;
|
||||
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