From 101f17cedf2e58cb09122adf7c494a69b4323523 Mon Sep 17 00:00:00 2001 From: Tom Brooks Date: Tue, 8 Sep 2020 11:42:53 +0100 Subject: [PATCH 1/2] Add util functions for matching reconstructed PMT hits to truth --- src/util/CMakeLists.txt | 1 + src/util/include/RAT/PMTMatching.hh | 35 ++++++++++++ src/util/src/PMTMatching.cc | 87 +++++++++++++++++++++++++++++ 3 files changed, 123 insertions(+) create mode 100644 src/util/include/RAT/PMTMatching.hh create mode 100644 src/util/src/PMTMatching.cc diff --git a/src/util/CMakeLists.txt b/src/util/CMakeLists.txt index 499ae966..2be15b3e 100644 --- a/src/util/CMakeLists.txt +++ b/src/util/CMakeLists.txt @@ -4,6 +4,7 @@ add_library(util OBJECT src/GLG4StringUtil.cc src/GaussianRatioPDF.cc + src/PMTMatching.cc src/MultiChargeDist.cc src/PolygonOrientation.cc src/ReadFile.cc diff --git a/src/util/include/RAT/PMTMatching.hh b/src/util/include/RAT/PMTMatching.hh new file mode 100644 index 00000000..ac1720ec --- /dev/null +++ b/src/util/include/RAT/PMTMatching.hh @@ -0,0 +1,35 @@ +#ifndef __RAT_PMTMatching__ +#define __RAT_PMTMatching__ + +#include +#include +#include +#include +#include + +namespace RAT{ + + /* + * Match between trigger PMT and MCPMT/MCPhoton objects + */ + + // Returns all MCPMTs with matching IDs + std::vector MCPMTsFromPMT(DS::PMT& pmt, DS::MC& mc); + + // Returns all MCPhotons that contributed to PMT hit + std::vector MCPhotonsFromPMT(DS::PMT& pmt, DS::MC& mc); + + // Returns key with maximum value in map + std::string MaxProcess(std::map processMap); + + // Returns the photon production process that contributed the most energy to the PMT hit + std::string ProcessFromEnergy(DS::PMT& pmt, DS::MC& mc, bool rollupScint=true); + std::string ProcessFromEnergy(DS::MCPMT& pmt, bool rollupScint=true); + + // Returns the photon production process that contributed the most photons to the PMT hit + std::string ProcessFromNHit(DS::PMT& pmt, DS::MC& mc, bool rollupScint=true); + std::string ProcessFromNHit(DS::MCPMT& pmt, bool rollupScint=true); + +} + +#endif diff --git a/src/util/src/PMTMatching.cc b/src/util/src/PMTMatching.cc new file mode 100644 index 00000000..7172ee49 --- /dev/null +++ b/src/util/src/PMTMatching.cc @@ -0,0 +1,87 @@ +#include + +namespace RAT{ + + // Returns all MCPMTs with matching IDs + std::vector MCPMTsFromPMT(DS::PMT& pmt, DS::MC& mc) { + int id = pmt.GetID(); + double time = pmt.GetTime(); + std::vector matchedPmts; + for(int i = 0; i < mc.GetMCPMTCount(); i++){ + int mcId = mc.GetMCPMT(i)->GetID(); + if(id != mcId) continue; + matchedPmts.push_back(mc.GetMCPMT(i)); + } + return matchedPmts; + } + + // Returns all MCPhotons that contributed to PMT hit + std::vector MCPhotonsFromPMT(DS::PMT& pmt, DS::MC& mc){ + int id = pmt.GetID(); + double time = pmt.GetTime(); + std::vector matchedPhotons; + for(int i = 0; i < mc.GetMCPMTCount(); i++){ + int mcId = mc.GetMCPMT(i)->GetID(); + if(id != mcId) continue; + if(mc.GetMCPMT(i)->GetMCPhotonCount() == 0) continue; + for(int j = 0; j < mc.GetMCPMT(i)->GetMCPhotonCount(); j++){ + double mcTime = mc.GetMCPMT(i)->GetMCPhoton(j)->GetFrontEndTime(); + if(std::abs(time - mcTime) > 800) continue; + matchedPhotons.push_back(mc.GetMCPMT(i)->GetMCPhoton(j)); + } + } + return matchedPhotons; + } + + // Returns key with maximum value in map + std::string MaxProcess(std::map processMap){ + double max = 0; + std::string maxProcess = "unknown"; + for(auto const& kv : processMap){ + if(kv.second > max){ + max = kv.second; + maxProcess = kv.first; + } + } + return maxProcess; + } + + // Returns the photon production process that contributed the most energy to the PMT hit + std::string ProcessFromEnergy(DS::PMT& pmt, DS::MC& mc){ + std::vector mcphotons = MCPhotonsFromPMT(pmt, mc); + if(mcphotons.size() == 0) return "unknown"; + std::map processEdepMap; + for(auto& mcphoton : mcphotons){ + processEdepMap[mcphoton->GetProcess()] += mcphoton->GetCharge(); + } + return MaxProcess(processEdepMap); + } + + std::string ProcessFromEnergy(DS::MCPMT& pmt){ + std::map processEdepMap; + for(int i = 0; i < pmt.GetMCPhotonCount(); i++){ + processEdepMap[pmt.GetMCPhoton(i)->GetProcess()] += pmt.GetMCPhoton(i)->GetCharge(); + } + return MaxProcess(processEdepMap); + } + + // Returns the photon production process that contributed the most photons to the PMT hit + std::string ProcessFromNHit(DS::PMT& pmt, DS::MC& mc){ + std::vector mcphotons = MCPhotonsFromPMT(pmt, mc); + if(mcphotons.size() == 0) return "unknown"; + std::map processHitMap; + for(auto& mcphoton : mcphotons){ + processHitMap[mcphoton->GetProcess()] += 1; + } + return MaxProcess(processHitMap); + } + + std::string ProcessFromNHit(DS::MCPMT& pmt){ + std::map processHitMap; + for(int i = 0; i < pmt.GetMCPhotonCount(); i++){ + processHitMap[pmt.GetMCPhoton(i)->GetProcess()] += 1; + } + return MaxProcess(processHitMap); + } + +} From 6d1b36eeb01ce9fe9eb56728a8d94e4ba23bf703 Mon Sep 17 00:00:00 2001 From: Tom Brooks Date: Wed, 9 Sep 2020 16:27:19 +0100 Subject: [PATCH 2/2] Add helper functions for matching PMTs to photon production processes --- src/CMakeLists.txt | 1 + src/ds/CMakeLists.txt | 3 +- src/util/CMakeLists.txt | 2 +- src/util/include/RAT/PMTMatching.hh | 30 +++++-- src/util/src/PMTMatching.cc | 134 ++++++++++++++++++++-------- 5 files changed, 124 insertions(+), 46 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 14bc6e22..ac991e78 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -54,6 +54,7 @@ add_library(RATEvent SHARED stlplus/src/string_utilities.cc stlplus/src/stringio.cc stlplus/src/textio.cc + util/src/PMTMatching.cc util/src/ReadFile.cc) target_include_directories(RATEvent PUBLIC ${RATPAC_INCDIR}) target_link_libraries(RATEvent PUBLIC diff --git a/src/ds/CMakeLists.txt b/src/ds/CMakeLists.txt index 976a16f5..0eb65107 100644 --- a/src/ds/CMakeLists.txt +++ b/src/ds/CMakeLists.txt @@ -50,8 +50,9 @@ root_generate_dictionary(G__RATDict RAT/Log.hh RAT/ObjInt.hh RAT/ObjDbl.hh + RAT/PMTMatching.hh LINKDEF include/RAT/DS/LinkDef.hh - DEPENDENCIES core ds db io) + DEPENDENCIES core ds db io util) add_library(RATDict OBJECT G__RATDict.cxx) target_include_directories(RATDict PUBLIC diff --git a/src/util/CMakeLists.txt b/src/util/CMakeLists.txt index 2be15b3e..d639af48 100644 --- a/src/util/CMakeLists.txt +++ b/src/util/CMakeLists.txt @@ -4,8 +4,8 @@ add_library(util OBJECT src/GLG4StringUtil.cc src/GaussianRatioPDF.cc - src/PMTMatching.cc src/MultiChargeDist.cc + src/PMTMatching.cc src/PolygonOrientation.cc src/ReadFile.cc src/Sampling.cc diff --git a/src/util/include/RAT/PMTMatching.hh b/src/util/include/RAT/PMTMatching.hh index ac1720ec..290ef89c 100644 --- a/src/util/include/RAT/PMTMatching.hh +++ b/src/util/include/RAT/PMTMatching.hh @@ -2,7 +2,9 @@ #define __RAT_PMTMatching__ #include +#include #include +#include #include #include #include @@ -10,25 +12,39 @@ namespace RAT{ /* - * Match between trigger PMT and MCPMT/MCPhoton objects + * Helper functions for matching between triggered PMT hits and + * MCPMT/MCPhoton objects, and matching PMT/MCPMT hits to the + * dominant photon production process by either total energy or + * total number of photons. + * Can be used in ROOT macro: + * > for(int i = 0; i < event->GetMC()->GetMCPMTCount(); i++){ + * > RAT::DS* mcpmt = event->GetMC()->GetMCPMT(i); + * > std::string process = RAT::ProcessFromEnergy(mcpmt); + * > } */ // Returns all MCPMTs with matching IDs - std::vector MCPMTsFromPMT(DS::PMT& pmt, DS::MC& mc); + std::vector MCPMTsFromPMT(DS::PMT* pmt, DS::MC* mc); // Returns all MCPhotons that contributed to PMT hit - std::vector MCPhotonsFromPMT(DS::PMT& pmt, DS::MC& mc); + std::vector MCPhotonsFromPMT(DS::PMT* pmt, DS::EV* trig, DS::MC* mc); + std::vector MCPhotonsFromPMT(DS::PMT* pmt, DS::Root* event); // Returns key with maximum value in map std::string MaxProcess(std::map processMap); // Returns the photon production process that contributed the most energy to the PMT hit - std::string ProcessFromEnergy(DS::PMT& pmt, DS::MC& mc, bool rollupScint=true); - std::string ProcessFromEnergy(DS::MCPMT& pmt, bool rollupScint=true); + std::string ProcessFromEnergy(DS::PMT* pmt, DS::EV* trig, DS::MC* mc, bool rollupScint=true); + std::string ProcessFromEnergy(DS::PMT* pmt, DS::Root* event, bool rollupScint=true); + std::string ProcessFromEnergy(DS::MCPMT* pmt, bool rollupScint=true); // Returns the photon production process that contributed the most photons to the PMT hit - std::string ProcessFromNHit(DS::PMT& pmt, DS::MC& mc, bool rollupScint=true); - std::string ProcessFromNHit(DS::MCPMT& pmt, bool rollupScint=true); + std::string ProcessFromNHit(DS::PMT* pmt, DS::EV* trig, DS::MC* mc, bool rollupScint=true); + std::string ProcessFromNHit(DS::PMT* pmt, DS::Root* event, bool rollupScint=true); + std::string ProcessFromNHit(DS::MCPMT* pmt, bool rollupScint=true); + + // Returns the photon production process that contributed the most of a certain metric (energy or hits) to a PMT hit + std::string ProcessFrom(std::vector mcphotons, std::string metric, bool rollupScint=true); } diff --git a/src/util/src/PMTMatching.cc b/src/util/src/PMTMatching.cc index 7172ee49..d95bff6c 100644 --- a/src/util/src/PMTMatching.cc +++ b/src/util/src/PMTMatching.cc @@ -1,38 +1,66 @@ #include +#include namespace RAT{ + // Returns all MCPMTs with matching IDs - std::vector MCPMTsFromPMT(DS::PMT& pmt, DS::MC& mc) { - int id = pmt.GetID(); - double time = pmt.GetTime(); + std::vector MCPMTsFromPMT(DS::PMT* pmt, DS::MC* mc) { + int id = pmt->GetID(); std::vector matchedPmts; - for(int i = 0; i < mc.GetMCPMTCount(); i++){ - int mcId = mc.GetMCPMT(i)->GetID(); + + for(int i = 0; i < mc->GetMCPMTCount(); i++){ + int mcId = mc->GetMCPMT(i)->GetID(); if(id != mcId) continue; - matchedPmts.push_back(mc.GetMCPMT(i)); + matchedPmts.push_back(mc->GetMCPMT(i)); } + return matchedPmts; } + // Returns all MCPhotons that contributed to PMT hit - std::vector MCPhotonsFromPMT(DS::PMT& pmt, DS::MC& mc){ - int id = pmt.GetID(); - double time = pmt.GetTime(); + std::vector MCPhotonsFromPMT(DS::PMT* pmt, DS::EV* trig, DS::MC* mc){ + int id = pmt->GetID(); + // Get the global time of the PMT hit + double time = pmt->GetTime() + trig->GetCalibratedTriggerTime(); std::vector matchedPhotons; - for(int i = 0; i < mc.GetMCPMTCount(); i++){ - int mcId = mc.GetMCPMT(i)->GetID(); + + for(int i = 0; i < mc->GetMCPMTCount(); i++){ + int mcId = mc->GetMCPMT(i)->GetID(); + // Find MCPMTs with the same ID if(id != mcId) continue; - if(mc.GetMCPMT(i)->GetMCPhotonCount() == 0) continue; - for(int j = 0; j < mc.GetMCPMT(i)->GetMCPhotonCount(); j++){ - double mcTime = mc.GetMCPMT(i)->GetMCPhoton(j)->GetFrontEndTime(); + if(mc->GetMCPMT(i)->GetMCPhotonCount() == 0) continue; + for(int j = 0; j < mc->GetMCPMT(i)->GetMCPhotonCount(); j++){ + // Return MCPhotons within trigger window + // TODO this is using the LessSimpleDAQ trigger window, would be good to have a way to set this from the event + double mcTime = mc->GetMCPMT(i)->GetMCPhoton(j)->GetFrontEndTime(); if(std::abs(time - mcTime) > 800) continue; - matchedPhotons.push_back(mc.GetMCPMT(i)->GetMCPhoton(j)); + matchedPhotons.push_back(mc->GetMCPMT(i)->GetMCPhoton(j)); } } + return matchedPhotons; } + std::vector MCPhotonsFromPMT(DS::PMT* pmt, DS::Root* event){ + // Get the MC and trigger info by matching the PMTs + DS::MC* mc = event->GetMC(); + for(int i = 0; i < event->GetEVCount(); i++){ + DS::EV* trig = event->GetEV(i); + for(int j = 0; j < trig->GetPMTCount(); j++){ + if(trig->GetPMT(j)->GetID() == pmt->GetID() && trig->GetPMT(j)->GetTime() == pmt->GetTime()){ + return MCPhotonsFromPMT(pmt, trig, mc); + } + } + } + + // If no matching PMTs found return an empty vector + std::vector empty; + return empty; + } + + // Returns key with maximum value in map std::string MaxProcess(std::map processMap){ double max = 0; @@ -46,42 +74,74 @@ namespace RAT{ return maxProcess; } + // Returns the photon production process that contributed the most energy to the PMT hit - std::string ProcessFromEnergy(DS::PMT& pmt, DS::MC& mc){ - std::vector mcphotons = MCPhotonsFromPMT(pmt, mc); - if(mcphotons.size() == 0) return "unknown"; - std::map processEdepMap; - for(auto& mcphoton : mcphotons){ - processEdepMap[mcphoton->GetProcess()] += mcphoton->GetCharge(); - } - return MaxProcess(processEdepMap); + std::string ProcessFromEnergy(DS::PMT* pmt, DS::EV* trig, DS::MC* mc, bool rollupScint){ + std::vector mcphotons = MCPhotonsFromPMT(pmt, trig, mc); + return ProcessFrom(mcphotons, "energy", rollupScint); + } + + std::string ProcessFromEnergy(DS::PMT* pmt, DS::Root* event, bool rollupScint){ + std::vector mcphotons = MCPhotonsFromPMT(pmt, event); + return ProcessFrom(mcphotons, "energy", rollupScint); } - std::string ProcessFromEnergy(DS::MCPMT& pmt){ + std::string ProcessFromEnergy(DS::MCPMT* pmt, bool rollupScint){ std::map processEdepMap; - for(int i = 0; i < pmt.GetMCPhotonCount(); i++){ - processEdepMap[pmt.GetMCPhoton(i)->GetProcess()] += pmt.GetMCPhoton(i)->GetCharge(); + for(int i = 0; i < pmt->GetMCPhotonCount(); i++){ + if(rollupScint && pmt->GetMCPhoton(i)->GetProcess()=="Reemission"){ + processEdepMap["Scintillation"] += pmt->GetMCPhoton(i)->GetCharge(); + } + else{ + processEdepMap[pmt->GetMCPhoton(i)->GetProcess()] += pmt->GetMCPhoton(i)->GetCharge(); + } } return MaxProcess(processEdepMap); } + // Returns the photon production process that contributed the most photons to the PMT hit - std::string ProcessFromNHit(DS::PMT& pmt, DS::MC& mc){ - std::vector mcphotons = MCPhotonsFromPMT(pmt, mc); - if(mcphotons.size() == 0) return "unknown"; - std::map processHitMap; - for(auto& mcphoton : mcphotons){ - processHitMap[mcphoton->GetProcess()] += 1; - } - return MaxProcess(processHitMap); + std::string ProcessFromNHit(DS::PMT* pmt, DS::EV* trig, DS::MC* mc, bool rollupScint){ + std::vector mcphotons = MCPhotonsFromPMT(pmt, trig, mc); + return ProcessFrom(mcphotons, "hits", rollupScint); + } + + std::string ProcessFromNHit(DS::PMT* pmt, DS::Root* event, bool rollupScint){ + std::vector mcphotons = MCPhotonsFromPMT(pmt, event); + return ProcessFrom(mcphotons, "hits", rollupScint); } - std::string ProcessFromNHit(DS::MCPMT& pmt){ + std::string ProcessFromNHit(DS::MCPMT* pmt, bool rollupScint){ std::map processHitMap; - for(int i = 0; i < pmt.GetMCPhotonCount(); i++){ - processHitMap[pmt.GetMCPhoton(i)->GetProcess()] += 1; + for(int i = 0; i < pmt->GetMCPhotonCount(); i++){ + if(rollupScint && pmt->GetMCPhoton(i)->GetProcess()=="Reemission"){ + processHitMap["Scintillation"] += 1; + } + else{ + processHitMap[pmt->GetMCPhoton(i)->GetProcess()] += 1; + } } return MaxProcess(processHitMap); } + + // Returns the photon production process that contributed the most of a certain metric (energy or hits) to a PMT hit + std::string ProcessFrom(std::vector mcphotons, std::string metric, bool rollupScint){ + if(mcphotons.size() == 0) return "unknown"; + std::map processEdepMap; + for(auto& mcphoton : mcphotons){ + double value = 1; + if(metric == "energy"){ + value = mcphoton->GetCharge(); + } + if(rollupScint && mcphoton->GetProcess()=="Reemission"){ + processEdepMap["Scintillation"] += value; + } + else{ + processEdepMap[mcphoton->GetProcess()] += value; + } + } + return MaxProcess(processEdepMap); + } + }