Merge branch 'main' of github.com:DUNE/ND_CAFMaker into feature/Pando…

…raLArRecoND
DUNE · Oct 3, 2024 · 89c4003 · 89c4003
2 parents 8f9264f + 3b93cfe
commit 89c4003
Show file tree

Hide file tree

Showing 20 changed files with 992 additions and 202 deletions.
diff --git a/Makefile b/Makefile
@@ -4,6 +4,7 @@ export ROOTFLAGS = $(shell root-config --cflags)
 export INCLUDE = -I$(HDF5_INC)
 INCLUDE += -I$(GENIE_INC)/GENIE
 INCLUDE += -I$(LOG4CPP_INC)
+INCLUDE += -I$(EDEPSIM_INC)/EDepSim
 #INCLUDE += -I$(NUSYST) -I$(NUSYST)/build/systematicstools/src/systematicstools
 #INCLUDE += -I$(NUSYST)/build/Linux/include/
 INCLUDE += -I$(BOOST_INC)
@@ -12,7 +13,6 @@ INCLUDE += -I$(CETLIB_EXCEPT_INC)
 INCLUDE += -I$(FHICLCPP_INC)
 INCLUDE += -I$(DUNEANAOBJ_INC)
 INCLUDE += -I$src
-
 export LDLIBS += -L$(LOG4CPP_LIB) -llog4cpp
 LDLIBS += -L$(TBB_LIB) -ltbb
 LDLIBS += -L$(LIBXML2_FQ_DIR)/lib -lxml2
@@ -32,7 +32,7 @@ LDLIBS += -lGeom -lEGPythia6 -lGenVector
 LDLIBS += -L$(BOOST_LIB) -lboost_program_options
 LDLIBS += -L$(CETLIB_LIB) -L$(CETLIB_EXCEPT_LIB) -lcetlib -lcetlib_except
 LDLIBS += -L$(FHICLCPP_LIB) -lfhiclcpp -lfhiclcpp_types
-
+LDLIBS += -L$(EDEPSIM_LIB) -ledepsim 
 export LIBDIR = $(PWD)/lib
 export BINDIR = $(PWD)/bin
 

diff --git a/build.sh b/build.sh
@@ -8,6 +8,10 @@ if [[ $# == 1 && x$1 == x-f ]]; then FORCE=yes; fi
 # set up software
 source ndcaf_setup.sh
 
+# 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`
+
 # Just use the edep-sim UPS product, don't clone master branch off repos!
 # Get edep-sim and build it
 #if [ $FORCE == yes ]; then rm -rf edep-sim; fi
@@ -57,7 +61,7 @@ cd ${TOPDIR}
 # Add pyGeoEff to pythonpath
 export PYTHONPATH=${PYTHONPATH}:${TOPDIR}/DUNE_ND_GeoEff/lib/
 
+
 # make tarballs of edep-sim and nusystematics for grid jobs
-#tar -zcf edep-sim.tar.gz edep-sim
 #tar -zcf nusystematics.tar.gz nusystematics
 tar -zcf DUNE_ND_GeoEff.tar.gz DUNE_ND_GeoEff
diff --git a/ndcaf_setup.sh b/ndcaf_setup.sh
@@ -1,19 +1,19 @@
 source /cvmfs/dune.opensciencegrid.org/products/dune/setup_dune.sh
 setup cmake v3_22_2
-setup gcc v12_1_0
+setup gcc v9_3_0
 setup pycurl
 setup ifdhc
-setup geant4 v4_10_6_p01g -q debug:e26
-setup dk2nugenie   v01_10_01q -q e26:prof
+setup geant4 v4_11_0_p01c -q e20:debug
+setup dk2nugenie   v01_10_01k -q debug:e20
 setup genie_xsec   v3_04_00 -q AR2320i00000:e1000:k250
 setup genie_phyopt v3_04_00 -q dkcharmtau
 setup jobsub_client
 setup eigen v3_3_5
-setup duneanaobj v03_05_00 -q debug:e26
+setup duneanaobj v03_06_01b -q e20:prof
 setup hdf5 v1_10_5a -q e20
-setup fhiclcpp v4_18_04 -q debug:e26
-setup root v6_28_12 -q e26:p3915:prof
-setup genie v3_04_02 -q e26:prof
+setup fhiclcpp v4_15_03 -q debug:e20
+setup edepsim v3_2_0c -q debug:e20
+setup root v6_26_06b -q e20:p3913:prof
 
 # edep-sim needs to know where a certain GEANT .cmake file is...
 G4_cmake_file=`find ${GEANT4_FQ_DIR}/lib64 -name 'Geant4Config.cmake'`
@@ -36,7 +36,6 @@ export PYTHONPATH=${PYTHONPATH}:${PWD}/DUNE_ND_GeoEff/lib/
 
 # duneananobj needs to be in the libs too
 export LD_LIBRARY_PATH=${DUNEANAOBJ_LIB}:$LD_LIBRARY_PATH
-
 # finally, add our lib & bin to the paths
 mydir="$( cd "$( dirname "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )"
 export LD_LIBRARY_PATH=$mydir/lib:$LD_LIBRARY_PATH

diff --git a/src/Params.h b/src/Params.h
@@ -32,7 +32,8 @@ namespace cafmaker
   {
     // these are mandatory and have no default values
     fhicl::OptionalSequence<std::string> GHEPFiles     { fhicl::Name{"GHEPFiles"},   fhicl::Comment("Input .ghep (GENIE) file(s) for truth matching") };
-    fhicl::Atom<std::string>             outputFile    { fhicl::Name{"OutputFile"},  fhicl::Comment("Filename for output CAF") };
+    fhicl::OptionalAtom<std::string> edepsimFile     { fhicl::Name{"EdepsimFile"},   fhicl::Comment("Input .root (EDPSIM) file for truth matching") };
+    fhicl::Atom<std::string> outputFile    { fhicl::Name{"OutputFile"},  fhicl::Comment("Filename for output CAF") };
 
     // this one is mandatory but has a default.  (the 'fhicl.fcl' file is provided in the 'sim_inputs' directory).
     // fixme: this file is currently not used for anything, but will be once DIRT-II is done and re-enables the interaction systematics
@@ -45,6 +46,9 @@ namespace cafmaker
     fhicl::OptionalAtom<std::string> minervaRecoFile  { fhicl::Name{"MINERVARecoFile"}, fhicl::Comment("Input MINERVA reco .root file") };
     fhicl::OptionalAtom<std::string> pandoraLArRecoNDFile  { fhicl::Name{"PandoraLArRecoNDFile"}, fhicl::Comment("Input Pandora LArRecoND .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 };
 
@@ -56,6 +60,17 @@ 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 };
+
+
+    // Track matching criteria (default values for 2x2 based on DOCDB 31970)
+    fhicl::Atom<float> trackMatchExtrapolatedZ { fhicl::Name("TrackMatchExtrapolatedZ"), fhicl::Comment("Z position where the track transversal displacement is computed when doing matching"), -70};
+    fhicl::Atom<float> trackMatchdX { fhicl::Name("TrackMatchDeltaX"), fhicl::Comment("Maximum displacement authorised in X [cm]"), 17};
+    fhicl::Atom<float> trackMatchdY { fhicl::Name("TrackMatchDeltaY"), fhicl::Comment("Maximum displacement authorised in Y [cm]"), 19};
+    fhicl::Atom<float> trackMatchdThetaX { fhicl::Name("TrackMatchDeltaThetaX"), fhicl::Comment("Maximum angle difference with respect to X axis [rad]"), .08};
+    fhicl::Atom<float> trackMatchdThetaY { fhicl::Name("TrackMatchDeltaThetaY"), fhicl::Comment("Maximum angle difference with respect to Y axis [rad]"), .09};
+
     // options are VERBOSE, DEBUG, INFO, WARNING, ERROR, FATAL
     fhicl::Atom<std::string> verbosity { fhicl::Name("Verbosity"), fhicl::Comment("Verbosity level of output"), "WARNING" };
   };

diff --git a/src/makeCAF.C b/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"
@@ -23,10 +29,12 @@
 #include "Params.h"
 #include "reco/MLNDLArRecoBranchFiller.h"
 #include "reco/TMSRecoBranchFiller.h"
-#include "reco/NDLArTMSMatchRecoFiller.h"
-#include "reco/SANDRecoBranchFiller.h"
 #include "reco/MINERvARecoBranchFiller.h"
+
+#include "reco/NDLArTMSMatchRecoFiller.h"
+#include "reco/NDLArMINERvAMatchRecoFiller.h"
 #include "reco/PandoraLArRecoNDBranchFiller.h"
+#include "reco/SANDRecoBranchFiller.h"
 #include "truth/FillTruth.h"
 #include "util/GENIEQuiet.h"
 #include "util/Logger.h"
@@ -135,6 +143,7 @@ std::vector<std::unique_ptr<cafmaker::IRecoBranchFiller>> getRecoFillers(const c
     recoFillers.emplace_back(std::make_unique<cafmaker::SANDRecoBranchFiller>(sandFile));
     std::cout << "   SAND\n";
   }
+
   // Pandora LArRecoND
   std::string pandoraFile;
   if (par().cafmaker().pandoraLArRecoNDFile(pandoraFile))
@@ -166,6 +175,11 @@ std::vector<std::unique_ptr<cafmaker::IRecoBranchFiller>> getRecoFillers(const c
     std::cout << "   ND-LAr + TMS matching\n";
   }
 
+  if (!ndlarFile.empty() && !minervaFile.empty())
+  {
+    recoFillers.emplace_back(std::make_unique<cafmaker::NDLArMINERvAMatchRecoFiller>(par().cafmaker().trackMatchExtrapolatedZ(), par().cafmaker().trackMatchdX(), par().cafmaker().trackMatchdY(), par().cafmaker().trackMatchdThetaX(), par().cafmaker().trackMatchdThetaY()));
+    std::cout << "   ND-LAr + MINERvA matching\n";
+  }
   // for now all the fillers get the same threshold.
   // if we decide we need to do it differently later
   // we can adjust the FCL params...
@@ -178,7 +192,8 @@ 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);
+  const unsigned long int s_to_ns = 1e9;
+  return ( (std::max(t1.triggerTime_s, t2.triggerTime_s) - std::min(t1.triggerTime_s, t2.triggerTime_s)) * s_to_ns + std::max(t1.triggerTime_ns, t2.triggerTime_ns) - std::min(t1.triggerTime_ns, t2.triggerTime_ns) ) < dT;
 }
 
 struct triggerTimeCmp
@@ -332,23 +347,125 @@ 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,
           cafmaker::Params &par,
           const std::vector<std::string> & ghepFilenames,
+          string edepsimFilename,
           const std::vector<std::unique_ptr<cafmaker::IRecoBranchFiller>> &recoFillers)
 {
   // if this is a data file, there won't be any truth, of course,
   // but the TruthMatching knows not to try to do anything with a null gtree
   cafmaker::Logger::THRESHOLD thresh = cafmaker::Logger::parseStringThresh(par().cafmaker().verbosity());
-  cafmaker::TruthMatcher truthMatcher(ghepFilenames, caf.mcrec,
+  cafmaker::TruthMatcher truthMatcher(ghepFilenames, edepsimFilename , caf.mcrec,
                                       [&caf](const genie::NtpMCEventRecord* mcrec){ return caf.StoreGENIEEvent(mcrec); });
   truthMatcher.SetLogThrehsold(thresh);
   // figure out which triggers we need to loop over between the various reco fillers
   std::map<const cafmaker::IRecoBranchFiller*, std::deque<cafmaker::Trigger>> triggersByRBF;
   for (const std::unique_ptr<cafmaker::IRecoBranchFiller>& filler : recoFillers)
+  {
+    if (filler->FillerType() != cafmaker::RecoFillerType::BaseReco) continue; //We don't want to store a trigger from a Matcher algorithm
     triggersByRBF.insert({filler.get(), filler->GetTriggers()});
+  }
   std::vector<std::vector<std::pair<const cafmaker::IRecoBranchFiller*, cafmaker::Trigger>>>
     groupedTriggers = buildTriggerList(triggersByRBF, par().cafmaker().trigMatchDT());
 
@@ -381,12 +498,6 @@ 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])
@@ -395,6 +506,22 @@ void loop(CAF &caf,
       fillerTrigPair.first->FillRecoBranches(fillerTrigPair.second, caf.sr, par, &truthMatcher);
     }
 
+    // Once all the reco fillers have been called, let's call the matching fillers
+    for (const std::unique_ptr<cafmaker::IRecoBranchFiller>& filler : recoFillers)
+    {
+      if (filler->FillerType() == cafmaker::RecoFillerType::Matcher)
+      {
+        filler->FillRecoBranches(groupedTriggers[ii][0].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();
   }
   progBar.Done();
@@ -418,11 +545,13 @@ int main( int argc, char const *argv[] )
   cafmaker::QuietGENIE();  // the GENIE events were already made earlier, we don't need more warnings about them
 
   std::vector<std::string> GHEPFiles;
+  std::string edepsimFile;
   par().cafmaker().GHEPFiles(GHEPFiles);  // fills the vector in if the key is found
+  par().cafmaker().edepsimFile(edepsimFile);  // fills the vector in if the key is found
 
   CAF caf(par().cafmaker().outputFile(), par().cafmaker().nusystsFcl(), par().cafmaker().makeFlatCAF(), !GHEPFiles.empty());
 
-  loop(caf, par, GHEPFiles, getRecoFillers(par, logThresh));
+  loop(caf, par, GHEPFiles, edepsimFile, getRecoFillers(par, logThresh));
 
   caf.version = 5;
   printf( "Run %d POT %g\n", caf.meta_run, caf.pot );