Skip to content

Commit

Permalink
Merge pull request #89 from DUNE/pot_info
Browse files Browse the repository at this point in the history
Interface for IFbeam
  • Loading branch information
sindhu-ku authored Oct 28, 2024
2 parents 78fce86 + 2cebf69 commit 9d00ee2
Show file tree
Hide file tree
Showing 8 changed files with 294 additions and 111 deletions.
5 changes: 4 additions & 1 deletion ndcaf_setup.sh
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ source /cvmfs/dune.opensciencegrid.org/products/dune/setup_dune.sh
setup cmake v3_22_2
setup gcc v9_3_0
setup pycurl
setup curl
setup ifdhc
setup geant4 v4_11_0_p01c -q e20:debug
setup dk2nugenie v01_10_01k -q debug:e20
Expand All @@ -14,7 +15,9 @@ setup hdf5 v1_10_5a -q e20
setup fhiclcpp v4_15_03 -q debug:e20
setup edepsim v3_2_0c -q debug:e20

# edep-sim needs to know where a certain GEANT .cmake file is...
export LD_LIBRARY_PATH=$CURL_ROOT/lib:$LD_LIBRARY_PATH

#edep-sim needs to know where a certain GEANT .cmake file is...
G4_cmake_file=`find ${GEANT4_FQ_DIR}/lib64 -name 'Geant4Config.cmake'`
export Geant4_DIR=`dirname $G4_cmake_file`

Expand Down
3 changes: 2 additions & 1 deletion src/Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@ SOURCES = $(shell find $(SRCDIR) -name '*.cxx')
HEADERS = $(shell find $(SRCDIR) -name '*.h')
OBJS = $(patsubst %.cxx, %.o, $(SOURCES))

INCLUDE += -I$(SRCDIR) -I$(SRPROXY_INC)
INCLUDE += -I$(SRCDIR) -I$(SRPROXY_INC) -I/usr/include/curl
LDLIBS += -L$(SANDRECO_LIB) -lStruct -L/usr/lib64 -lcurl

# enable TMS unless otherwise specified
# (set ENABLE_TMS=0 in the make command line)
Expand Down
6 changes: 3 additions & 3 deletions src/Params.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ namespace cafmaker

// fixme: placeholder. won't work for data, which will need to either pass through from upstream or interface to POT database here
fhicl::Atom<float> POTPerSpill { fhicl::Name("POTPerSpill"), fhicl::Comment("Fixed POT per spill (units of 10^13).")};

fhicl::Atom<bool> fhc { fhicl::Name("IsFHC"), fhicl::Comment("Is this an FHC run?"), true };
fhicl::Atom<bool> IsGasTPC { fhicl::Name("IsGasTPC"), fhicl::Comment("Was GArTPC geometry used? (If not, TMS)"), false };

Expand All @@ -45,12 +45,12 @@ 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 };

fhicl::Atom<bool> ForceDisableIFBeam { fhicl::Name("ForceDisableIFBeam"), fhicl::Comment("Forcefully disable IFBeam interface"), false}; //Disable IFBeam interface when needed (use case: running simulation without GENIE/edepsim)

// these are optional and have defaults
fhicl::Atom<int> first { fhicl::Name("FirstEvt"), fhicl::Comment("Start processing from this event number"), 0 };
fhicl::Atom<int> numevts { fhicl::Name("NumEvts"), fhicl::Comment("Number of events to process (-1 means 'all')"), -1 };
Expand Down
158 changes: 158 additions & 0 deletions src/beam/IFBeam.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,158 @@
/// \file IFBeam.cxx
///
/// IFBeam class to query beam information
///
/// \author S. Kumaran <[email protected]>
/// \date Oct. 2024

#include "IFBeam.h"
#include "util/Logger.h"
#include "util/IFBeamUtils.h"
#include <curl/curl.h>
#include <regex>
#include <sstream>
#include <limits>
#include <nlohmann/json.hpp>

using json = nlohmann::json;

namespace cafmaker
{


IFBeam::IFBeam(const std::vector<TriggerGroup>& groupedTriggers, bool is_data) {
if (is_data)
loadBeamSpills(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::ostringstream url_stream;
url_stream << "https://dbdata3vm.fnal.gov:9443/ifbeam/data/data?v=" << potDevice
<< "&e=" << "e,a9"
<< "&t0=" << min_time_iso
<< "&t1=" << max_time_iso
<< "&f=json";
return url_stream.str();
}

double IFBeam::unitToFactor(const std::string& unit) { //because the exponent part in IFBeam is stored as a string under "value" :/
std::regex re("E(\\d+)");
std::smatch match;

if (std::regex_search(unit, match, re) && match.size() > 1) {
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
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) {
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);
}
}
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);


CURL* curl;
CURLcode res;
std::string readBuffer;

std::cout << "Fetching beam information from: " << url << "\n";
curl = curl_easy_init();
if (curl) {
curl_easy_setopt(curl, CURLOPT_URL, url.c_str());
curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, util::WriteCallback);
curl_easy_setopt(curl, CURLOPT_WRITEDATA, &readBuffer);
res = curl_easy_perform(curl);
if (res != CURLE_OK) {
cafmaker::LOG_S("Fetch beam information").ERROR() << "curl_easy_perform() failed: " << curl_easy_strerror(res) << "\n";
std::abort();
}
curl_easy_cleanup(curl);
try {
auto json_data = json::parse(readBuffer);
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;
}
} catch (const json::exception& e) {
cafmaker::LOG_S("Fetch beam information").ERROR() << "Failed to parse JSON data: " << e.what() << "\n";
std::abort();
}
}
}

double IFBeam::getPOT(const cafmaker::Params& par, const TriggerGroup& groupedTrigger, int ii) {
double pot = 0.0;

auto it = std::find_if(beamSpills.begin(), beamSpills.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;
}
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) {
bool matched = std::any_of(beamSpills.cbegin(), beamSpills.cend(),
[par, &trig](const auto& spill) {
return std::abs(util::getTriggerTime(trig.second) - spill.first) < par().cafmaker().beamMatchDT();
});

if (matched) {
any_matched = true;
matched_triggers.push_back(trig);
} else {
unmatched_triggers.push_back(trig);
}
}

std::stringstream log_message;

if (any_matched) {
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() << " " << util::getTriggerTime(trig.second) << "\n";

log_message << "Unmatched triggers: \n";
for (auto trig : unmatched_triggers)
log_message << std::fixed << trig.first->GetName() << " " << util::getTriggerTime(trig.second) << "\n";

cafmaker::LOG_S("Fill beam POT information").ERROR() << log_message.str() << "\n";
std::abort();
}
else {
log_message << "No matching spill found for trigger group " << ii << " with triggers: \n";
for (auto trig : groupedTrigger)
log_message << std::fixed << trig.first->GetName() << " " << util::getTriggerTime(trig.second) << "\n";
cafmaker::LOG_S("Fill beam POT information").WARNING() << log_message.str() << "\n";
}
}

return pot;
}
}
40 changes: 40 additions & 0 deletions src/beam/IFBeam.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
/// \file IFBeam.h
///
/// IFBeam class to query beam information
///
/// \author S. Kumaran <[email protected]>
/// \date Oct. 2024

#ifndef IFBEAM_H
#define IFBEAM_H

#include <Params.h>
#include "reco/IRecoBranchFiller.h"
#include <string>
#include <map>
#include <vector>


namespace cafmaker
{
class IFBeam {
public:
using BeamSpills = std::map<double, double>;
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);


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);
double unitToFactor(const std::string& unit);
};

}
#endif // IFBEAM_H
Loading

0 comments on commit 9d00ee2

Please sign in to comment.