From 19ccc3610db9efc63ef36b4ab5ceaeb8471684d7 Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Thu, 26 May 2022 14:12:50 -0400 Subject: [PATCH] Fleshed out project. Needs testing --- src/CalDict/cal_dict.cxx_tmp_3535 | 36 +++++++ src/CalDict/cal_dict.cxx_tmp_3544 | 36 +++++++ src/CalDict/cal_dict.cxx_tmp_3665 | 36 +++++++ src/CalibrationMap.cpp | 52 ++++++++++ src/CalibrationMap.h | 29 ++++++ src/Calibrator.cpp | 160 +++++++++++++++++++++++++++++ src/Calibrator.h | 36 +++++++ src/DataOrganizer.cpp | 163 ++++++++++++++++++++++++++++++ src/DataOrganizer.h | 58 +++++++++++ src/GainMap.cpp | 51 ++++++++++ src/GainMap.h | 30 ++++++ src/GainMatcher.cpp | 8 +- src/main.cpp | 136 +++++++++++++++++++++++++ 13 files changed, 827 insertions(+), 4 deletions(-) create mode 100644 src/CalDict/cal_dict.cxx_tmp_3535 create mode 100644 src/CalDict/cal_dict.cxx_tmp_3544 create mode 100644 src/CalDict/cal_dict.cxx_tmp_3665 create mode 100644 src/CalibrationMap.cpp create mode 100644 src/CalibrationMap.h create mode 100644 src/Calibrator.cpp create mode 100644 src/Calibrator.h create mode 100644 src/DataOrganizer.cpp create mode 100644 src/DataOrganizer.h create mode 100644 src/GainMap.cpp create mode 100644 src/GainMap.h create mode 100644 src/main.cpp diff --git a/src/CalDict/cal_dict.cxx_tmp_3535 b/src/CalDict/cal_dict.cxx_tmp_3535 new file mode 100644 index 0000000..5acc222 --- /dev/null +++ b/src/CalDict/cal_dict.cxx_tmp_3535 @@ -0,0 +1,36 @@ +// Do NOT change. Changes will be lost next time file is generated + +#define R__DICTIONARY_FILENAME srcdICalDictdIcal_dict +#define R__NO_DEPRECATION + +/*******************************************************************/ +#include +#include +#include +#include +#include +#define G__DICTIONARY +#include "RConfig.h" +#include "TClass.h" +#include "TDictAttributeMap.h" +#include "TInterpreter.h" +#include "TROOT.h" +#include "TBuffer.h" +#include "TMemberInspector.h" +#include "TInterpreter.h" +#include "TVirtualMutex.h" +#include "TError.h" + +#ifndef G__ROOT +#define G__ROOT +#endif + +#include "RtypesImp.h" +#include "TIsAProxy.h" +#include "TFileMergeInfo.h" +#include +#include "TCollectionProxyInfo.h" +/*******************************************************************/ + +#include "TDataMember.h" + diff --git a/src/CalDict/cal_dict.cxx_tmp_3544 b/src/CalDict/cal_dict.cxx_tmp_3544 new file mode 100644 index 0000000..5acc222 --- /dev/null +++ b/src/CalDict/cal_dict.cxx_tmp_3544 @@ -0,0 +1,36 @@ +// Do NOT change. Changes will be lost next time file is generated + +#define R__DICTIONARY_FILENAME srcdICalDictdIcal_dict +#define R__NO_DEPRECATION + +/*******************************************************************/ +#include +#include +#include +#include +#include +#define G__DICTIONARY +#include "RConfig.h" +#include "TClass.h" +#include "TDictAttributeMap.h" +#include "TInterpreter.h" +#include "TROOT.h" +#include "TBuffer.h" +#include "TMemberInspector.h" +#include "TInterpreter.h" +#include "TVirtualMutex.h" +#include "TError.h" + +#ifndef G__ROOT +#define G__ROOT +#endif + +#include "RtypesImp.h" +#include "TIsAProxy.h" +#include "TFileMergeInfo.h" +#include +#include "TCollectionProxyInfo.h" +/*******************************************************************/ + +#include "TDataMember.h" + diff --git a/src/CalDict/cal_dict.cxx_tmp_3665 b/src/CalDict/cal_dict.cxx_tmp_3665 new file mode 100644 index 0000000..5acc222 --- /dev/null +++ b/src/CalDict/cal_dict.cxx_tmp_3665 @@ -0,0 +1,36 @@ +// Do NOT change. Changes will be lost next time file is generated + +#define R__DICTIONARY_FILENAME srcdICalDictdIcal_dict +#define R__NO_DEPRECATION + +/*******************************************************************/ +#include +#include +#include +#include +#include +#define G__DICTIONARY +#include "RConfig.h" +#include "TClass.h" +#include "TDictAttributeMap.h" +#include "TInterpreter.h" +#include "TROOT.h" +#include "TBuffer.h" +#include "TMemberInspector.h" +#include "TInterpreter.h" +#include "TVirtualMutex.h" +#include "TError.h" + +#ifndef G__ROOT +#define G__ROOT +#endif + +#include "RtypesImp.h" +#include "TIsAProxy.h" +#include "TFileMergeInfo.h" +#include +#include "TCollectionProxyInfo.h" +/*******************************************************************/ + +#include "TDataMember.h" + diff --git a/src/CalibrationMap.cpp b/src/CalibrationMap.cpp new file mode 100644 index 0000000..0f56766 --- /dev/null +++ b/src/CalibrationMap.cpp @@ -0,0 +1,52 @@ +#include "CalibrationMap.h" +#include + +namespace SabreCal { + + CalibrationMap::CalibrationMap() : + m_isValid(false) + { + } + + CalibrationMap::CalibrationMap(const std::string& filename) : + m_isValid(false) + { + Init(filename); + } + + CalibrationMap::~CalibrationMap() {} + + void CalibrationMap::Init(const std::string &filename) + { + std::ifstream input(filename); + if(!input.is_open()) + { + m_isValid = false; + return; + } + + std::string junk; + int gchan; + Parameters params; + + std::getline(input, junk); + std::getline(input, junk); + + while(input>>gchan) + { + input>>params.slope>>params.offset; + m_map[gchan] = params; + } + + m_isValid = true; + } + + const Parameters& CalibrationMap::GetParameters(int gchan) const + { + auto iter = m_map.find(gchan); + if(iter != m_map.end()) + return iter->second; + else + return m_dummyParams; + } +} diff --git a/src/CalibrationMap.h b/src/CalibrationMap.h new file mode 100644 index 0000000..3a09aab --- /dev/null +++ b/src/CalibrationMap.h @@ -0,0 +1,29 @@ +#ifndef CALIBRATION_MAP_H +#define CALIBRATION_MAP_H + +#include +#include +#include "CalibrationStructs.h" + +namespace SabreCal { + + class CalibrationMap + { + public: + CalibrationMap(); + CalibrationMap(const std::string& filename); + ~CalibrationMap(); + + void Init(const std::string& filename); + + const Parameters& GetParameters(int gchan) const; + inline const bool IsValid() const { return m_isValid; } + + private: + bool m_isValid; + std::unordered_map m_map; + Parameters m_dummyParams; //for null result + }; +} + +#endif diff --git a/src/Calibrator.cpp b/src/Calibrator.cpp new file mode 100644 index 0000000..945094c --- /dev/null +++ b/src/Calibrator.cpp @@ -0,0 +1,160 @@ +#include "Calibrator.h" +#include "TFile.h" +#include "TTree.h" +#include "DataStructs.h" +#include +#include +#include "TSpectrum.h" + +namespace SabreCal { + + Calibrator::Calibrator(const std::string& gainfile) : + m_gainMap(gainfile) + { + TH1::AddDirectory(kFALSE); + } + + Calibrator::~Calibrator() {} + + void Calibrator::Run(const std::string& datafile, const std::string& outputfile, const std::string& plotfile) + { + if(!m_gainMap.IsValid()) + { + std::cerr<<"ERR -- Bad gain map at Calibrator::Run!"<IsOpen()) + { + std::cerr<<"ERR -- Bad data file "<Get("Data"); + uint16_t energy, board, channel; + tree->SetBranchAddress("Energy", &energy); + tree->SetBranchAddress("Board", &board); + tree->SetBranchAddress("Channel", &channel); + + int gchan; + + uint64_t nentries = tree->GetEntries(); + float flush_percent = 0.1f; + uint64_t count=0, flush_count=0, flush_val=nentries*flush_percent; + + std::vector histograms; + for(int i=0; iGetEntry(i); + count++; + if(count == flush_val) + { + flush_count++; + count = 0; + std::cout<<"\rPercent of data processed: "< 100.0) //Reject electronic noise + histograms[i].Fill(params.slope*(energy*s_gainFactor) + params.offset); + } + std::cout<Close(); + + std::cout<<"Finding calibration parameters..."<& grams) + { + TSpectrum finder; + double threshold = 0.5; + double sigma = 5.0; + + Parameters these_params; + double peakMean; + + for(auto& histo : grams) + { + finder.Search(&histo, sigma, "nobackground", threshold); + if(finder.GetNPeaks() == 0) + { + std::cerr<<"No peaks found in histogram "< 1) + { + peakMean = 0.0; + double peakAmp = 0.0; + std::cerr<<"Multiple peaks found in histogram "< peakAmp) + { + peakMean = mean; + peakAmp = amp; + } + } + these_params.slope = s_241Am_alphaKE/peakMean; + these_params.offset = 0.0; + m_params.push_back(these_params); + } + else + { + peakMean = finder.GetPositionX()[0]; + these_params.slope = s_241Am_alphaKE/peakMean; + these_params.offset = 0.0; + m_params.push_back(these_params); + } + } + } + + void Calibrator::WriteHistograms(const std::vector &grams, const std::string &plotfile) + { + TFile* output = TFile::Open(plotfile.c_str(), "RECREATE"); + if(!output->IsOpen()) + { + std::cerr<<"Unable to open file "<Close(); + } + + void Calibrator::WriteParameters(const std::string &outputfile) + { + std::ofstream output(outputfile); + if(!output.is_open()) + { + std::cerr<<"Unable to open file "< +#include +#include "ChannelMap.h" +#include "GainMap.h" +#include "CalibrationStructs.h" + +#include "TH1.h" + +namespace SabreCal { + + class Calibrator + { + public: + Calibrator(const std::string& gainfile); + ~Calibrator(); + + void Run(const std::string& datafile, const std::string& outputfile, const std::string& plotfile=""); + + private: + void FindParameters(const std::vector& grams); + void WriteHistograms(const std::vector& grams, const std::string& plotfile); + void WriteParameters(const std::string& outputfile); + + GainMap m_gainMap; + std::vector m_params; + + static constexpr int s_totalChannels = 127; //Total channels in SABRE + static constexpr double s_241Am_alphaKE = 5.468; //MeV + static constexpr double s_gainFactor = 1.5; //For if you're dumb and the electronics gain needs adjusted + }; +} + +#endif /* Calibrator_hpp */ diff --git a/src/DataOrganizer.cpp b/src/DataOrganizer.cpp new file mode 100644 index 0000000..ed0c3f2 --- /dev/null +++ b/src/DataOrganizer.cpp @@ -0,0 +1,163 @@ +#include "DataOrganizer.h" +#include "TFile.h" +#include "TTree.h" + +#include + +namespace SabreCal { + + DataOrganizer::DataOrganizer(const std::string& channelfile, const std::string& gainfile, + const std::string& calfile) : + m_channelMap(channelfile), m_gainMap(gainfile), m_calMap(calfile) + { + } + + DataOrganizer::~DataOrganizer() {} + + void DataOrganizer::Run(const std::string& inputdata, const std::string& outputdata) + { + if(!(m_channelMap.IsValid() && m_gainMap.IsValid() && m_calMap.IsValid())) + { + std::cerr<<"ERR -- Invalid map files at DataOrganizer::Run()"<IsOpen()) + { + std::cerr<<"ERR -- Unable to open input data file "<Get("SPSTree"); + + TFile* output = TFile::Open(outputdata.c_str(), "RECREATE"); + if(!output->IsOpen()) + { + std::cerr<<"ERR -- Unable to open output data file "<SetBranchAddress("event", &inevent); + + CalEvent outevent; + outtree->Branch("event", &outevent); + + std::vector rings, wedges; + double eTemp; + double ratio; + SabrePair this_pair; + + uint64_t nevents = intree->GetEntries(); + m_stats.totalEvents = nevents; + double flush_percent = 0.05; + uint64_t count=0, flush_count=0, flush_val=nevents*flush_percent; + + for(uint64_t i=0; iGetEvent(i); + count++; + if(count == flush_val) + { + flush_count++; + count = 0; + std::cout<<"\rPercent of data processed: "<xavg; + outevent.scintE = inevent->scintLeft; + outevent.cathodeE = inevent->cathode; + outevent.anodeFrontE = inevent->anodeFront; + outevent.anodeBackE = inevent->anodeBack; + outevent.scintT = inevent->scintLeftTime; + + for(int j=0; j<5; j++) + { + rings.clear(); + wedges.clear(); + auto& this_sabre = inevent->sabreArray[j]; + + int64_t diff = this_sabre.rings.size() - this_sabre.wedges.size(); + m_stats.sabreRingFrags += diff > 0 ? diff : 0; + m_stats.sabreWedgeFrags += diff < 0 ? -1*diff : 0; + for(auto& ring : this_sabre.rings) + { + m_stats.totalSabreHits++; + auto& ringGains = m_gainMap.GetParameters(ring.Ch); + auto& ringCals = m_calMap.GetParameters(ring.Ch); + eTemp = ringCals.slope*(ringGains.slope*ring.Long + ringGains.offset) + ringCals.offset; + if(eTemp > s_sabreThreshold) + { + m_stats.goodSabreHits++; + rings.emplace_back(eTemp, ring.Time, ring.Ch, false); + } + else + m_stats.sabreBelowThreshold++; + } + for(auto& wedge : this_sabre.wedges) + { + m_stats.totalSabreHits++; + auto& wedgeGains = m_gainMap.GetParameters(wedge.Ch); + auto& wedgeCals = m_calMap.GetParameters(wedge.Ch); + eTemp = wedgeCals.slope*(wedgeGains.slope*wedge.Long + wedgeGains.offset) + wedgeCals.offset; + if(eTemp > s_sabreThreshold) + { + m_stats.goodSabreHits++; + wedges.emplace_back(eTemp, wedge.Time, wedge.Ch, false); + } + else + m_stats.sabreBelowThreshold++; + } + + //Now match rings and wedges + for(auto& ring : rings) + { + for(auto& wedge : wedges) + { + if(wedge.used || ring.used) + continue; + ratio = std::abs(1.0 - ring.energy/wedge.energy); + if(ratio < s_sabreMatchCond) + { + auto& ringdata = m_channelMap.GetChannel(ring.gchan); + auto& wedgedata = m_channelMap.GetChannel(wedge.gchan); + this_pair.ringch = ring.gchan; + this_pair.wedgech = wedge.gchan; + this_pair.local_ring = ringdata.localChannel; + this_pair.local_wedge = wedgedata.localChannel; + this_pair.detID = wedgedata.detID; + this_pair.ringE = ring.energy; + this_pair.wedgeE = wedge.energy; + outevent.sabre.push_back(this_pair); + m_stats.pairsMade++; + break; + } + } + } + } + outtree->Fill(); + } + std::cout<Close(); + output->cd(); + outtree->Write(outtree->GetName(), TObject::kOverwrite); + output->Close(); + } + + void DataOrganizer::ShowStats() + { + std::cout<<"Total number of events: "< + +namespace SabreCal { + + struct AnalysisStats + { + uint64_t totalEvents = 0; + uint64_t totalSabreHits = 0; + uint64_t goodSabreHits = 0; //Above threshold, good number for matching + uint64_t pairsMade = 0; + uint64_t sabreRingFrags = 0; + uint64_t sabreWedgeFrags = 0; + uint64_t sabreBelowThreshold = 0; + }; + + struct SabreTemp + { + SabreTemp(double e, double t, int g, bool u) : + energy(e), time(t), gchan(g), used(u) + {} + double energy; + double time; + int gchan; + bool used = false; + }; + + class DataOrganizer + { + public: + DataOrganizer(const std::string& channelfile, const std::string& gainfile, const std::string& calfile); + ~DataOrganizer(); + + void Run(const std::string& inputdata, const std::string& outputdata); + void ShowStats(); + + private: + ChannelMap m_channelMap; + GainMap m_gainMap; + CalibrationMap m_calMap; + + AnalysisStats m_stats; + + static constexpr double s_sabreMatchCond = 0.2; //Ring/wedge %diff + static constexpr double s_sabreThreshold = 0.2; //MeV + static constexpr double s_sabreTimeCond = 100.0; //ns + }; +} + +#endif diff --git a/src/GainMap.cpp b/src/GainMap.cpp new file mode 100644 index 0000000..2a8cdf6 --- /dev/null +++ b/src/GainMap.cpp @@ -0,0 +1,51 @@ +#include "GainMap.h" +#include + +namespace SabreCal { + + GainMap::GainMap() : + m_isValid(false) + { + } + + GainMap::GainMap(const std::string& filename) : + m_isValid(false) + { + Init(filename); + } + + GainMap::~GainMap() {} + + void GainMap::Init(const std::string& filename) + { + std::ifstream input(filename); + if(!input.is_open()) + { + m_isValid = false; + return; + } + std::string junk; + int gchan; + Parameters pars; + + std::getline(input, junk); + std::getline(input, junk); + + while(input>>gchan) + { + input>>pars.slope>>pars.offset; + m_map[gchan] = pars; + } + + m_isValid = true; + } + + const Parameters& GainMap::GetParameters(int gchan) const + { + auto iter = m_map.find(gchan); + if(iter != m_map.end()) + return iter->second; + else + return m_dummyParams; + } +} diff --git a/src/GainMap.h b/src/GainMap.h new file mode 100644 index 0000000..328bbd0 --- /dev/null +++ b/src/GainMap.h @@ -0,0 +1,30 @@ +#ifndef GAIN_MAP_H +#define GAIN_MAP_H + +#include +#include +#include "CalibrationStructs.h" + +namespace SabreCal { + + class GainMap + { + public: + GainMap(); + GainMap(const std::string& filename); + ~GainMap(); + + void Init(const std::string& filename); + + const Parameters& GetParameters(int gchan) const; + inline const bool IsValid() const { return m_isValid; } + + private: + std::unordered_map m_map; + Parameters m_dummyParams; //for null results + + bool m_isValid; + }; +} + +#endif diff --git a/src/GainMatcher.cpp b/src/GainMatcher.cpp index 8e1bed0..065271f 100644 --- a/src/GainMatcher.cpp +++ b/src/GainMatcher.cpp @@ -50,7 +50,7 @@ namespace SabreCal { { count = 0; flushes++; - std::cout<<"\rPercent of data processed: "< +#include +#include + +int main(int argc, const char** argv) +{ + std::string inputfile; + std::string option; + if(argc < 2) + { + std::cerr<<"ERR -- SabreCal requires an input configuration file."< "<>junk>>gaindata; + input>>junk>>gainplots; + input>>junk>>caldata; + input>>junk>>calplots; + input>>junk>>inputdata; + input>>junk>>outputdata; + input>>junk>>gainfile; + input>>junk>>calfile; + input>>junk>>channelfile; + input>>junk>>ring_to_match>>junk>>wedge_to_match; + + std::cout<<"----------SABRE Calibration----------"< "<