Merge current changes

src/Params.h

@@ -44,6 +44,9 @@ namespace cafmaker
     fhicl::OptionalAtom<std::string> sandRecoFile  { fhicl::Name{"SANDRecoFile"}, fhicl::Comment("Input SAND reco .root file") };
     fhicl::OptionalAtom<std::string> minervaRecoFile  { fhicl::Name{"MINERVARecoFile"}, fhicl::Comment("Input MINERVA reco .root file") };
 
+    // fixme: temporary hack for data.  should be removed when interface to IFDB is complete
+    fhicl::OptionalAtom<std::string> POTFile { fhicl::Name{"POTFile"}, fhicl::Comment("Input txt file with beam spill information") };
+
     // this is optional by way of the default value. Will result in an extra output file if enabled
     fhicl::Atom<bool> makeFlatCAF { fhicl::Name{"MakeFlatCAF"}, fhicl::Comment("Make 'flat' CAF in addition to structured CAF?"), true };
 
@@ -55,6 +58,9 @@ namespace cafmaker
     // 100 us is default
     fhicl::Atom<unsigned int>  trigMatchDT { fhicl::Name("TriggerMatchDeltaT"), fhicl::Comment("Maximum time difference, in ns, between triggers to be considered a match"), 100000 };
 
+    // 0.1 s is default
+    fhicl::Atom<float>  beamMatchDT { fhicl::Name("BeamMatchDeltaT"), fhicl::Comment("Maximum time difference, in s, between triggers and beam"), 0.1 };
+
     // options are VERBOSE, DEBUG, INFO, WARNING, ERROR, FATAL
     fhicl::Atom<std::string> verbosity { fhicl::Name("Verbosity"), fhicl::Comment("Verbosity level of output"), "WARNING" };
   };

src/makeCAF.C

@@ -1,5 +1,11 @@
 #include <cstdio>
 #include <numeric>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <map>
+#include <cmath>
+#include <algorithm>
 
 #include "boost/program_options/options_description.hpp"
 #include "boost/program_options/parsers.hpp"
@@ -170,7 +176,7 @@ std::vector<std::unique_ptr<cafmaker::IRecoBranchFiller>> getRecoFillers(const c
 // -------------------------------------------------
 bool doTriggersMatch(const cafmaker::Trigger& t1, const cafmaker::Trigger& t2, unsigned int dT)
 {
-  return (t1.triggerTime_s == t2.triggerTime_s && abs(int(t1.triggerTime_ns - t2.triggerTime_ns)) < dT);
+  return ( (std::max(t1.triggerTime_s, t2.triggerTime_s) - std::min(t1.triggerTime_s, t2.triggerTime_s)) * 10000000 + std::max(t1.triggerTime_ns, t2.triggerTime_ns) - std::min(t1.triggerTime_ns, t2.triggerTime_ns) ) < dT;
 }
 
 struct triggerTimeCmp
@@ -324,6 +330,104 @@ buildTriggerList(std::map<const cafmaker::IRecoBranchFiller*, std::deque<cafmake
   return ret;
 }
 
+// -------------------------------------------------
+//Temporary hack to fill beam POT info for data
+
+//Get the spills from an input text file
+using BeamSpills = std::map<double, double>;
+
+bool loadBeamSpills(const std::string& filename, BeamSpills& beam_spills) {
+  std::ifstream file(filename);
+  if (!file.is_open()) {
+    std::cerr << "Could not open the file " << filename << "\n";
+    return false;
+  }
+
+  std::string line;
+  while (std::getline(file, line)) {
+    std::istringstream iss(line);
+    double time;
+    double pot;
+    if (iss >> time >> pot)
+      beam_spills[time] = pot;
+  }
+  file.close();
+  return true;
+}
+
+double GetTriggerTime(const cafmaker::Trigger& trigger) {
+    return trigger.triggerTime_s + 1e-9 * trigger.triggerTime_ns;
+}
+
+//Return pot from text file if all of the trigger times match the beam times within a certain dt, else fill given pot
+double getPOT(const cafmaker::Params& par, std::vector<std::pair<const cafmaker::IRecoBranchFiller*, cafmaker::Trigger>>& groupedTrigger, int ii) {
+  std::string potFile;
+  BeamSpills beam_spills;
+
+  double pot = 0.0;
+
+  if (par().cafmaker().POTFile(potFile) && loadBeamSpills(potFile, beam_spills)) { //If there is a POT file in config that is readable, check if all trigger times match any of the beam times
+    auto it = std::find_if(beam_spills.begin(), beam_spills.end(),                
+                           [par, &groupedTrigger](const auto& spill) {
+      return std::all_of(groupedTrigger.cbegin(), groupedTrigger.cend(),
+                         [par, &spill](const auto& groupedTrigger) {
+        return std::abs(GetTriggerTime(groupedTrigger.second) - spill.first) < par().cafmaker().beamMatchDT(); //fixme: shouldn't be abs after 2x2 trigger times are fixed
+        });
+      });
+
+    if (it != beam_spills.end()) {
+      pot = it->second;
+    }     
+    else { //Check if any of the triggers match the beam time if there is no beam time match in the previous search
+      bool any_matched = false;
+      std::vector<std::pair<const cafmaker::IRecoBranchFiller*, cafmaker::Trigger>> matched_triggers;
+      std::vector<std::pair<const cafmaker::IRecoBranchFiller*, cafmaker::Trigger>> unmatched_triggers;
+
+      for (auto trig : groupedTrigger) {
+        bool matched = std::any_of(beam_spills.cbegin(), beam_spills.cend(),
+                                   [par, &trig](const auto& spill) {
+          return std::abs(GetTriggerTime(trig.second) - spill.first) < par().cafmaker().beamMatchDT(); //fixme: shouldn't be abs after 2x2 trigger times are fixed
+        });
+
+        if (matched) {
+          any_matched = true;
+          matched_triggers.push_back(trig);
+        } 
+        else {
+          unmatched_triggers.push_back(trig);
+        }
+      }
+
+      auto LOG = [&]() -> const cafmaker::Logger & { return cafmaker::LOG_S("Beam spill matching"); };
+      std::stringstream log_message;
+
+      if (any_matched) { //If only some of the trigger times in a grouped trigger match the beam spill, abort!
+         log_message << "Only some triggers match beam spill for trigger group " << ii << ":\n"
+                       << "Matched triggers: \n";
+        for (auto trig : matched_triggers)
+          log_message << std::fixed << trig.first->GetName() << " " << GetTriggerTime(trig.second) << "\n";
+
+        log_message << "Unmatched triggers: \n";
+        for (auto trig : unmatched_triggers)
+          log_message << std::fixed << trig.first->GetName() << " " << GetTriggerTime(trig.second) << "\n";
+
+        LOG().ERROR() << log_message.str() << "\n";
+        std::abort();
+      } 
+      else { //If none of the trigger times match give a warning
+        log_message << "No matching spill found for trigger group " << ii << " with triggers: \n";
+        for (auto trig : groupedTrigger)
+          log_message << std::fixed << trig.first->GetName() << " " << GetTriggerTime(trig.second) << "\n";
+        LOG().WARNING() << log_message.str() << "\n";
+      } 
+    } 
+  } 
+  else { //else get POT from config
+    pot = par().runInfo().POTPerSpill() * 1e13;
+  }
+
+  return pot;
+}
 // -------------------------------------------------
 // main loop function
 void loop(CAF &caf,
@@ -373,19 +477,20 @@ void loop(CAF &caf,
     // reset (the default constructor initializes its variables)
     caf.setToBS();
 
-    double pot = par().runInfo().POTPerSpill() * 1e13;
-    if (std::isnan(caf.pot))
-      caf.pot = 0;
-    caf.pot += pot;
-    caf.sr.beam.pulsepot = pot;
-    caf.sr.beam.ismc = true;  // when we have data, we won't be able to use this "POT from config" approach anyway
 
     // hand off to the correct reco filler(s).
     for (const auto & fillerTrigPair : groupedTriggers[ii])
     {
       cafmaker::LOG_S("loop()").INFO() << "Global trigger idx : " << ii << ", reco filler: '" << fillerTrigPair.first->GetName() << "', reco trigger eventID: " << fillerTrigPair.second.evtID << "\n";
       fillerTrigPair.first->FillRecoBranches(fillerTrigPair.second, caf.sr, par, &truthMatcher);
     }
+    //Fill POT
+    double pot = getPOT(par, groupedTriggers[ii], ii);
+    if (std::isnan(caf.pot))
+      caf.pot = 0;
+    caf.pot += pot;
+    caf.sr.beam.pulsepot = pot;
+    caf.sr.beam.ismc = par().cafmaker().POTFile.hasValue();  // fixme: when we have proper IFDB interface, should use the same mechanism as however we decide when to use that
 
     caf.fill();
   }