Feature/ingratta add beam quality info

src/beam/IFBeam.cxx

@@ -9,6 +9,7 @@
 #include "util/Logger.h"
 #include "util/IFBeamUtils.h"
 #include <curl/curl.h>
+#include <iostream>
 #include <regex>
 #include <sstream>
 #include <limits>
@@ -20,12 +21,12 @@ namespace cafmaker
 {
 
 
-  IFBeam::IFBeam(const std::vector<TriggerGroup>& groupedTriggers, bool is_data) {
+  IFBeam::IFBeam(const cafmaker::Params& par, const std::vector<TriggerGroup>& groupedTriggers, bool is_data) {
       if (is_data)
-  	loadBeamSpills(groupedTriggers); // Load all beam spills upon instantiation if data
+  	loadBeamSpills(par, groupedTriggers); // Load all beam spills upon instantiation if data
   }
 
-  std::string IFBeam::createUrl(const std::string& min_time_iso, const std::string& max_time_iso) {
+  std::string IFBeam::createUrl(const std::string potDevice, const std::string& min_time_iso, const std::string& max_time_iso) {
       std::ostringstream url_stream;
       url_stream << "https://dbdata3vm.fnal.gov:9443/ifbeam/data/data?v=" << potDevice
                  << "&e=" << "e,a9"
@@ -43,42 +44,23 @@ namespace cafmaker
           int exponent = std::stoi(match.str(1));
           return std::pow(10, exponent);
       } else {
-          cafmaker::LOG_S("Fetch beam information").WARNING() << "Unknown unit in beam database: " << unit << "\n";
           return 1.0;
       }
   }
-
-  void IFBeam::loadBeamSpills(const std::vector<TriggerGroup>& groupedTriggers) { //todo: this should return other information as well like horn current, position, etc., querying all devices and storing in a map
-      beamSpills.clear();
-      double dt = 5.0; //time window to query before and after first and last spill, respectively
+
+  IFBeam::BeamInfo IFBeam::retrieveInfoFromDataBase(const std::string url) {
+
+      cafmaker::LOG_S("Fetchig beam info from url : ").VERBOSE() << url << "\n";
+
       double ms_to_s = 1e-3;
-      double min_time = std::numeric_limits<double>::max();
-      double max_time = std::numeric_limits<double>::lowest();
-      for (const auto& group : groupedTriggers) {
-          for (const auto& trig : group) {
-              if (!trig.first->IsBeamTrigger(trig.second.triggerType)) continue;
-              double trigger_time = util::getTriggerTime(trig.second);
 
-              min_time = std::min(min_time, trigger_time - dt);
-              max_time = std::max(max_time, trigger_time + dt);
-          }
-      }
-      if (min_time == std::numeric_limits<double>::max() || max_time == std::numeric_limits<double>::lowest())
-      {
-        std::cerr<<"No beam trigger loaded, are you sure you want to use IFBeam database? \n";
-        return;
-      }
-      std::string min_time_iso = util::toISO8601(min_time);
-      std::string max_time_iso = util::toISO8601(max_time);
-      std::string url = createUrl(min_time_iso, max_time_iso);
-
-
+      BeamInfo data;
+
       CURL* curl;
       CURLcode res;
-      std::string readBuffer;
-
-      std::cout << "Fetching beam information from: " << url << "\n";
+      std::string readBuffer; 
       curl = curl_easy_init();
+
       if (curl) {
           curl_easy_setopt(curl, CURLOPT_URL, url.c_str());
           curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, util::WriteCallback);
@@ -94,38 +76,127 @@ namespace cafmaker
               for (const auto& spill : json_data["rows"]) {
                   double time = spill["clock"].get<double>() * ms_to_s;
                   std::string unit = spill["units"];
-                  double pot = spill["value"].get<double>() * unitToFactor(unit);
-                  beamSpills[time] = pot;
+                  auto& c_array = spill["c"];
+
+                  std::vector<double> values;
+                  for (const auto& c_entry : c_array) 
+                  {
+                    if (c_entry["v"].is_number())
+                    {
+                      double this_v = c_entry["v"].get<double>() * unitToFactor(unit);
+                      values.push_back(this_v);
+                    }
+                  }
+                  data[time] = values;
               }
           } catch (const json::exception& e) {
               cafmaker::LOG_S("Fetch beam information").ERROR() << "Failed to parse JSON data: " << e.what() << "\n";
   	    std::abort();
           }
       }
+      return data;
   }
-
-  double IFBeam::getPOT(const cafmaker::Params& par, const TriggerGroup& groupedTrigger, int ii) {
-      double pot = 0.0;
-      if (!(groupedTrigger.front().first->IsBeamTrigger(groupedTrigger.front().second.triggerType))) return 0.0;
-      auto it = std::find_if(beamSpills.begin(), beamSpills.end(),
+
+  void IFBeam::loadBeamSpills(const cafmaker::Params& par, const std::vector<TriggerGroup>& groupedTriggers) { 
+
+      double dt = 5.0; //time window to query before and after first and last spill, respectively
+      double min_time = std::numeric_limits<double>::max();
+      double max_time = std::numeric_limits<double>::lowest();
+      for (const auto& group : groupedTriggers) {
+          for (const auto& trig : group) {
+              if (!trig.first->IsBeamTrigger(trig.second.triggerType)) continue;
+              double trigger_time = util::getTriggerTime(trig.second);
+
+              min_time = std::min(min_time, trigger_time - dt);
+              max_time = std::max(max_time, trigger_time + dt);
+          }
+      }
+      if (min_time == std::numeric_limits<double>::max() || max_time == std::numeric_limits<double>::lowest())
+      {
+        std::cerr<<"No beam trigger loaded, are you sure you want to use IFBeam database? \n";
+        return;
+      }
+      std::string min_time_iso = util::toISO8601(min_time);
+      std::string max_time_iso = util::toISO8601(max_time);
+
+      for (auto& device : deviceMap)
+      {
+        device.second.clear();
+        std::string url_device = createUrl(device.first, min_time_iso, max_time_iso);
+        device.second = retrieveInfoFromDataBase(url_device);
+      }
+
+      // see: https://cdcvs.fnal.gov/redmine/projects/novaart/repository/entry/trunk/IFDBSpillInfo/IFDBSpillInfo_module.cc#L685
+      // the horn current is calculated as I_horn = (E:NSLINA-(+0.01))/0.9951 + (E:NSLINB-(-0.14))/0.9957 + (E:NSLINC-(-0.05))/0.9965 + (E:NSLIND-(-0.07))/0.9945
+
+      double currentA=0.0, currentB=0.0, currentC=0.0, currentD=0.0;
+
+      for(const auto& pairA : deviceMap["E:NSLINA"]) {
+        auto timeA = pairA.first;
+        auto currentA = pairA.second.at(0); 
+
+        if(currentA == 0.) continue;
+
+        for(const auto& pairB : deviceMap["E:NSLINB"]) {
+          auto timeB = pairB.first;
+          if (abs(timeA - timeB) <= par().cafmaker().beamMatchDT()) {
+            currentB = pairB.second.at(0);
+            break;
+          }
+        }
+
+        if(currentB == 0.) continue;
+
+        for(const auto& pairC : deviceMap["E:NSLINC"]) {
+          auto timeC = pairC.first;
+          if (abs(timeA - timeC) <= par().cafmaker().beamMatchDT()) 
+          {
+            currentC = pairC.second.at(0);
+            break;
+          }
+        }
+
+        if(currentC == 0.) continue;
+
+        for(const auto& pairD : deviceMap["E:NSLIND"]) {
+          auto timeD = pairD.first;
+          if (abs(timeA - timeD) <= par().cafmaker().beamMatchDT()){
+            currentD = pairD.second.at(0);
+            break;
+          }
+        }
+
+        if(currentD == 0.) continue;
+
+        deviceMap["E:NSLIN"][timeA] = {(currentA-0.01)/0.9951 + (currentB+0.14)/0.9957 + (currentC+0.05)/0.9965 + (currentD+0.07)/0.9945};
+      }
+
+  }
+
+  std::vector<double> IFBeam::getData(const cafmaker::Params& par, const TriggerGroup& groupedTrigger, int ii, const BeamInfo& data) {
+
+      std::vector<double> values;
+      if (!(groupedTrigger.front().first->IsBeamTrigger(groupedTrigger.front().second.triggerType))) return {0.0};
+
+      auto it = std::find_if(data.begin(), data.end(),
           [par, &groupedTrigger](const auto& spill) {
               return std::all_of(groupedTrigger.cbegin(), groupedTrigger.cend(),
                   [par, &spill](const auto& groupedTrigger) {
                       return std::abs(util::getTriggerTime(groupedTrigger.second) - spill.first) < par().cafmaker().beamMatchDT();
                   });
           });
-  
-      if (it != beamSpills.end()) {
-          pot = it->second;
+
+      if (it != data.end()) {
+          values = it->second;
       } 
       else {
           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) {
-              //Only fetch pot for beam trigger. 
-              if (!trig.first->IsBeamTrigger(trig.second.triggerType)) return 0;
-              bool matched = std::any_of(beamSpills.cbegin(), beamSpills.cend(),
+              //Only fetch values for beam trigger. 
+              if (!trig.first->IsBeamTrigger(trig.second.triggerType)) return {0.0};
+              bool matched = std::any_of(data.cbegin(), data.cend(),
                   [par, &trig](const auto& spill) {
                       return std::abs(util::getTriggerTime(trig.second) - spill.first) < par().cafmaker().beamMatchDT();
                   });
@@ -160,6 +231,16 @@ namespace cafmaker
           }
       }
 
-     return pot;
+     return values;
+  }
+
+  std::vector<double> IFBeam::getData(const cafmaker::Params& par, const TriggerGroup& groupedTrigger, int ii, const std::string& deviceName) {
+
+    if(deviceMap.count(deviceName)>0)
+      return getData(par, groupedTrigger, ii, deviceMap[deviceName]);
+
+    cafmaker::LOG_S("Trying retrieving info from unkwnown device ").WARNING() << deviceName << ", returning empty vector\n";
+
+    return {};
   }
 }

src/beam/IFBeam.h

@@ -9,7 +9,7 @@
 #define IFBEAM_H
 
 #include <Params.h>
-#include "reco/IRecoBranchFiller.h"
+#include "reco/IRecoBranchFiller.h"     
 #include <string>
 #include <map>
 #include <vector>
@@ -19,20 +19,43 @@ namespace cafmaker
 {
   class IFBeam {
   public:
-      using BeamSpills = std::map<double, double>;
+      using BeamInfo = std::map<double, std::vector<double>>;
+      using DeviceMap = std::map<const std::string, BeamInfo>;
       using TriggerGroup = std::vector<std::pair<const cafmaker::IRecoBranchFiller*, cafmaker::Trigger>>;
 
-      IFBeam(const std::vector<TriggerGroup>& groupedTriggers, bool is_data);   
-  
-      double getPOT(const cafmaker::Params& par, const TriggerGroup & groupedTrigger, int ii);
-
+      IFBeam(const cafmaker::Params& par, const std::vector<TriggerGroup>& groupedTriggers, bool is_data);   
+
+      std::vector<double> getData(const cafmaker::Params& par, const TriggerGroup& groupedTrigger, int ii, const BeamInfo& data);
+      std::vector<double> getData(const cafmaker::Params& par, const TriggerGroup& groupedTrigger, int ii,const std::string& deviceName);
 
   private:
-      const std::string potDevice = "E:TRTGTD";
-      BeamSpills beamSpills;
-
-      void loadBeamSpills(const std::vector<TriggerGroup>& groupedTriggers);
-      std::string createUrl(const std::string& min_time_iso, const std::string& max_time_iso);
+      DeviceMap deviceMap = {
+      //_________POT _______________________
+        {"E:TRTGTD", {}},
+        {"E:TOR101", {}},
+        {"E:TR101D", {}},
+      //_________Horn Current_______________
+        {"E:NSLINA", {}},
+        {"E:NSLINB", {}},
+        {"E:NSLINC", {}},
+        {"E:NSLIND", {}},
+        {"E:NSLIN", {}},
+      //________Horn polarity_______________
+        {"E:HRNDIR", {}},
+      //__________Beam Position_____________
+        {"E:HPTGT[]", {}},
+        {"E:HITGT[]", {}},
+        {"E:HP121[]", {}},
+        {"E:VPTGT[]", {}},
+        {"E:VITGT[]", {}},
+        {"E:VP121[]", {}},
+      //__________Beam Width________________
+        {"E:MTGTDS[]", {}}
+      };
+
+      BeamInfo retrieveInfoFromDataBase(const std::string url);
+      void loadBeamSpills(const cafmaker::Params& par, const std::vector<TriggerGroup>& groupedTriggers);
+      std::string createUrl(const std::string potDevice, const std::string& min_time_iso, const std::string& max_time_iso);
       double unitToFactor(const std::string& unit);
   };
 

src/makeCAF.C

@@ -413,7 +413,7 @@ void loop(CAF &caf,
   bool useIFBeam = false;
   if (ghepFilenames.empty() && edepsimFilename.empty() && !par().cafmaker().ForceDisableIFBeam()) useIFBeam = true;
 
-  cafmaker::IFBeam beamManager(groupedTriggers, useIFBeam); //initialize IFBeam manager if data and when IFBeam is not force disabled
+  cafmaker::IFBeam beamManager(par, groupedTriggers, useIFBeam); //initialize IFBeam manager if data and when IFBeam is not force disabled
   // Main event loop
   cafmaker::Progress progBar("Processing " + std::to_string(N - start) + " triggers");
   for( int ii = start; ii < start + N; ++ii )
@@ -444,21 +444,56 @@ void loop(CAF &caf,
       }
     }
 
-    //Fill POT
-    double pot = 0.0;
+    //Fill Beam info
+    std::vector<double> potTRTGTD;
+    std::vector<double> potTOR101;
+    std::vector<double> potTR101D;
+    std::vector<double> hornI;
+    std::vector<double> hornDir;
+    std::vector<double> horizontalposTGT;
+    std::vector<double> horizontalintTGT;
+    std::vector<double> horizontalpos121;
+    std::vector<double> verticalposTGT;
+    std::vector<double> verticalintTGT;
+    std::vector<double> verticalpos121;
+    std::vector<double> multiwireInfo;
     if (useIFBeam)
     {
-        pot = beamManager.getPOT(par, groupedTriggers[ii], ii);
+        potTRTGTD = beamManager.getData(par, groupedTriggers[ii], ii, "E:TRTGTD");
+        potTOR101 = beamManager.getData(par, groupedTriggers[ii], ii, "E:TOR101");
+        potTR101D = beamManager.getData(par, groupedTriggers[ii], ii, "E:TR101D");
+        hornI = beamManager.getData(par, groupedTriggers[ii], ii, "E:NSLIN");
+        hornDir = beamManager.getData(par, groupedTriggers[ii], ii, "E:HRNDIR");
+        horizontalposTGT = beamManager.getData(par, groupedTriggers[ii], ii, "E:HPTGT[]");
+        horizontalintTGT = beamManager.getData(par, groupedTriggers[ii], ii, "E:HITGT[]");
+        horizontalpos121 = beamManager.getData(par, groupedTriggers[ii], ii, "E:HP121[]");
+        verticalposTGT = beamManager.getData(par, groupedTriggers[ii], ii, "E:VPTGT[]");
+        verticalintTGT = beamManager.getData(par, groupedTriggers[ii], ii, "E:VITGT[]");
+        verticalpos121 = beamManager.getData(par, groupedTriggers[ii], ii, "E:VP121[]");
+        multiwireInfo = beamManager.getData(par, groupedTriggers[ii], ii, "E:MTGTDS[]");
     }
     else
     {
-	pot = par().runInfo().POTPerSpill() * 1e13;
+	potTRTGTD = {par().runInfo().POTPerSpill() * 1e13};
         caf.sr.beam.ismc = true;
     }
     if (std::isnan(caf.pot))
       caf.pot = 0;
-    caf.pot += pot;
-    caf.sr.beam.pulsepot = pot;
+    caf.pot += (potTRTGTD.size()==0) ? 0. : potTRTGTD.at(0);
+    caf.sr.beam.pulsepot = (potTRTGTD.size()==0) ? 0. : potTRTGTD.at(0);
+    caf.sr.beam.potTOR101 = (potTOR101.size()==0) ? 0. : potTOR101.at(0);
+    caf.sr.beam.potTR101D = (potTR101D.size()==0) ? 0. : potTR101D.at(0);
+    caf.sr.beam.hornI = (hornI.size()==0) ? 0. : hornI.at(0);
+    caf.sr.beam.hornDir = (hornDir.size()==0) ? 0. : hornDir.at(0);
+    caf.sr.beam.horizontalposTGT = horizontalposTGT;
+    caf.sr.beam.horizontalintTGT = horizontalintTGT;
+    caf.sr.beam.horizontalpos121 = horizontalpos121;
+    caf.sr.beam.verticalposTGT = verticalposTGT;
+    caf.sr.beam.verticalintTGT = verticalintTGT;
+    caf.sr.beam.verticalpos121 = verticalpos121;
+    caf.sr.beam.multiwireInfo = multiwireInfo;
+
+    cafmaker::LOG_S("").INFO() << "Filling caf tree \n";
     caf.fill();
   }
   progBar.Done();

src/reco/PandoraLArRecoNDBranchFiller.cxx

@@ -533,7 +533,7 @@ namespace cafmaker
         // unix_ts trigger time (seconds)
         trig.triggerTime_s = m_unixTime;
         // unix_time_usec ticks (microseconds) converted to nanoseconds
-        trig.triggerTime_ns = m_unixTimeUsec * 1000;
+        trig.triggerTime_ns = m_startTime * 100;
 
         LOG.VERBOSE() << "  added trigger: evtID = " << trig.evtID
                       << ", triggerType = " << trig.triggerType