Helper functions for PMT to photon production process matching

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 

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

src/util/CMakeLists.txt

@@ -5,6 +5,7 @@ add_library(util OBJECT
         src/GLG4StringUtil.cc
         src/GaussianRatioPDF.cc
         src/MultiChargeDist.cc
+        src/PMTMatching.cc
         src/PolygonOrientation.cc
         src/ReadFile.cc
         src/Sampling.cc

src/util/include/RAT/PMTMatching.hh

@@ -0,0 +1,51 @@
+#ifndef __RAT_PMTMatching__
+#define __RAT_PMTMatching__
+
+#include <TObject.h>
+#include <RAT/DS/Root.hh>
+#include <RAT/DS/MC.hh>
+#include <RAT/DS/EV.hh>
+#include <RAT/DS/PMT.hh>
+#include <RAT/DS/MCPMT.hh>
+#include <RAT/DS/MCPhoton.hh>
+
+namespace RAT{
+
+  /*
+   * 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<DS::MCPMT*> MCPMTsFromPMT(DS::PMT* pmt, DS::MC* mc);
+
+  // Returns all MCPhotons that contributed to PMT hit
+  std::vector<DS::MCPhoton*> MCPhotonsFromPMT(DS::PMT* pmt, DS::EV* trig, DS::MC* mc);
+  std::vector<DS::MCPhoton*> MCPhotonsFromPMT(DS::PMT* pmt, DS::Root* event);
+
+  // Returns key with maximum value in map
+  std::string MaxProcess(std::map<std::string, double> processMap);
+
+  // Returns the photon production process that contributed the most energy to the PMT hit
+  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::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<DS::MCPhoton*> mcphotons, std::string metric, bool rollupScint=true);
+
+}
+
+#endif

src/util/src/PMTMatching.cc

@@ -0,0 +1,147 @@
+#include <RAT/PMTMatching.hh>
+#include <iostream>
+
+namespace RAT{
+
+
+  // Returns all MCPMTs with matching IDs
+  std::vector<DS::MCPMT*> MCPMTsFromPMT(DS::PMT* pmt, DS::MC* mc) {
+    int id = pmt->GetID();
+    std::vector<DS::MCPMT*> 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<DS::MCPhoton*> 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<DS::MCPhoton*> matchedPhotons;
+
+    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++){
+        // 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));
+      }
+    }
+
+    return matchedPhotons;
+  }
+
+  std::vector<DS::MCPhoton*> 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<DS::MCPhoton*> empty;
+    return empty;
+  }
+
+
+  // Returns key with maximum value in map
+  std::string MaxProcess(std::map<std::string, double> 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::EV* trig, DS::MC* mc, bool rollupScint){
+    std::vector<DS::MCPhoton*> mcphotons = MCPhotonsFromPMT(pmt, trig, mc);
+    return ProcessFrom(mcphotons, "energy", rollupScint);
+  }
+
+  std::string ProcessFromEnergy(DS::PMT* pmt, DS::Root* event, bool rollupScint){
+    std::vector<DS::MCPhoton*> mcphotons = MCPhotonsFromPMT(pmt, event);
+    return ProcessFrom(mcphotons, "energy", rollupScint);
+  }
+
+  std::string ProcessFromEnergy(DS::MCPMT* pmt, bool rollupScint){
+    std::map<std::string, double> processEdepMap;
+    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::EV* trig, DS::MC* mc, bool rollupScint){
+    std::vector<DS::MCPhoton*> mcphotons = MCPhotonsFromPMT(pmt, trig, mc);
+    return ProcessFrom(mcphotons, "hits", rollupScint);
+  }
+
+  std::string ProcessFromNHit(DS::PMT* pmt, DS::Root* event, bool rollupScint){
+    std::vector<DS::MCPhoton*> mcphotons = MCPhotonsFromPMT(pmt, event);
+    return ProcessFrom(mcphotons, "hits", rollupScint);
+  }
+
+  std::string ProcessFromNHit(DS::MCPMT* pmt, bool rollupScint){
+   std::map<std::string, double> processHitMap;
+   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<DS::MCPhoton*> mcphotons, std::string metric, bool rollupScint){
+    if(mcphotons.size() == 0) return "unknown";
+    std::map<std::string, double> 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);
+  }
+
+}