From c497707b8d71bc6d78f884c40c9d2f2bb0264f3b Mon Sep 17 00:00:00 2001 From: dominikpanek <3adominikpanek@gmail.com> Date: Wed, 31 Aug 2022 10:46:19 +0200 Subject: [PATCH 1/8] Updated version of tripple coincidences code, accoring to Pawel remarks --- EventAnalysis.cpp | 363 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 363 insertions(+) create mode 100644 EventAnalysis.cpp diff --git a/EventAnalysis.cpp b/EventAnalysis.cpp new file mode 100644 index 0000000..9f8fb72 --- /dev/null +++ b/EventAnalysis.cpp @@ -0,0 +1,363 @@ +#include +#include +#include +#include + +#include "EventAnalysis.h" +#include "Hits.h" + +using namespace std; + +//================================================================================ +// MATHEMATICAL AND ADDITIONAL FUNCTIONS +//================================================================================ + +// struct used for sorting vectors of hits (functor) + +namespace event_analysis{ + +struct EntityComp { + + string property; + EntityComp(string property) {this->property = property;} + bool operator()(const Hit& h1, const Hit& h2) const { + bool val = false; + if(property == "TIME") + val = h1.time < h2.time; + else if (property == "EDEP") + val = h1.edep < h2.edep; + else if (property == "EVENTID") + val = h1.eventID < h2.eventID; + return val; + } + +}; + +void sort_hits(vector &hits, string key) { + + sort(hits.begin(), hits.end(), EntityComp(key)); + +} + +Hit merge_hits(const std::vector &hits, const AveragingMethod winner = kCentroidWinnerEnergyWeightedFirstTime) { + const unsigned int nHits = hits.size(); + /// WK: temporary checks for debugs to be switched off + assert(nHits > 0); + if (nHits > 0) { + auto eID = hits[0].eventID; + auto vID = hits[0].volumeID; + assert(std::all_of(hits.begin(), hits.end(), [eID](const Hit& h)->bool { return h.eventID == eID; })); + assert(std::all_of(hits.begin(), hits.end(), [vID](const Hit& h)->bool { return h.volumeID == vID; })); + } + /// WK: end + Hit h; + h.eventID = hits[0].eventID; + h.volumeID = hits[0].volumeID; + for (unsigned int i=0; i times; + for (unsigned int i = 0; i < nHits; i++) + times.push_back(hits[i].time); + unsigned int min_index = std::distance( + times.begin(), std::min_element(times.begin(), times.end())); + h.time = hits[min_index].time; + for (unsigned int i = 0; i < nHits; i++) { + h.posX += hits[i].posX * hits[i].edep; + h.posY += hits[i].posY * hits[i].edep; + h.posZ += hits[i].posZ * hits[i].edep; + } + h.posX /= h.edep; + h.posY /= h.edep; + h.posZ /= h.edep; + break; + } + + case kEnergyWinner: { + std::vector energies; + // for (unsigned int i = 0; i> select_coincident_hits(const vector &hits, double compton_energy_threshold) +{ + std::vector selected_hits; + + int nb_above_noise_treshold = hits.size(); + for (unsigned int i=0; i compton_energy_threshold) selected_hits.push_back(hits[i]); + } + int nb_above_compton_threshold = selected_hits.size(); + return std::make_tuple(nb_above_noise_treshold, nb_above_compton_threshold, selected_hits); +} + +std::tuple> select_coincident_singles(const std::vector &hits, double compton_energy_threshold) { + + const string systemType = string(getenv("GOJA_SYSTEM_TYPE")); + + map> singles_tmp; + for (unsigned int i=0; i tmp; + tmp.push_back(hits[i]); + singles_tmp[ID] = tmp; + } else { + singles_tmp[ID].push_back(hits[i]); + } + } + + /// WK: Why it is always centroid here? + vector singles; + map>::iterator it_tmp = singles_tmp.begin(); + while(it_tmp != singles_tmp.end()) { + singles.push_back(merge_hits(it_tmp->second, kCentroidWinnerEnergyWeightedFirstTime)); + it_tmp++; + } + return select_coincident_hits(singles, compton_energy_threshold); +} + +std::tuple> check_gamma_type(const vector &hits, double prompt_energy_threshold) +{ + std::vector triple_hits; + + for (unsigned int i=0; i prompt_energy_threshold) { + triple_hits.push_back(hits[i]); + prompt_multiplicity++; + } + } + if(prompt_multiplicity == 1) { + bool isGood = true; + } + else { + bool isGood = false; + } + return std::make_tuple(isGood, triple_hits); + } + +} + +EventType verify_type_of_triple_coincidence(const Hit &h1, const Hit &h2, const Hit &h3) { + + EventType t = kUnspecified; + + if (h1.eventID==h2.eventID and h1.eventID==h3.eventID) { //true, phantom-scattered and detector-scattered + if (h1.nPhantomCompton==0 and h2.nPhantomCompton==0 and h3.nPhantomCompton==0) { + if (h1.nCrystalCompton==1 and h2.nCrystalCompton==1 and h3.nCrystalCompton==1) { //true + t = kTrue; + } + else { //detector-scattered + t = kDetectorScattered; + } + } + else { //phantom-scattered + t = kPhantomScattered; + } + } + else { //accidental + t = kAccidental; + } + + return t; + +} + + +//================================================================================ +// PRINTING +//================================================================================ + +void print_triple_coincidences(const std::vector& hits, PROMPT_E_TH) { + + assert(hits.size() ==3); + + std::vector triple_hits; + bool isGood; + std::tie(isGood, triple_hits) = check_gamma_type(hits, PROMPT_E_TH); + + if(isGood){ + + cout.setf(ios::fixed); + + Hit h1 = triple_hits[0]; + Hit h2 = triple_hits[1]; + Hit h3 = triple_hits[2]; + + cout.precision(2); + cout << h1.posX/10. << "\t" << h1.posY/10. << "\t" << h1.posZ/10. << "\t"; + cout.precision(1); + cout << h1.time << "\t"; + + cout.precision(2); + cout << h2.posX/10. << "\t" << h2.posY/10. << "\t" << h2.posZ/10. << "\t"; + cout.precision(1); + cout << h2.time << "\t"; + + cout << h1.volumeID << "\t"; + cout << h2.volumeID << "\t"; + + cout.precision(2); + cout << h1.edep << "\t"; + cout << h2.edep << "\t"; + + cout << verify_type_of_triple_coincidence(h1, h2, h3) << "\t"; + + cout.precision(2); + cout << h1.sourcePosX/10. << "\t" << h1.sourcePosY/10. << "\t" << h1.sourcePosZ/10. << "\t"; + cout << h2.sourcePosX/10. << "\t" << h2.sourcePosY/10. << "\t" << h2.sourcePosZ/10. << endl; + + cout.precision(2); + cout << h3.posX/10. << "\t" << h3.posY/10. << "\t" << h3.posZ/10. << "\t"; + cout.precision(1); + cout << h3.time << "\t"; + cout << h3.volumeID << "\t"; + cout.precision(2); + cout << h3.edep << "\t"; + cout.precision(2); + cout << h3.sourcePosX/10. << "\t" << h3.sourcePosY/10. << "\t" << h3.sourcePosZ/10. << "\t"; + +} + +} + +//================================================================================ +// MAIN ANALYSIS FUNCTION +//================================================================================ + +void analyze_event(vector &hits, bool hits_are_singles, bool triple_coincidence) +{ + + double COMPTON_E_TH = atof(getenv("GOJA_COMPTON_E_TH"))*1e3; + double PROMPT_E_TH = atof(getenv(“GOJA_PROMPT_E_TH”))*1e3; + int MAX_N = int(atof(getenv("GOJA_MAX_N"))); + int MAX_N0 = int(atof(getenv("GOJA_MAX_N0"))); + int N =0; + int N0 =0; + sort_hits(hits, "TIME"); + + std::vector selected_hits; + if (hits_are_singles) { + std::tie(N0, N, selected_hits) = select_coincident_singles(hits, COMPTON_E_TH); + } + else { + std::tie(N0, N, selected_hits) = select_coincident_hits(hits, COMPTON_E_TH); + } + + + if (triple_coincidences) { + if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); + } + else { + if (N==MAX_N and N0<=MAX_N0) print_triple_coincidences(selected_hits); + } + + if (DEBUG) { + cout.setf(ios::fixed); + cout.precision(1); + for (unsigned int i=0; i0) { + cout.precision(1); + cout << "\t" << (hits[i].time-hits[i-1].time); + } + if (i==(hits.size()-1)) { + cout.precision(1); + cout << "\t" << (hits[i].time-hits[0].time) << endl; + cout << N << "\t" << N0 << endl; + } + cout << endl; + } + } + +} + +void RunTests() +{ + const double epsilon = 0.000001; + Hit h1; + h1.edep = 10; + Hit h2; + h2.edep = 80; + Hit h3; + h2.edep = 50; + std::vector hits = {h1, h2, h3}; + for (auto& h : hits) + { + h.eventID = 10; + h.volumeID = 9; + } + auto mergedHit = merge_hits(hits, kEnergyWinner); + std::cout << mergedHit.edep << std::endl; + // WK: this test breaks I guess the energy of the merged hit is wrongly assigned + /// It should be the energy of the highest hit, so 80 + // assert(std::abs(mergedHit.edep - 80) < epsilon); + assert(mergedHit.eventID == 10); + assert(mergedHit.volumeID == 9); +} +} \ No newline at end of file From d2a572e2cd7d572c7909779656e4d169fee9a5f4 Mon Sep 17 00:00:00 2001 From: dominikpanek <3adominikpanek@gmail.com> Date: Wed, 31 Aug 2022 10:49:18 +0200 Subject: [PATCH 2/8] New verison moved to the proper place --- EventAnalysis.cpp | 363 ----------------------------------------- goja/EventAnalysis.cpp | 208 ++++++++++------------- 2 files changed, 84 insertions(+), 487 deletions(-) delete mode 100644 EventAnalysis.cpp mode change 100755 => 100644 goja/EventAnalysis.cpp diff --git a/EventAnalysis.cpp b/EventAnalysis.cpp deleted file mode 100644 index 9f8fb72..0000000 --- a/EventAnalysis.cpp +++ /dev/null @@ -1,363 +0,0 @@ -#include -#include -#include -#include - -#include "EventAnalysis.h" -#include "Hits.h" - -using namespace std; - -//================================================================================ -// MATHEMATICAL AND ADDITIONAL FUNCTIONS -//================================================================================ - -// struct used for sorting vectors of hits (functor) - -namespace event_analysis{ - -struct EntityComp { - - string property; - EntityComp(string property) {this->property = property;} - bool operator()(const Hit& h1, const Hit& h2) const { - bool val = false; - if(property == "TIME") - val = h1.time < h2.time; - else if (property == "EDEP") - val = h1.edep < h2.edep; - else if (property == "EVENTID") - val = h1.eventID < h2.eventID; - return val; - } - -}; - -void sort_hits(vector &hits, string key) { - - sort(hits.begin(), hits.end(), EntityComp(key)); - -} - -Hit merge_hits(const std::vector &hits, const AveragingMethod winner = kCentroidWinnerEnergyWeightedFirstTime) { - const unsigned int nHits = hits.size(); - /// WK: temporary checks for debugs to be switched off - assert(nHits > 0); - if (nHits > 0) { - auto eID = hits[0].eventID; - auto vID = hits[0].volumeID; - assert(std::all_of(hits.begin(), hits.end(), [eID](const Hit& h)->bool { return h.eventID == eID; })); - assert(std::all_of(hits.begin(), hits.end(), [vID](const Hit& h)->bool { return h.volumeID == vID; })); - } - /// WK: end - Hit h; - h.eventID = hits[0].eventID; - h.volumeID = hits[0].volumeID; - for (unsigned int i=0; i times; - for (unsigned int i = 0; i < nHits; i++) - times.push_back(hits[i].time); - unsigned int min_index = std::distance( - times.begin(), std::min_element(times.begin(), times.end())); - h.time = hits[min_index].time; - for (unsigned int i = 0; i < nHits; i++) { - h.posX += hits[i].posX * hits[i].edep; - h.posY += hits[i].posY * hits[i].edep; - h.posZ += hits[i].posZ * hits[i].edep; - } - h.posX /= h.edep; - h.posY /= h.edep; - h.posZ /= h.edep; - break; - } - - case kEnergyWinner: { - std::vector energies; - // for (unsigned int i = 0; i> select_coincident_hits(const vector &hits, double compton_energy_threshold) -{ - std::vector selected_hits; - - int nb_above_noise_treshold = hits.size(); - for (unsigned int i=0; i compton_energy_threshold) selected_hits.push_back(hits[i]); - } - int nb_above_compton_threshold = selected_hits.size(); - return std::make_tuple(nb_above_noise_treshold, nb_above_compton_threshold, selected_hits); -} - -std::tuple> select_coincident_singles(const std::vector &hits, double compton_energy_threshold) { - - const string systemType = string(getenv("GOJA_SYSTEM_TYPE")); - - map> singles_tmp; - for (unsigned int i=0; i tmp; - tmp.push_back(hits[i]); - singles_tmp[ID] = tmp; - } else { - singles_tmp[ID].push_back(hits[i]); - } - } - - /// WK: Why it is always centroid here? - vector singles; - map>::iterator it_tmp = singles_tmp.begin(); - while(it_tmp != singles_tmp.end()) { - singles.push_back(merge_hits(it_tmp->second, kCentroidWinnerEnergyWeightedFirstTime)); - it_tmp++; - } - return select_coincident_hits(singles, compton_energy_threshold); -} - -std::tuple> check_gamma_type(const vector &hits, double prompt_energy_threshold) -{ - std::vector triple_hits; - - for (unsigned int i=0; i prompt_energy_threshold) { - triple_hits.push_back(hits[i]); - prompt_multiplicity++; - } - } - if(prompt_multiplicity == 1) { - bool isGood = true; - } - else { - bool isGood = false; - } - return std::make_tuple(isGood, triple_hits); - } - -} - -EventType verify_type_of_triple_coincidence(const Hit &h1, const Hit &h2, const Hit &h3) { - - EventType t = kUnspecified; - - if (h1.eventID==h2.eventID and h1.eventID==h3.eventID) { //true, phantom-scattered and detector-scattered - if (h1.nPhantomCompton==0 and h2.nPhantomCompton==0 and h3.nPhantomCompton==0) { - if (h1.nCrystalCompton==1 and h2.nCrystalCompton==1 and h3.nCrystalCompton==1) { //true - t = kTrue; - } - else { //detector-scattered - t = kDetectorScattered; - } - } - else { //phantom-scattered - t = kPhantomScattered; - } - } - else { //accidental - t = kAccidental; - } - - return t; - -} - - -//================================================================================ -// PRINTING -//================================================================================ - -void print_triple_coincidences(const std::vector& hits, PROMPT_E_TH) { - - assert(hits.size() ==3); - - std::vector triple_hits; - bool isGood; - std::tie(isGood, triple_hits) = check_gamma_type(hits, PROMPT_E_TH); - - if(isGood){ - - cout.setf(ios::fixed); - - Hit h1 = triple_hits[0]; - Hit h2 = triple_hits[1]; - Hit h3 = triple_hits[2]; - - cout.precision(2); - cout << h1.posX/10. << "\t" << h1.posY/10. << "\t" << h1.posZ/10. << "\t"; - cout.precision(1); - cout << h1.time << "\t"; - - cout.precision(2); - cout << h2.posX/10. << "\t" << h2.posY/10. << "\t" << h2.posZ/10. << "\t"; - cout.precision(1); - cout << h2.time << "\t"; - - cout << h1.volumeID << "\t"; - cout << h2.volumeID << "\t"; - - cout.precision(2); - cout << h1.edep << "\t"; - cout << h2.edep << "\t"; - - cout << verify_type_of_triple_coincidence(h1, h2, h3) << "\t"; - - cout.precision(2); - cout << h1.sourcePosX/10. << "\t" << h1.sourcePosY/10. << "\t" << h1.sourcePosZ/10. << "\t"; - cout << h2.sourcePosX/10. << "\t" << h2.sourcePosY/10. << "\t" << h2.sourcePosZ/10. << endl; - - cout.precision(2); - cout << h3.posX/10. << "\t" << h3.posY/10. << "\t" << h3.posZ/10. << "\t"; - cout.precision(1); - cout << h3.time << "\t"; - cout << h3.volumeID << "\t"; - cout.precision(2); - cout << h3.edep << "\t"; - cout.precision(2); - cout << h3.sourcePosX/10. << "\t" << h3.sourcePosY/10. << "\t" << h3.sourcePosZ/10. << "\t"; - -} - -} - -//================================================================================ -// MAIN ANALYSIS FUNCTION -//================================================================================ - -void analyze_event(vector &hits, bool hits_are_singles, bool triple_coincidence) -{ - - double COMPTON_E_TH = atof(getenv("GOJA_COMPTON_E_TH"))*1e3; - double PROMPT_E_TH = atof(getenv(“GOJA_PROMPT_E_TH”))*1e3; - int MAX_N = int(atof(getenv("GOJA_MAX_N"))); - int MAX_N0 = int(atof(getenv("GOJA_MAX_N0"))); - int N =0; - int N0 =0; - sort_hits(hits, "TIME"); - - std::vector selected_hits; - if (hits_are_singles) { - std::tie(N0, N, selected_hits) = select_coincident_singles(hits, COMPTON_E_TH); - } - else { - std::tie(N0, N, selected_hits) = select_coincident_hits(hits, COMPTON_E_TH); - } - - - if (triple_coincidences) { - if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); - } - else { - if (N==MAX_N and N0<=MAX_N0) print_triple_coincidences(selected_hits); - } - - if (DEBUG) { - cout.setf(ios::fixed); - cout.precision(1); - for (unsigned int i=0; i0) { - cout.precision(1); - cout << "\t" << (hits[i].time-hits[i-1].time); - } - if (i==(hits.size()-1)) { - cout.precision(1); - cout << "\t" << (hits[i].time-hits[0].time) << endl; - cout << N << "\t" << N0 << endl; - } - cout << endl; - } - } - -} - -void RunTests() -{ - const double epsilon = 0.000001; - Hit h1; - h1.edep = 10; - Hit h2; - h2.edep = 80; - Hit h3; - h2.edep = 50; - std::vector hits = {h1, h2, h3}; - for (auto& h : hits) - { - h.eventID = 10; - h.volumeID = 9; - } - auto mergedHit = merge_hits(hits, kEnergyWinner); - std::cout << mergedHit.edep << std::endl; - // WK: this test breaks I guess the energy of the merged hit is wrongly assigned - /// It should be the energy of the highest hit, so 80 - // assert(std::abs(mergedHit.edep - 80) < epsilon); - assert(mergedHit.eventID == 10); - assert(mergedHit.volumeID == 9); -} -} \ No newline at end of file diff --git a/goja/EventAnalysis.cpp b/goja/EventAnalysis.cpp old mode 100755 new mode 100644 index 9f7a7ec..9f8fb72 --- a/goja/EventAnalysis.cpp +++ b/goja/EventAnalysis.cpp @@ -53,7 +53,6 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner = kCen Hit h; h.eventID = hits[0].eventID; h.volumeID = hits[0].volumeID; - // Energy of single4 == sum of energies of hits: for (unsigned int i=0; i &hits, const AveragingMethod winner = kCen case kEnergyWinner: { std::vector energies; - for (unsigned int i = 0; i &hits, const AveragingMethod winner = kCen return h; } + + std::tuple> select_coincident_hits(const vector &hits, double compton_energy_threshold) { std::vector selected_hits; @@ -142,8 +143,9 @@ std::tuple> select_coincident_hits(const vector return std::make_tuple(nb_above_noise_treshold, nb_above_compton_threshold, selected_hits); } -std::tuple> select_coincident_singles (const std::vector &hits, - double compton_energy_threshold, const std::string& systemType, const AveragingMethod averagingMethod = kCentroidWinnerEnergyWeightedFirstTime) { +std::tuple> select_coincident_singles(const std::vector &hits, double compton_energy_threshold) { + + const string systemType = string(getenv("GOJA_SYSTEM_TYPE")); map> singles_tmp; for (unsigned int i=0; i> select_coincident_singles (const std::vec ID = std::to_string(hits[i].volumeID); else if (systemType == "cylindricalPET") ID = std::to_string(hits[i].rsectorID) + '_' + std::to_string(hits[i].layerID); - if (singles_tmp.find(ID) == singles_tmp.end()) { + if (singles_tmp.find(ID) == singles_tmp.end() ) { vector tmp; tmp.push_back(hits[i]); singles_tmp[ID] = tmp; @@ -161,22 +163,48 @@ std::tuple> select_coincident_singles (const std::vec } } + /// WK: Why it is always centroid here? vector singles; map>::iterator it_tmp = singles_tmp.begin(); while(it_tmp != singles_tmp.end()) { - singles.push_back(merge_hits(it_tmp->second, averagingMethod)); + singles.push_back(merge_hits(it_tmp->second, kCentroidWinnerEnergyWeightedFirstTime)); it_tmp++; } return select_coincident_hits(singles, compton_energy_threshold); } -EventType verify_type_of_coincidence(const Hit &h1, const Hit &h2) { +std::tuple> check_gamma_type(const vector &hits, double prompt_energy_threshold) +{ + std::vector triple_hits; + + for (unsigned int i=0; i prompt_energy_threshold) { + triple_hits.push_back(hits[i]); + prompt_multiplicity++; + } + } + if(prompt_multiplicity == 1) { + bool isGood = true; + } + else { + bool isGood = false; + } + return std::make_tuple(isGood, triple_hits); + } + +} + +EventType verify_type_of_triple_coincidence(const Hit &h1, const Hit &h2, const Hit &h3) { EventType t = kUnspecified; - if (h1.eventID==h2.eventID) { //true, phantom-scattered and detector-scattered - if (h1.nPhantomCompton==0 and h2.nPhantomCompton==0) { - if (h1.nCrystalCompton==1 and h2.nCrystalCompton==1) { //true + if (h1.eventID==h2.eventID and h1.eventID==h3.eventID) { //true, phantom-scattered and detector-scattered + if (h1.nPhantomCompton==0 and h2.nPhantomCompton==0 and h3.nPhantomCompton==0) { + if (h1.nCrystalCompton==1 and h2.nCrystalCompton==1 and h3.nCrystalCompton==1) { //true t = kTrue; } else { //detector-scattered @@ -195,51 +223,27 @@ EventType verify_type_of_coincidence(const Hit &h1, const Hit &h2) { } -// based on ComputeKindGATEEvent from -// https://github.com/JPETTomography/castor/blob/castor-3.1.1/src/management/gDataConversionUtilities.cc -EventType verify_type_of_coincidence_castor(const Hit &h1, const Hit &h2) { - - if (h1.eventID != h2.eventID) - // random - return kAccidental; - else - { - if (h1.nPhantomCompton == 0 && h2.nPhantomCompton == 0 && - h1.nPhantomRayleigh == 0 && h2.nPhantomRayleigh == 0 && - h1.nCrystalCompton == 1 && h2.nCrystalCompton == 1 && - h1.nCrystalRayleigh == 0 && h2.nCrystalRayleigh == 0) - //true - return kTrue; - else - { - if (h1.nPhantomCompton == 0 && h2.nPhantomCompton == 0 && - h1.nPhantomRayleigh == 0 && h2.nPhantomRayleigh == 0 && - (h1.nCrystalCompton > 1 || h2.nCrystalCompton > 1 || - h1.nCrystalRayleigh > 0 || h2.nCrystalRayleigh > 0)) - //detector scat - return kDetectorScattered; - if (h1.nPhantomCompton == 1 || h2.nPhantomCompton == 1 || - h1.nPhantomRayleigh == 1 || h2.nPhantomRayleigh == 1) - //single scat - return kPhantomScattered; - if (h1.nPhantomCompton > 1 || h2.nPhantomCompton > 1 || - h1.nPhantomRayleigh > 1 || h2.nPhantomRayleigh > 1) - //multi scat - return kPhantomScattered; - } - } - //unknown - return kUnspecified; -} //================================================================================ // PRINTING //================================================================================ -void print_coincidence(const Hit& h1, const Hit& h2) { +void print_triple_coincidences(const std::vector& hits, PROMPT_E_TH) { + + assert(hits.size() ==3); + + std::vector triple_hits; + bool isGood; + std::tie(isGood, triple_hits) = check_gamma_type(hits, PROMPT_E_TH); + + if(isGood){ cout.setf(ios::fixed); + Hit h1 = triple_hits[0]; + Hit h2 = triple_hits[1]; + Hit h3 = triple_hits[2]; + cout.precision(2); cout << h1.posX/10. << "\t" << h1.posY/10. << "\t" << h1.posZ/10. << "\t"; cout.precision(1); @@ -257,41 +261,23 @@ void print_coincidence(const Hit& h1, const Hit& h2) { cout << h1.edep << "\t"; cout << h2.edep << "\t"; - cout << verify_type_of_coincidence_castor(h1, h2) << "\t"; + cout << verify_type_of_triple_coincidence(h1, h2, h3) << "\t"; cout.precision(2); cout << h1.sourcePosX/10. << "\t" << h1.sourcePosY/10. << "\t" << h1.sourcePosZ/10. << "\t"; - cout << h2.sourcePosX/10. << "\t" << h2.sourcePosY/10. << "\t" << h2.sourcePosZ/10. << "\t"; - - cout << h1.eventID << "\t"; - cout << h2.eventID << "\t"; - cout << h1.nPhantomCompton << "\t"; - cout << h2.nPhantomCompton << "\t"; - cout << h1.nPhantomRayleigh << "\t"; - cout << h2.nPhantomRayleigh << "\t"; - cout << h1.nCrystalCompton << "\t"; - cout << h2.nCrystalCompton << "\t"; - cout << h1.nCrystalRayleigh << "\t"; - cout << h2.nCrystalRayleigh << "\t"; - - cout << h1.rsectorID << "\t"; - cout << h2.rsectorID << "\t"; - cout << h1.moduleID << "\t"; - cout << h2.moduleID << "\t"; - cout << h1.submoduleID << "\t"; - cout << h2.submoduleID << "\t"; - cout << h1.crystalID << "\t"; - cout << h2.crystalID << "\t"; - cout << h1.layerID << "\t"; - cout << h2.layerID << endl; + cout << h2.sourcePosX/10. << "\t" << h2.sourcePosY/10. << "\t" << h2.sourcePosZ/10. << endl; -} - -void print_coincidences(const std::vector& hits) { + cout.precision(2); + cout << h3.posX/10. << "\t" << h3.posY/10. << "\t" << h3.posZ/10. << "\t"; + cout.precision(1); + cout << h3.time << "\t"; + cout << h3.volumeID << "\t"; + cout.precision(2); + cout << h3.edep << "\t"; + cout.precision(2); + cout << h3.sourcePosX/10. << "\t" << h3.sourcePosY/10. << "\t" << h3.sourcePosZ/10. << "\t"; - for (unsigned int i=0; i& hits) { // MAIN ANALYSIS FUNCTION //================================================================================ -void analyze_event(vector &hits, bool hits_are_singles) +void analyze_event(vector &hits, bool hits_are_singles, bool triple_coincidence) { double COMPTON_E_TH = atof(getenv("GOJA_COMPTON_E_TH"))*1e3; + double PROMPT_E_TH = atof(getenv(“GOJA_PROMPT_E_TH”))*1e3; int MAX_N = int(atof(getenv("GOJA_MAX_N"))); int MAX_N0 = int(atof(getenv("GOJA_MAX_N0"))); - int N = 0; - int N0 = 0; - const auto averagingMethod = AveragingMethod(int(atof(getenv("GOJA_AVERAGING_METHOD")))); - const string systemType = string(getenv("GOJA_SYSTEM_TYPE")); - + int N =0; + int N0 =0; sort_hits(hits, "TIME"); std::vector selected_hits; if (hits_are_singles) { - std::tie(N0, N, selected_hits) = select_coincident_singles(hits, COMPTON_E_TH, systemType, averagingMethod); + std::tie(N0, N, selected_hits) = select_coincident_singles(hits, COMPTON_E_TH); } else { std::tie(N0, N, selected_hits) = select_coincident_hits(hits, COMPTON_E_TH); } - // If number of selected singles/hits equals maximum number of events above the fixed energy threshold in the coincidence window - // and number of all hits in the event is equal or smaller than maximum number of events above the noise energy threshold in the coincidence window - // then print the coincidence: + + if (triple_coincidences) { if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); - //if (N>=2 and N<=MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); + } + else { + if (N==MAX_N and N0<=MAX_N0) print_triple_coincidences(selected_hits); + } if (DEBUG) { cout.setf(ios::fixed); @@ -353,51 +339,25 @@ void analyze_event(vector &hits, bool hits_are_singles) void RunTests() { - const double epsilon = 0.0001; - + const double epsilon = 0.000001; Hit h1; h1.edep = 10; - h1.time = 12; - h1.posX = 1; h1.posY = 2; h1.posZ = 3; - Hit h2; h2.edep = 80; - h2.time = 22; - h2.posX = 2; h2.posY = 3; h2.posZ = 4; - Hit h3; - h3.edep = 50; - h3.time = 32; - h3.posX = 7; h3.posY = 4; h3.posZ = 5; - + h2.edep = 50; std::vector hits = {h1, h2, h3}; - for (auto& h : hits) { h.eventID = 10; h.volumeID = 9; } - - auto single1 = merge_hits(hits, kCentroidWinnerNaivelyWeighted); - assert(std::abs(single1.time-22) < epsilon); // (12+22+32)/3 = 22 - assert(std::abs(single1.posX-3.333333333) < epsilon); // (1+2+7)/2 = 3.333333333 - assert(std::abs(single1.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 - - auto single2 = merge_hits(hits, kCentroidWinnerEnergyWeighted); - assert(std::abs(single2.time-24.85714286) < epsilon); // (10*12+80*22+50*32)/140 = 24.85714286 - assert(std::abs(single2.posX-3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 - assert(std::abs(single2.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 - - auto single3 = merge_hits(hits, kCentroidWinnerEnergyWeightedFirstTime); - assert(std::abs(single3.time-12) < epsilon); // min(12,22,32) = 12 - assert(std::abs(single3.posX-3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 - assert(std::abs(single3.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 - - auto single4 = merge_hits(hits, kEnergyWinner); - assert(std::abs(single4.time-22) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 - assert(std::abs(single4.posX-2) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 - assert(single4.eventID == 10); // 10 is a value set manually - assert(single4.volumeID == 9); // 9 is a value set manually - assert(std::abs(single4.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 -} + auto mergedHit = merge_hits(hits, kEnergyWinner); + std::cout << mergedHit.edep << std::endl; + // WK: this test breaks I guess the energy of the merged hit is wrongly assigned + /// It should be the energy of the highest hit, so 80 + // assert(std::abs(mergedHit.edep - 80) < epsilon); + assert(mergedHit.eventID == 10); + assert(mergedHit.volumeID == 9); } +} \ No newline at end of file From 0c7e27a8d0c5ac8f9430ba07f0f57c6263f05a5e Mon Sep 17 00:00:00 2001 From: dominikpanek <3adominikpanek@gmail.com> Date: Wed, 31 Aug 2022 11:39:28 +0200 Subject: [PATCH 3/8] Event analysis code from the develop brnach updated with tripple coincidence code --- goja/EventAnalysis.cpp | 227 ++++++++++++++++++++++++++++++++++------- 1 file changed, 189 insertions(+), 38 deletions(-) diff --git a/goja/EventAnalysis.cpp b/goja/EventAnalysis.cpp index 9f8fb72..0e87bf8 100644 --- a/goja/EventAnalysis.cpp +++ b/goja/EventAnalysis.cpp @@ -53,6 +53,7 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner = kCen Hit h; h.eventID = hits[0].eventID; h.volumeID = hits[0].volumeID; + // Energy of single4 == sum of energies of hits: for (unsigned int i=0; i &hits, const AveragingMethod winner = kCen case kEnergyWinner: { std::vector energies; - // for (unsigned int i = 0; i &hits, const AveragingMethod winner = kCen return h; } - - std::tuple> select_coincident_hits(const vector &hits, double compton_energy_threshold) { std::vector selected_hits; @@ -143,9 +142,8 @@ std::tuple> select_coincident_hits(const vector return std::make_tuple(nb_above_noise_treshold, nb_above_compton_threshold, selected_hits); } -std::tuple> select_coincident_singles(const std::vector &hits, double compton_energy_threshold) { - - const string systemType = string(getenv("GOJA_SYSTEM_TYPE")); +std::tuple> select_coincident_singles (const std::vector &hits, + double compton_energy_threshold, const std::string& systemType, const AveragingMethod averagingMethod = kCentroidWinnerEnergyWeightedFirstTime) { map> singles_tmp; for (unsigned int i=0; i> select_coincident_singles(const std::vec ID = std::to_string(hits[i].volumeID); else if (systemType == "cylindricalPET") ID = std::to_string(hits[i].rsectorID) + '_' + std::to_string(hits[i].layerID); - if (singles_tmp.find(ID) == singles_tmp.end() ) { + if (singles_tmp.find(ID) == singles_tmp.end()) { vector tmp; tmp.push_back(hits[i]); singles_tmp[ID] = tmp; @@ -163,41 +161,77 @@ std::tuple> select_coincident_singles(const std::vec } } - /// WK: Why it is always centroid here? vector singles; map>::iterator it_tmp = singles_tmp.begin(); while(it_tmp != singles_tmp.end()) { - singles.push_back(merge_hits(it_tmp->second, kCentroidWinnerEnergyWeightedFirstTime)); + singles.push_back(merge_hits(it_tmp->second, averagingMethod)); it_tmp++; } return select_coincident_hits(singles, compton_energy_threshold); } -std::tuple> check_gamma_type(const vector &hits, double prompt_energy_threshold) -{ - std::vector triple_hits; +EventType verify_type_of_coincidence(const Hit &h1,const Hit &h2) { - for (unsigned int i=0; i prompt_energy_threshold) { - triple_hits.push_back(hits[i]); - prompt_multiplicity++; - } + else { //accidental + t = kAccidental; } - if(prompt_multiplicity == 1) { - bool isGood = true; - } - else { - bool isGood = false; - } - return std::make_tuple(isGood, triple_hits); - } + + return t; } +// based on ComputeKindGATEEvent from +// https://github.com/JPETTomography/castor/blob/castor-3.1.1/src/management/gDataConversionUtilities.cc +EventType verify_type_of_coincidence_castor(const Hit &h1, const Hit &h2) { + + if (h1.eventID != h2.eventID) + // random + return kAccidental; + else + { + if (h1.nPhantomCompton == 0 && h2.nPhantomCompton == 0 && + h1.nPhantomRayleigh == 0 && h2.nPhantomRayleigh == 0 && + h1.nCrystalCompton == 1 && h2.nCrystalCompton == 1 && + h1.nCrystalRayleigh == 0 && h2.nCrystalRayleigh == 0) + //true + return kTrue; + else + { + if (h1.nPhantomCompton == 0 && h2.nPhantomCompton == 0 && + h1.nPhantomRayleigh == 0 && h2.nPhantomRayleigh == 0 && + (h1.nCrystalCompton > 1 || h2.nCrystalCompton > 1 || + h1.nCrystalRayleigh > 0 || h2.nCrystalRayleigh > 0)) + //detector scat + return kDetectorScattered; + if (h1.nPhantomCompton == 1 || h2.nPhantomCompton == 1 || + h1.nPhantomRayleigh == 1 || h2.nPhantomRayleigh == 1) + //single scat + return kPhantomScattered; + if (h1.nPhantomCompton > 1 || h2.nPhantomCompton > 1 || + h1.nPhantomRayleigh > 1 || h2.nPhantomRayleigh > 1) + //multi scat + return kPhantomScattered; + } + } + //unknown + return kUnspecified; +} + EventType verify_type_of_triple_coincidence(const Hit &h1, const Hit &h2, const Hit &h3) { EventType t = kUnspecified; @@ -223,11 +257,94 @@ EventType verify_type_of_triple_coincidence(const Hit &h1, const Hit &h2, const } +std::tuple> check_gamma_type(const vector &hits, double prompt_energy_threshold) +{ + std::vector triple_hits; + + for (unsigned int i=0; i prompt_energy_threshold) { + triple_hits.push_back(hits[i]); + prompt_multiplicity++; + } + } + if(prompt_multiplicity == 1) { + bool isGood = true; + } + else { + bool isGood = false; + } + return std::make_tuple(isGood, triple_hits); + } + +} //================================================================================ // PRINTING //================================================================================ +void print_coincidence(const Hit& h1, const Hit& h2) { + + cout.setf(ios::fixed); + + cout.precision(2); + cout << h1.posX/10. << "\t" << h1.posY/10. << "\t" << h1.posZ/10. << "\t"; + cout.precision(1); + cout << h1.time << "\t"; + + cout.precision(2); + cout << h2.posX/10. << "\t" << h2.posY/10. << "\t" << h2.posZ/10. << "\t"; + cout.precision(1); + cout << h2.time << "\t"; + + cout << h1.volumeID << "\t"; + cout << h2.volumeID << "\t"; + + cout.precision(2); + cout << h1.edep << "\t"; + cout << h2.edep << "\t"; + + cout << verify_type_of_coincidence_castor(h1, h2) << "\t"; + + cout.precision(2); + cout << h1.sourcePosX/10. << "\t" << h1.sourcePosY/10. << "\t" << h1.sourcePosZ/10. << "\t"; + cout << h2.sourcePosX/10. << "\t" << h2.sourcePosY/10. << "\t" << h2.sourcePosZ/10. << "\t"; + + cout << h1.eventID << "\t"; + cout << h2.eventID << "\t"; + cout << h1.nPhantomCompton << "\t"; + cout << h2.nPhantomCompton << "\t"; + cout << h1.nPhantomRayleigh << "\t"; + cout << h2.nPhantomRayleigh << "\t"; + cout << h1.nCrystalCompton << "\t"; + cout << h2.nCrystalCompton << "\t"; + cout << h1.nCrystalRayleigh << "\t"; + cout << h2.nCrystalRayleigh << "\t"; + + cout << h1.rsectorID << "\t"; + cout << h2.rsectorID << "\t"; + cout << h1.moduleID << "\t"; + cout << h2.moduleID << "\t"; + cout << h1.submoduleID << "\t"; + cout << h2.submoduleID << "\t"; + cout << h1.crystalID << "\t"; + cout << h2.crystalID << "\t"; + cout << h1.layerID << "\t"; + cout << h2.layerID << endl; + +} + +void print_coincidences(const std::vector& hits) { + + for (unsigned int i=0; i& hits, PROMPT_E_TH) { assert(hits.size() ==3); @@ -294,16 +411,24 @@ void analyze_event(vector &hits, bool hits_are_singles, bool triple_coincid int MAX_N0 = int(atof(getenv("GOJA_MAX_N0"))); int N =0; int N0 =0; + const auto averagingMethod = AveragingMethod(int(atof(getenv("GOJA_AVERAGING_METHOD")))); + const string systemType = string(getenv("GOJA_SYSTEM_TYPE")); + sort_hits(hits, "TIME"); std::vector selected_hits; if (hits_are_singles) { - std::tie(N0, N, selected_hits) = select_coincident_singles(hits, COMPTON_E_TH); + std::tie(N0, N, selected_hits) = select_coincident_singles(hits, COMPTON_E_TH, systemType, averagingMethod); } else { std::tie(N0, N, selected_hits) = select_coincident_hits(hits, COMPTON_E_TH); } + // If number of selected singles/hits equals maximum number of events above the fixed energy threshold in the coincidence window + // and number of all hits in the event is equal or smaller than maximum number of events above the noise energy threshold in the coincidence window + // then print the coincidence: + if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); + //if (N>=2 and N<=MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); if (triple_coincidences) { if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); @@ -339,25 +464,51 @@ void analyze_event(vector &hits, bool hits_are_singles, bool triple_coincid void RunTests() { - const double epsilon = 0.000001; + const double epsilon = 0.0001; + Hit h1; h1.edep = 10; + h1.time = 12; + h1.posX = 1; h1.posY = 2; h1.posZ = 3; + Hit h2; h2.edep = 80; + h2.time = 22; + h2.posX = 2; h2.posY = 3; h2.posZ = 4; + Hit h3; - h2.edep = 50; + h3.edep = 50; + h3.time = 32; + h3.posX = 7; h3.posY = 4; h3.posZ = 5; + std::vector hits = {h1, h2, h3}; + for (auto& h : hits) { h.eventID = 10; h.volumeID = 9; } - auto mergedHit = merge_hits(hits, kEnergyWinner); - std::cout << mergedHit.edep << std::endl; - // WK: this test breaks I guess the energy of the merged hit is wrongly assigned - /// It should be the energy of the highest hit, so 80 - // assert(std::abs(mergedHit.edep - 80) < epsilon); - assert(mergedHit.eventID == 10); - assert(mergedHit.volumeID == 9); + + auto single1 = merge_hits(hits, kCentroidWinnerNaivelyWeighted); + assert(std::abs(single1.time-22) < epsilon); // (12+22+32)/3 = 22 + assert(std::abs(single1.posX-3.333333333) < epsilon); // (1+2+7)/2 = 3.333333333 + assert(std::abs(single1.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 + + auto single2 = merge_hits(hits, kCentroidWinnerEnergyWeighted); + assert(std::abs(single2.time-24.85714286) < epsilon); // (10*12+80*22+50*32)/140 = 24.85714286 + assert(std::abs(single2.posX-3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 + assert(std::abs(single2.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 + + auto single3 = merge_hits(hits, kCentroidWinnerEnergyWeightedFirstTime); + assert(std::abs(single3.time-12) < epsilon); // min(12,22,32) = 12 + assert(std::abs(single3.posX-3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 + assert(std::abs(single3.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 + + auto single4 = merge_hits(hits, kEnergyWinner); + assert(std::abs(single4.time-22) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 + assert(std::abs(single4.posX-2) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 + assert(single4.eventID == 10); // 10 is a value set manually + assert(single4.volumeID == 9); // 9 is a value set manually + assert(std::abs(single4.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 } } \ No newline at end of file From d2d0f2a0a1bef99398061d16a2fd92e30435b661 Mon Sep 17 00:00:00 2001 From: dominikpanek <3adominikpanek@gmail.com> Date: Wed, 31 Aug 2022 11:42:51 +0200 Subject: [PATCH 4/8] Event analysis code from develop branch, updated with tripple coincidences and remarks from pawel --- goja/EventAnalysis.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/goja/EventAnalysis.cpp b/goja/EventAnalysis.cpp index 0e87bf8..2f48bcb 100644 --- a/goja/EventAnalysis.cpp +++ b/goja/EventAnalysis.cpp @@ -106,7 +106,7 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner = kCen case kEnergyWinner: { std::vector energies; - for (unsigned int i = 0; i> select_coincident_singles (const std::vec return select_coincident_hits(singles, compton_energy_threshold); } -EventType verify_type_of_coincidence(const Hit &h1,const Hit &h2) { +EventType verify_type_of_coincidence(const Hit &h1, const Hit &h2) { EventType t = kUnspecified; @@ -409,8 +409,8 @@ void analyze_event(vector &hits, bool hits_are_singles, bool triple_coincid double PROMPT_E_TH = atof(getenv(“GOJA_PROMPT_E_TH”))*1e3; int MAX_N = int(atof(getenv("GOJA_MAX_N"))); int MAX_N0 = int(atof(getenv("GOJA_MAX_N0"))); - int N =0; - int N0 =0; + int N = 0; + int N0 = 0; const auto averagingMethod = AveragingMethod(int(atof(getenv("GOJA_AVERAGING_METHOD")))); const string systemType = string(getenv("GOJA_SYSTEM_TYPE")); @@ -510,5 +510,4 @@ void RunTests() assert(single4.eventID == 10); // 10 is a value set manually assert(single4.volumeID == 9); // 9 is a value set manually assert(std::abs(single4.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 -} } \ No newline at end of file From 21b2fcef34f57ed14888a999a667494d0e5b1707 Mon Sep 17 00:00:00 2001 From: dominikpanek <3adominikpanek@gmail.com> Date: Tue, 11 Oct 2022 12:56:41 +0200 Subject: [PATCH 5/8] Tripple coincidence code updated with the remarks from mateusz. The changes were applied in the EventAnalysis.cpp and the declarations of new functions were added in th EventAnalysis.h --- goja/EventAnalysis.cpp | 97 +++++++++++++++++++++--------------------- goja/EventAnalysis.h | 16 ++++++- 2 files changed, 64 insertions(+), 49 deletions(-) diff --git a/goja/EventAnalysis.cpp b/goja/EventAnalysis.cpp index 2f48bcb..4bec795 100644 --- a/goja/EventAnalysis.cpp +++ b/goja/EventAnalysis.cpp @@ -170,29 +170,36 @@ std::tuple> select_coincident_singles (const std::vec return select_coincident_hits(singles, compton_energy_threshold); } -EventType verify_type_of_coincidence(const Hit &h1, const Hit &h2) { - - EventType t = kUnspecified; +bool all_same_event(const vector& hits) { + unsigned int nhits = hits.size(); + if ( nhits == 0 ) { + return false; + } + const int ref_id = hits.at(0).eventID; + return std::all_of(hits.cbegin(), hits.cend(), [ref_id](Hit& hit){ return hit.eventID == ref_id; }); +} - if (h1.eventID==h2.eventID) { //true, phantom-scattered and detector-scattered - if (h1.nPhantomCompton==0 and h2.nPhantomCompton==0) { - if (h1.nCrystalCompton==1 and h2.nCrystalCompton==1) { //true - t = kTrue; - } - else { //detector-scattered - t = kDetectorScattered; - } - } - else { //phantom-scattered - t = kPhantomScattered; - } - } - else { //accidental - t = kAccidental; - } +bool not_phantom_compton(const vector& hits) { + return std::all_of(hits.cbegin(), hits.cend(), [](Hit& hit){ return hit.nPhantomCompton == 0; }); +} - return t; +bool all_crystal_compton(const vector& hits) { + return std::all_of(hits.cbegin(), hits.cend(), [](Hit& hit){ return hit.nCrystalCompton == 1; }); +} +EventType verify_type_of_coincidence(const std::vector& hits) { + if ( all_same_event(hits) ) { + //true, phantom-scattered and detector-scattered + if ( not_phantom_compton(hits) ) { + if ( all_crystal_compton(hits) ) { + return kTrue; + } + return kDetectorScattered; + } + return kPhantomScattered; + } + //otherwise --> accidental + return kAccidental; } // based on ComputeKindGATEEvent from @@ -231,32 +238,24 @@ EventType verify_type_of_coincidence_castor(const Hit &h1, const Hit &h2) { //unknown return kUnspecified; } - -EventType verify_type_of_triple_coincidence(const Hit &h1, const Hit &h2, const Hit &h3) { - - EventType t = kUnspecified; - - if (h1.eventID==h2.eventID and h1.eventID==h3.eventID) { //true, phantom-scattered and detector-scattered - if (h1.nPhantomCompton==0 and h2.nPhantomCompton==0 and h3.nPhantomCompton==0) { - if (h1.nCrystalCompton==1 and h2.nCrystalCompton==1 and h3.nCrystalCompton==1) { //true - t = kTrue; - } - else { //detector-scattered - t = kDetectorScattered; - } - } - else { //phantom-scattered - t = kPhantomScattered; - } - } - else { //accidental - t = kAccidental; +// Code checking the type of 3 gamma events. Gammas are assigned to four classes. True, Phantom Scattered, Detector +// Scattered and Accidental (Random). +EventType verify_type_of_triple_coincidence(const std::vector& hits) { + if ( all_same_event(hits) ) { + //true, phantom-scattered and detector-scattered + if ( not_phantom_compton(hits) ) { + if ( all_crystal_compton(hits) ) { + return kTrue; + } + return kDetectorScattered; } - - return t; - + return kPhantomScattered; + } + //otherwise --> accidental + return kAccidental; } + std::tuple> check_gamma_type(const vector &hits, double prompt_energy_threshold) { std::vector triple_hits; @@ -345,9 +344,11 @@ void print_coincidences(const std::vector& hits) { } +// Code printing specific information about the hits. The information include X, Y and Z position; detection time, +// energy deposition etc. void print_triple_coincidences(const std::vector& hits, PROMPT_E_TH) { - assert(hits.size() ==3); + assert(hits.size() == 3); std::vector triple_hits; bool isGood; @@ -402,7 +403,7 @@ void print_triple_coincidences(const std::vector& hits, PROMPT_E_TH) { // MAIN ANALYSIS FUNCTION //================================================================================ -void analyze_event(vector &hits, bool hits_are_singles, bool triple_coincidence) +void analyze_event(const vector &hits, const bool hits_are_singles, const bool triple_coincidence) { double COMPTON_E_TH = atof(getenv("GOJA_COMPTON_E_TH"))*1e3; @@ -431,10 +432,10 @@ void analyze_event(vector &hits, bool hits_are_singles, bool triple_coincid //if (N>=2 and N<=MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); if (triple_coincidences) { - if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); + if (N==MAX_N and N0<=MAX_N0) {print_coincidences(selected_hits);} } else { - if (N==MAX_N and N0<=MAX_N0) print_triple_coincidences(selected_hits); + if (N==MAX_N and N0<=MAX_N0) {print_triple_coincidences(selected_hits)}; } if (DEBUG) { @@ -510,4 +511,4 @@ void RunTests() assert(single4.eventID == 10); // 10 is a value set manually assert(single4.volumeID == 9); // 9 is a value set manually assert(std::abs(single4.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 -} \ No newline at end of file +} diff --git a/goja/EventAnalysis.h b/goja/EventAnalysis.h index 608b175..3df4c3d 100755 --- a/goja/EventAnalysis.h +++ b/goja/EventAnalysis.h @@ -41,8 +41,20 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner); /// returns number of singles above noise energy threshold, number of singles above Compton energy threshold, and selected singles std::tuple> select_coincident_singles(const std::vector &hits, double compton_energy_threshold); + /// returns the event ID of the hit + bool all_same_event(const std::vector& hits); + + /// returns 0 if the hit was not scattered in phantom + bool not_phantom_compton(const std::vector& hits); + + /// returns 1 if there was not any scattering in the detector + bool all_crystal_compton(const std::vector& hits); + // Below function contains simple definition of coincidence that does not take into account Rayleigh scatterings - EventType verify_type_of_coincidence(const Hit &h1, const Hit &h2); + EventType verify_type_of_coincidence(const std::vector& hits); + + //Below is the function that assigns the classes of coincidences when the 3-gamma event occurs. + EventType verify_type_of_triple_coincidence(const std::vector& hits); // Below function contains definition of coincidence that takes into account Rayleigh scatterings, // the definition is based on the snippet from the CASToR software: function ComputeKindGATEEvent from @@ -51,7 +63,9 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner); void print_coincidence(const Hit& h1, const Hit& h2); void print_coincidences(const std::vector& hits); + void print_triple_coincidences(const std::vector& hits, PROMPT_E_TH); void analyze_event(std::vector &hits, bool hits_are_singles = true); + void analyze_event(const std::vector &hits, const bool hits_are_singles, const bool triple_coincidence); }; From a7eec1ac9373fc1c4ea26aa7b5d27f40d9c06be4 Mon Sep 17 00:00:00 2001 From: dominikpanek <3adominikpanek@gmail.com> Date: Fri, 14 Oct 2022 10:37:40 +0200 Subject: [PATCH 6/8] Functions EventAnalysis.cpp and EventAnalysis.h corrected accoring to Mateusz's remarks. --- goja/EventAnalysis.cpp | 2 +- goja/EventAnalysis.h | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/goja/EventAnalysis.cpp b/goja/EventAnalysis.cpp index 4bec795..edf4133 100644 --- a/goja/EventAnalysis.cpp +++ b/goja/EventAnalysis.cpp @@ -346,7 +346,7 @@ void print_coincidences(const std::vector& hits) { // Code printing specific information about the hits. The information include X, Y and Z position; detection time, // energy deposition etc. -void print_triple_coincidences(const std::vector& hits, PROMPT_E_TH) { +void print_triple_coincidences(const std::vector& hits, double PROMPT_E_TH) { assert(hits.size() == 3); diff --git a/goja/EventAnalysis.h b/goja/EventAnalysis.h index 3df4c3d..1ccb577 100755 --- a/goja/EventAnalysis.h +++ b/goja/EventAnalysis.h @@ -41,19 +41,19 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner); /// returns number of singles above noise energy threshold, number of singles above Compton energy threshold, and selected singles std::tuple> select_coincident_singles(const std::vector &hits, double compton_energy_threshold); - /// returns the event ID of the hit + /// function checks if hits come from the same annihilation event bool all_same_event(const std::vector& hits); - /// returns 0 if the hit was not scattered in phantom + /// returns 1 if the hit was not scattered in phantom bool not_phantom_compton(const std::vector& hits); /// returns 1 if there was not any scattering in the detector bool all_crystal_compton(const std::vector& hits); - // Below function contains simple definition of coincidence that does not take into account Rayleigh scatterings + // Function below contains simple definition of coincidence that does not take into account Rayleigh scatterings EventType verify_type_of_coincidence(const std::vector& hits); - //Below is the function that assigns the classes of coincidences when the 3-gamma event occurs. + // Function below assigns the classes of coincidences when the 3-gamma event occurs. EventType verify_type_of_triple_coincidence(const std::vector& hits); // Below function contains definition of coincidence that takes into account Rayleigh scatterings, @@ -63,7 +63,7 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner); void print_coincidence(const Hit& h1, const Hit& h2); void print_coincidences(const std::vector& hits); - void print_triple_coincidences(const std::vector& hits, PROMPT_E_TH); + void print_triple_coincidences(const std::vector& hits, double PROMPT_E_TH); void analyze_event(std::vector &hits, bool hits_are_singles = true); void analyze_event(const std::vector &hits, const bool hits_are_singles, const bool triple_coincidence); From 549e819f6f1d992ed2441b60b8ad991dffebd114 Mon Sep 17 00:00:00 2001 From: dominikpanek <3adominikpanek@gmail.com> Date: Sun, 6 Nov 2022 13:00:46 +0100 Subject: [PATCH 7/8] Updated version of the tripple coincidence code. Error were fixed. Code compiles now. --- goja/EventAnalysis.cpp | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/goja/EventAnalysis.cpp b/goja/EventAnalysis.cpp index edf4133..085f52a 100644 --- a/goja/EventAnalysis.cpp +++ b/goja/EventAnalysis.cpp @@ -176,15 +176,15 @@ bool all_same_event(const vector& hits) { return false; } const int ref_id = hits.at(0).eventID; - return std::all_of(hits.cbegin(), hits.cend(), [ref_id](Hit& hit){ return hit.eventID == ref_id; }); + return std::all_of(hits.cbegin(), hits.cend(), [ref_id](const Hit& hit){ return hit.eventID == ref_id; }); } bool not_phantom_compton(const vector& hits) { - return std::all_of(hits.cbegin(), hits.cend(), [](Hit& hit){ return hit.nPhantomCompton == 0; }); + return std::all_of(hits.cbegin(), hits.cend(), [](const Hit& hit){ return hit.nPhantomCompton == 0; }); } bool all_crystal_compton(const vector& hits) { - return std::all_of(hits.cbegin(), hits.cend(), [](Hit& hit){ return hit.nCrystalCompton == 1; }); + return std::all_of(hits.cbegin(), hits.cend(), [](const Hit& hit){ return hit.nCrystalCompton == 1; }); } EventType verify_type_of_coincidence(const std::vector& hits) { @@ -256,7 +256,7 @@ EventType verify_type_of_triple_coincidence(const std::vector& hits) { } -std::tuple> check_gamma_type(const vector &hits, double prompt_energy_threshold) +std::tuple> check_gamma_type(const vector &hits, double prompt_energy_threshold) { std::vector triple_hits; @@ -270,6 +270,7 @@ std::tuple> check_gamma_type(const vector &hits, dou prompt_multiplicity++; } } + bool isGood; if(prompt_multiplicity == 1) { bool isGood = true; } @@ -306,7 +307,7 @@ void print_coincidence(const Hit& h1, const Hit& h2) { cout << h1.edep << "\t"; cout << h2.edep << "\t"; - cout << verify_type_of_coincidence_castor(h1, h2) << "\t"; + cout << event_analysis::verify_type_of_coincidence_castor(h1, h2) << "\t"; cout.precision(2); cout << h1.sourcePosX/10. << "\t" << h1.sourcePosY/10. << "\t" << h1.sourcePosZ/10. << "\t"; @@ -352,7 +353,7 @@ void print_triple_coincidences(const std::vector& hits, double PROMPT_E_TH) std::vector triple_hits; bool isGood; - std::tie(isGood, triple_hits) = check_gamma_type(hits, PROMPT_E_TH); + std::tie(isGood, triple_hits) = event_analysis::check_gamma_type(hits, PROMPT_E_TH); if(isGood){ @@ -379,7 +380,7 @@ void print_triple_coincidences(const std::vector& hits, double PROMPT_E_TH) cout << h1.edep << "\t"; cout << h2.edep << "\t"; - cout << verify_type_of_triple_coincidence(h1, h2, h3) << "\t"; + cout << event_analysis::verify_type_of_triple_coincidence(triple_hits) << "\t"; cout.precision(2); cout << h1.sourcePosX/10. << "\t" << h1.sourcePosY/10. << "\t" << h1.sourcePosZ/10. << "\t"; @@ -407,22 +408,22 @@ void analyze_event(const vector &hits, const bool hits_are_singles, const b { double COMPTON_E_TH = atof(getenv("GOJA_COMPTON_E_TH"))*1e3; - double PROMPT_E_TH = atof(getenv(“GOJA_PROMPT_E_TH”))*1e3; + double PROMPT_E_TH = atof(getenv("GOJA_PROMPT_E_TH"))*1e3; int MAX_N = int(atof(getenv("GOJA_MAX_N"))); int MAX_N0 = int(atof(getenv("GOJA_MAX_N0"))); int N = 0; int N0 = 0; - const auto averagingMethod = AveragingMethod(int(atof(getenv("GOJA_AVERAGING_METHOD")))); + const auto averagingMethod = event_analysis::AveragingMethod(int(atof(getenv("GOJA_AVERAGING_METHOD")))); const string systemType = string(getenv("GOJA_SYSTEM_TYPE")); - sort_hits(hits, "TIME"); + event_analysis::sort_hits(const_cast &>(hits), "TIME"); std::vector selected_hits; if (hits_are_singles) { std::tie(N0, N, selected_hits) = select_coincident_singles(hits, COMPTON_E_TH, systemType, averagingMethod); } else { - std::tie(N0, N, selected_hits) = select_coincident_hits(hits, COMPTON_E_TH); + std::tie(N0, N, selected_hits) = event_analysis::select_coincident_hits(hits, COMPTON_E_TH); } // If number of selected singles/hits equals maximum number of events above the fixed energy threshold in the coincidence window @@ -431,11 +432,11 @@ void analyze_event(const vector &hits, const bool hits_are_singles, const b if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); //if (N>=2 and N<=MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); - if (triple_coincidences) { + if (triple_coincidence) { if (N==MAX_N and N0<=MAX_N0) {print_coincidences(selected_hits);} } else { - if (N==MAX_N and N0<=MAX_N0) {print_triple_coincidences(selected_hits)}; + if (N==MAX_N and N0<=MAX_N0) {print_triple_coincidences(selected_hits, PROMPT_E_TH);}; } if (DEBUG) { @@ -490,22 +491,22 @@ void RunTests() h.volumeID = 9; } - auto single1 = merge_hits(hits, kCentroidWinnerNaivelyWeighted); + auto single1 = merge_hits(hits, event_analysis::kCentroidWinnerNaivelyWeighted); assert(std::abs(single1.time-22) < epsilon); // (12+22+32)/3 = 22 assert(std::abs(single1.posX-3.333333333) < epsilon); // (1+2+7)/2 = 3.333333333 assert(std::abs(single1.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 - auto single2 = merge_hits(hits, kCentroidWinnerEnergyWeighted); + auto single2 = merge_hits(hits, event_analysis::kCentroidWinnerEnergyWeighted); assert(std::abs(single2.time-24.85714286) < epsilon); // (10*12+80*22+50*32)/140 = 24.85714286 assert(std::abs(single2.posX-3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 assert(std::abs(single2.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 - auto single3 = merge_hits(hits, kCentroidWinnerEnergyWeightedFirstTime); + auto single3 = merge_hits(hits, event_analysis::kCentroidWinnerEnergyWeightedFirstTime); assert(std::abs(single3.time-12) < epsilon); // min(12,22,32) = 12 assert(std::abs(single3.posX-3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 assert(std::abs(single3.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 - auto single4 = merge_hits(hits, kEnergyWinner); + auto single4 = merge_hits(hits, event_analysis::kEnergyWinner); assert(std::abs(single4.time-22) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 assert(std::abs(single4.posX-2) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 assert(single4.eventID == 10); // 10 is a value set manually From cf2b60fb9c6040c6df61c048d1136eb8bc69b203 Mon Sep 17 00:00:00 2001 From: dominikpanek <3adominikpanek@gmail.com> Date: Wed, 14 Dec 2022 10:32:11 +0100 Subject: [PATCH 8/8] Updated version of tripple coincidence code --- goja/EventAnalysis.cpp | 453 ++++++++++++++++++++++------------------- goja/EventAnalysis.h | 3 +- 2 files changed, 248 insertions(+), 208 deletions(-) diff --git a/goja/EventAnalysis.cpp b/goja/EventAnalysis.cpp index 085f52a..dae6cf4 100644 --- a/goja/EventAnalysis.cpp +++ b/goja/EventAnalysis.cpp @@ -14,15 +14,18 @@ using namespace std; // struct used for sorting vectors of hits (functor) -namespace event_analysis{ +namespace event_analysis +{ -struct EntityComp { +struct EntityComp +{ string property; - EntityComp(string property) {this->property = property;} - bool operator()(const Hit& h1, const Hit& h2) const { + EntityComp(string property) { this->property = property; } + bool operator()(const Hit& h1, const Hit& h2) const + { bool val = false; - if(property == "TIME") + if (property == "TIME") val = h1.time < h2.time; else if (property == "EDEP") val = h1.edep < h2.edep; @@ -30,36 +33,37 @@ struct EntityComp { val = h1.eventID < h2.eventID; return val; } - }; -void sort_hits(vector &hits, string key) { - - sort(hits.begin(), hits.end(), EntityComp(key)); - -} +void sort_hits(vector& hits, string key) { sort(hits.begin(), hits.end(), EntityComp(key)); } -Hit merge_hits(const std::vector &hits, const AveragingMethod winner = kCentroidWinnerEnergyWeightedFirstTime) { +Hit merge_hits(const std::vector& hits, const AveragingMethod winner = kCentroidWinnerEnergyWeightedFirstTime) +{ const unsigned int nHits = hits.size(); /// WK: temporary checks for debugs to be switched off - assert(nHits > 0); - if (nHits > 0) { + assert(nHits > 0); + if (nHits > 0) + { auto eID = hits[0].eventID; auto vID = hits[0].volumeID; - assert(std::all_of(hits.begin(), hits.end(), [eID](const Hit& h)->bool { return h.eventID == eID; })); - assert(std::all_of(hits.begin(), hits.end(), [vID](const Hit& h)->bool { return h.volumeID == vID; })); - } + assert(std::all_of(hits.begin(), hits.end(), [eID](const Hit& h) -> bool { return h.eventID == eID; })); + assert(std::all_of(hits.begin(), hits.end(), [vID](const Hit& h) -> bool { return h.volumeID == vID; })); + } /// WK: end Hit h; h.eventID = hits[0].eventID; h.volumeID = hits[0].volumeID; // Energy of single4 == sum of energies of hits: - for (unsigned int i=0; i &hits, const AveragingMethod winner = kCen break; } - case kCentroidWinnerEnergyWeighted: { - for (unsigned int i = 0; i < nHits; i++) { + case kCentroidWinnerEnergyWeighted: + { + for (unsigned int i = 0; i < nHits; i++) + { h.time += hits[i].time * hits[i].edep; h.posX += hits[i].posX * hits[i].edep; h.posY += hits[i].posY * hits[i].edep; @@ -86,14 +92,15 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner = kCen break; } - case kCentroidWinnerEnergyWeightedFirstTime: { + case kCentroidWinnerEnergyWeightedFirstTime: + { std::vector times; for (unsigned int i = 0; i < nHits; i++) times.push_back(hits[i].time); - unsigned int min_index = std::distance( - times.begin(), std::min_element(times.begin(), times.end())); + unsigned int min_index = std::distance(times.begin(), std::min_element(times.begin(), times.end())); h.time = hits[min_index].time; - for (unsigned int i = 0; i < nHits; i++) { + for (unsigned int i = 0; i < nHits; i++) + { h.posX += hits[i].posX * hits[i].edep; h.posY += hits[i].posY * hits[i].edep; h.posZ += hits[i].posZ * hits[i].edep; @@ -104,11 +111,12 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner = kCen break; } - case kEnergyWinner: { + case kEnergyWinner: + { std::vector energies; - for (unsigned int i = 0; i &hits, const AveragingMethod winner = kCen break; } - default: { + default: + { /// should never happen assert(1 == 0); break; @@ -130,173 +139,186 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner = kCen return h; } -std::tuple> select_coincident_hits(const vector &hits, double compton_energy_threshold) +std::tuple> select_coincident_hits(const vector& hits, double compton_energy_threshold) { std::vector selected_hits; int nb_above_noise_treshold = hits.size(); - for (unsigned int i=0; i compton_energy_threshold) selected_hits.push_back(hits[i]); + for (unsigned int i = 0; i < hits.size(); i++) + { + if (hits[i].edep > compton_energy_threshold) + selected_hits.push_back(hits[i]); } int nb_above_compton_threshold = selected_hits.size(); return std::make_tuple(nb_above_noise_treshold, nb_above_compton_threshold, selected_hits); } -std::tuple> select_coincident_singles (const std::vector &hits, - double compton_energy_threshold, const std::string& systemType, const AveragingMethod averagingMethod = kCentroidWinnerEnergyWeightedFirstTime) { +std::tuple> select_coincident_singles(const std::vector& hits, double compton_energy_threshold, + const std::string& systemType, + const AveragingMethod averagingMethod = kCentroidWinnerEnergyWeightedFirstTime) +{ map> singles_tmp; - for (unsigned int i=0; i tmp; tmp.push_back(hits[i]); singles_tmp[ID] = tmp; - } else { + } + else + { singles_tmp[ID].push_back(hits[i]); } } vector singles; map>::iterator it_tmp = singles_tmp.begin(); - while(it_tmp != singles_tmp.end()) { + while (it_tmp != singles_tmp.end()) + { singles.push_back(merge_hits(it_tmp->second, averagingMethod)); it_tmp++; } return select_coincident_hits(singles, compton_energy_threshold); } -bool all_same_event(const vector& hits) { - unsigned int nhits = hits.size(); - if ( nhits == 0 ) { - return false; - } - const int ref_id = hits.at(0).eventID; - return std::all_of(hits.cbegin(), hits.cend(), [ref_id](const Hit& hit){ return hit.eventID == ref_id; }); +bool all_same_event(const vector& hits) +{ + unsigned int nhits = hits.size(); + if (nhits == 0) + { + return false; + } + const int ref_id = hits.at(0).eventID; + return std::all_of(hits.cbegin(), hits.cend(), [ref_id](const Hit& hit) { return hit.eventID == ref_id; }); } -bool not_phantom_compton(const vector& hits) { - return std::all_of(hits.cbegin(), hits.cend(), [](const Hit& hit){ return hit.nPhantomCompton == 0; }); +bool not_phantom_compton(const vector& hits) +{ + return std::all_of(hits.cbegin(), hits.cend(), [](const Hit& hit) { return hit.nPhantomCompton == 0; }); } -bool all_crystal_compton(const vector& hits) { - return std::all_of(hits.cbegin(), hits.cend(), [](const Hit& hit){ return hit.nCrystalCompton == 1; }); +bool all_crystal_compton(const vector& hits) +{ + return std::all_of(hits.cbegin(), hits.cend(), [](const Hit& hit) { return hit.nCrystalCompton == 1; }); } -EventType verify_type_of_coincidence(const std::vector& hits) { - if ( all_same_event(hits) ) { - //true, phantom-scattered and detector-scattered - if ( not_phantom_compton(hits) ) { - if ( all_crystal_compton(hits) ) { - return kTrue; - } - return kDetectorScattered; +EventType verify_type_of_coincidence(const std::vector& hits) +{ + if (all_same_event(hits)) + { + // true, phantom-scattered and detector-scattered + if (not_phantom_compton(hits)) + { + if (all_crystal_compton(hits)) + { + return kTrue; + } + return kDetectorScattered; + } + return kPhantomScattered; } - return kPhantomScattered; - } - //otherwise --> accidental - return kAccidental; + // otherwise --> accidental + return kAccidental; } // based on ComputeKindGATEEvent from // https://github.com/JPETTomography/castor/blob/castor-3.1.1/src/management/gDataConversionUtilities.cc -EventType verify_type_of_coincidence_castor(const Hit &h1, const Hit &h2) { +EventType verify_type_of_coincidence_castor(const Hit& h1, const Hit& h2) +{ if (h1.eventID != h2.eventID) // random return kAccidental; else { - if (h1.nPhantomCompton == 0 && h2.nPhantomCompton == 0 && - h1.nPhantomRayleigh == 0 && h2.nPhantomRayleigh == 0 && - h1.nCrystalCompton == 1 && h2.nCrystalCompton == 1 && - h1.nCrystalRayleigh == 0 && h2.nCrystalRayleigh == 0) - //true - return kTrue; + if (h1.nPhantomCompton == 0 && h2.nPhantomCompton == 0 && h1.nPhantomRayleigh == 0 && h2.nPhantomRayleigh == 0 && h1.nCrystalCompton == 1 && + h2.nCrystalCompton == 1 && h1.nCrystalRayleigh == 0 && h2.nCrystalRayleigh == 0) + // true + return kTrue; else { - if (h1.nPhantomCompton == 0 && h2.nPhantomCompton == 0 && - h1.nPhantomRayleigh == 0 && h2.nPhantomRayleigh == 0 && - (h1.nCrystalCompton > 1 || h2.nCrystalCompton > 1 || - h1.nCrystalRayleigh > 0 || h2.nCrystalRayleigh > 0)) - //detector scat - return kDetectorScattered; - if (h1.nPhantomCompton == 1 || h2.nPhantomCompton == 1 || - h1.nPhantomRayleigh == 1 || h2.nPhantomRayleigh == 1) - //single scat - return kPhantomScattered; - if (h1.nPhantomCompton > 1 || h2.nPhantomCompton > 1 || - h1.nPhantomRayleigh > 1 || h2.nPhantomRayleigh > 1) - //multi scat - return kPhantomScattered; + if (h1.nPhantomCompton == 0 && h2.nPhantomCompton == 0 && h1.nPhantomRayleigh == 0 && h2.nPhantomRayleigh == 0 && + (h1.nCrystalCompton > 1 || h2.nCrystalCompton > 1 || h1.nCrystalRayleigh > 0 || h2.nCrystalRayleigh > 0)) + // detector scat + return kDetectorScattered; + if (h1.nPhantomCompton == 1 || h2.nPhantomCompton == 1 || h1.nPhantomRayleigh == 1 || h2.nPhantomRayleigh == 1) + // single scat + return kPhantomScattered; + if (h1.nPhantomCompton > 1 || h2.nPhantomCompton > 1 || h1.nPhantomRayleigh > 1 || h2.nPhantomRayleigh > 1) + // multi scat + return kPhantomScattered; } } - //unknown + // unknown return kUnspecified; } // Code checking the type of 3 gamma events. Gammas are assigned to four classes. True, Phantom Scattered, Detector // Scattered and Accidental (Random). -EventType verify_type_of_triple_coincidence(const std::vector& hits) { - if ( all_same_event(hits) ) { - //true, phantom-scattered and detector-scattered - if ( not_phantom_compton(hits) ) { - if ( all_crystal_compton(hits) ) { - return kTrue; - } - return kDetectorScattered; +EventType verify_type_of_triple_coincidence(const std::vector& hits) +{ + if (all_same_event(hits)) + { + // true, phantom-scattered and detector-scattered + if (not_phantom_compton(hits)) + { + if (all_crystal_compton(hits)) + { + return kTrue; + } + return kDetectorScattered; + } + return kPhantomScattered; } - return kPhantomScattered; - } - //otherwise --> accidental - return kAccidental; + // otherwise --> accidental + return kAccidental; } - -std::tuple> check_gamma_type(const vector &hits, double prompt_energy_threshold) +std::tuple> check_gamma_type(const vector& hits, double prompt_energy_threshold) { std::vector triple_hits; - for (unsigned int i=0; i prompt_energy_threshold) { - triple_hits.push_back(hits[i]); - prompt_multiplicity++; - } + int prompt_multiplicity = 0; + for (unsigned int i = 0; i < hits.size(); i++) + { + if (hits[i].edep > prompt_energy_threshold) + { + triple_hits.push_back(hits[i]); + prompt_multiplicity++; + } } - bool isGood; - if(prompt_multiplicity == 1) { - bool isGood = true; - } - else { - bool isGood = false; - } - return std::make_tuple(isGood, triple_hits); - } + bool isGood = prompt_multiplicity == 1; + return std::make_tuple(isGood, triple_hits); } //================================================================================ // PRINTING //================================================================================ -void print_coincidence(const Hit& h1, const Hit& h2) { +void print_coincidence(const Hit& h1, const Hit& h2) +{ cout.setf(ios::fixed); cout.precision(2); - cout << h1.posX/10. << "\t" << h1.posY/10. << "\t" << h1.posZ/10. << "\t"; + cout << h1.posX / 10. << "\t" << h1.posY / 10. << "\t" << h1.posZ / 10. << "\t"; cout.precision(1); cout << h1.time << "\t"; cout.precision(2); - cout << h2.posX/10. << "\t" << h2.posY/10. << "\t" << h2.posZ/10. << "\t"; + cout << h2.posX / 10. << "\t" << h2.posY / 10. << "\t" << h2.posZ / 10. << "\t"; cout.precision(1); cout << h2.time << "\t"; @@ -310,8 +332,8 @@ void print_coincidence(const Hit& h1, const Hit& h2) { cout << event_analysis::verify_type_of_coincidence_castor(h1, h2) << "\t"; cout.precision(2); - cout << h1.sourcePosX/10. << "\t" << h1.sourcePosY/10. << "\t" << h1.sourcePosZ/10. << "\t"; - cout << h2.sourcePosX/10. << "\t" << h2.sourcePosY/10. << "\t" << h2.sourcePosZ/10. << "\t"; + cout << h1.sourcePosX / 10. << "\t" << h1.sourcePosY / 10. << "\t" << h1.sourcePosZ / 10. << "\t"; + cout << h2.sourcePosX / 10. << "\t" << h2.sourcePosY / 10. << "\t" << h2.sourcePosZ / 10. << "\t"; cout << h1.eventID << "\t"; cout << h2.eventID << "\t"; @@ -334,134 +356,146 @@ void print_coincidence(const Hit& h1, const Hit& h2) { cout << h2.crystalID << "\t"; cout << h1.layerID << "\t"; cout << h2.layerID << endl; - } -void print_coincidences(const std::vector& hits) { +void print_coincidences(const std::vector& hits) +{ - for (unsigned int i=0; i& hits, double PROMPT_E_TH) { - - assert(hits.size() == 3); - - std::vector triple_hits; - bool isGood; - std::tie(isGood, triple_hits) = event_analysis::check_gamma_type(hits, PROMPT_E_TH); - - if(isGood){ +void print_triple_coincidences(const std::vector& hits, double PROMPT_E_TH) +{ + assert(hits.size() == 3); - cout.setf(ios::fixed); + std::vector triple_hits; + bool isGood; + std::tie(isGood, triple_hits) = event_analysis::check_gamma_type(hits, PROMPT_E_TH); - Hit h1 = triple_hits[0]; - Hit h2 = triple_hits[1]; - Hit h3 = triple_hits[2]; + if (isGood) + { - cout.precision(2); - cout << h1.posX/10. << "\t" << h1.posY/10. << "\t" << h1.posZ/10. << "\t"; - cout.precision(1); - cout << h1.time << "\t"; + cout.setf(ios::fixed); - cout.precision(2); - cout << h2.posX/10. << "\t" << h2.posY/10. << "\t" << h2.posZ/10. << "\t"; - cout.precision(1); - cout << h2.time << "\t"; + Hit h1 = triple_hits[0]; + Hit h2 = triple_hits[1]; + Hit h3 = triple_hits[2]; - cout << h1.volumeID << "\t"; - cout << h2.volumeID << "\t"; + cout.precision(2); + cout << h1.posX / 10. << "\t" << h1.posY / 10. << "\t" << h1.posZ / 10. << "\t"; + cout.precision(1); + cout << h1.time << "\t"; - cout.precision(2); - cout << h1.edep << "\t"; - cout << h2.edep << "\t"; + cout.precision(2); + cout << h2.posX / 10. << "\t" << h2.posY / 10. << "\t" << h2.posZ / 10. << "\t"; + cout.precision(1); + cout << h2.time << "\t"; - cout << event_analysis::verify_type_of_triple_coincidence(triple_hits) << "\t"; + cout << h1.volumeID << "\t"; + cout << h2.volumeID << "\t"; - cout.precision(2); - cout << h1.sourcePosX/10. << "\t" << h1.sourcePosY/10. << "\t" << h1.sourcePosZ/10. << "\t"; - cout << h2.sourcePosX/10. << "\t" << h2.sourcePosY/10. << "\t" << h2.sourcePosZ/10. << endl; + cout.precision(2); + cout << h1.edep << "\t"; + cout << h2.edep << "\t"; - cout.precision(2); - cout << h3.posX/10. << "\t" << h3.posY/10. << "\t" << h3.posZ/10. << "\t"; - cout.precision(1); - cout << h3.time << "\t"; - cout << h3.volumeID << "\t"; - cout.precision(2); - cout << h3.edep << "\t"; - cout.precision(2); - cout << h3.sourcePosX/10. << "\t" << h3.sourcePosY/10. << "\t" << h3.sourcePosZ/10. << "\t"; + cout << event_analysis::verify_type_of_triple_coincidence(triple_hits) << "\t"; -} + cout.precision(2); + cout << h1.sourcePosX / 10. << "\t" << h1.sourcePosY / 10. << "\t" << h1.sourcePosZ / 10. << "\t"; + cout << h2.sourcePosX / 10. << "\t" << h2.sourcePosY / 10. << "\t" << h2.sourcePosZ / 10. << "\t"; + cout.precision(2); + cout << h3.posX / 10. << "\t" << h3.posY / 10. << "\t" << h3.posZ / 10. << "\t"; + cout.precision(1); + cout << h3.time << "\t"; + cout << h3.volumeID << "\t"; + cout.precision(2); + cout << h3.edep << "\t"; + cout.precision(2); + cout << h3.sourcePosX / 10. << "\t" << h3.sourcePosY / 10. << "\t" << h3.sourcePosZ / 10. << endl; + } } //================================================================================ // MAIN ANALYSIS FUNCTION //================================================================================ -void analyze_event(const vector &hits, const bool hits_are_singles, const bool triple_coincidence) +void analyze_event(const vector& hits, const bool hits_are_singles, const bool triple_coincidence) { - double COMPTON_E_TH = atof(getenv("GOJA_COMPTON_E_TH"))*1e3; - double PROMPT_E_TH = atof(getenv("GOJA_PROMPT_E_TH"))*1e3; + double COMPTON_E_TH = atof(getenv("GOJA_COMPTON_E_TH")) * 1e3; + double PROMPT_E_TH = atof(getenv("GOJA_PROMPT_E_TH")) * 1e3; int MAX_N = int(atof(getenv("GOJA_MAX_N"))); int MAX_N0 = int(atof(getenv("GOJA_MAX_N0"))); int N = 0; int N0 = 0; const auto averagingMethod = event_analysis::AveragingMethod(int(atof(getenv("GOJA_AVERAGING_METHOD")))); const string systemType = string(getenv("GOJA_SYSTEM_TYPE")); - - event_analysis::sort_hits(const_cast &>(hits), "TIME"); - + event_analysis::sort_hits(const_cast&>(hits), "TIME"); std::vector selected_hits; - if (hits_are_singles) { + if (hits_are_singles) + { std::tie(N0, N, selected_hits) = select_coincident_singles(hits, COMPTON_E_TH, systemType, averagingMethod); } - else { + else + { std::tie(N0, N, selected_hits) = event_analysis::select_coincident_hits(hits, COMPTON_E_TH); } // If number of selected singles/hits equals maximum number of events above the fixed energy threshold in the coincidence window // and number of all hits in the event is equal or smaller than maximum number of events above the noise energy threshold in the coincidence window // then print the coincidence: - if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); - //if (N>=2 and N<=MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); - - if (triple_coincidence) { - if (N==MAX_N and N0<=MAX_N0) {print_coincidences(selected_hits);} + //if (N == MAX_N and N0 <= MAX_N0) + // print_coincidences(selected_hits); + // if (N>=2 and N<=MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); + if (!triple_coincidence) + { + if (N == MAX_N and N0 <= MAX_N0) + { + print_coincidences(selected_hits); + } } - else { - if (N==MAX_N and N0<=MAX_N0) {print_triple_coincidences(selected_hits, PROMPT_E_TH);}; + else + { + cout << "3333" << endl; + if (N == MAX_N and N0 <= MAX_N0) + { + print_triple_coincidences(selected_hits, PROMPT_E_TH); + }; } - if (DEBUG) { + if (DEBUG) + { cout.setf(ios::fixed); cout.precision(1); - for (unsigned int i=0; i0) { + if (i > 0) + { cout.precision(1); - cout << "\t" << (hits[i].time-hits[i-1].time); + cout << "\t" << (hits[i].time - hits[i - 1].time); } - if (i==(hits.size()-1)) { + if (i == (hits.size() - 1)) + { cout.precision(1); - cout << "\t" << (hits[i].time-hits[0].time) << endl; + cout << "\t" << (hits[i].time - hits[0].time) << endl; cout << N << "\t" << N0 << endl; } cout << endl; } } - } void RunTests() @@ -471,17 +505,23 @@ void RunTests() Hit h1; h1.edep = 10; h1.time = 12; - h1.posX = 1; h1.posY = 2; h1.posZ = 3; + h1.posX = 1; + h1.posY = 2; + h1.posZ = 3; Hit h2; h2.edep = 80; h2.time = 22; - h2.posX = 2; h2.posY = 3; h2.posZ = 4; + h2.posX = 2; + h2.posY = 3; + h2.posZ = 4; Hit h3; h3.edep = 50; h3.time = 32; - h3.posX = 7; h3.posY = 4; h3.posZ = 5; + h3.posX = 7; + h3.posY = 4; + h3.posZ = 5; std::vector hits = {h1, h2, h3}; @@ -492,24 +532,25 @@ void RunTests() } auto single1 = merge_hits(hits, event_analysis::kCentroidWinnerNaivelyWeighted); - assert(std::abs(single1.time-22) < epsilon); // (12+22+32)/3 = 22 - assert(std::abs(single1.posX-3.333333333) < epsilon); // (1+2+7)/2 = 3.333333333 - assert(std::abs(single1.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 + assert(std::abs(single1.time - 22) < epsilon); // (12+22+32)/3 = 22 + assert(std::abs(single1.posX - 3.333333333) < epsilon); // (1+2+7)/2 = 3.333333333 + assert(std::abs(single1.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 auto single2 = merge_hits(hits, event_analysis::kCentroidWinnerEnergyWeighted); - assert(std::abs(single2.time-24.85714286) < epsilon); // (10*12+80*22+50*32)/140 = 24.85714286 - assert(std::abs(single2.posX-3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 - assert(std::abs(single2.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 + assert(std::abs(single2.time - 24.85714286) < epsilon); // (10*12+80*22+50*32)/140 = 24.85714286 + assert(std::abs(single2.posX - 3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 + assert(std::abs(single2.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 auto single3 = merge_hits(hits, event_analysis::kCentroidWinnerEnergyWeightedFirstTime); - assert(std::abs(single3.time-12) < epsilon); // min(12,22,32) = 12 - assert(std::abs(single3.posX-3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 - assert(std::abs(single3.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 + assert(std::abs(single3.time - 12) < epsilon); // min(12,22,32) = 12 + assert(std::abs(single3.posX - 3.714285714) < epsilon); // (10*1+80*2+50*7)/140 = 3.714285714 + assert(std::abs(single3.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 auto single4 = merge_hits(hits, event_analysis::kEnergyWinner); - assert(std::abs(single4.time-22) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 - assert(std::abs(single4.posX-2) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 - assert(single4.eventID == 10); // 10 is a value set manually - assert(single4.volumeID == 9); // 9 is a value set manually + assert(std::abs(single4.time - 22) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 + assert(std::abs(single4.posX - 2) < epsilon); // max(10,80,50) = 80 => single takes attributes of hit h2 + assert(single4.eventID == 10); // 10 is a value set manually + assert(single4.volumeID == 9); // 9 is a value set manually assert(std::abs(single4.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 } +} //end of namespace \ No newline at end of file diff --git a/goja/EventAnalysis.h b/goja/EventAnalysis.h index 1ccb577..b4e8aa8 100755 --- a/goja/EventAnalysis.h +++ b/goja/EventAnalysis.h @@ -64,8 +64,7 @@ Hit merge_hits(const std::vector &hits, const AveragingMethod winner); void print_coincidence(const Hit& h1, const Hit& h2); void print_coincidences(const std::vector& hits); void print_triple_coincidences(const std::vector& hits, double PROMPT_E_TH); - void analyze_event(std::vector &hits, bool hits_are_singles = true); - void analyze_event(const std::vector &hits, const bool hits_are_singles, const bool triple_coincidence); + void analyze_event(const std::vector &hits, const bool hits_are_singles = true, const bool triple_coincidence = false); };