HcalDigiPipelineTest.cxx

#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers.hpp>
 
using Catch::Approx;
 
#include "DetDescr/HcalID.h"  //creating unique hcal IDs
#include "Framework/Configure/Python.h"
#include "Framework/EventProcessor.h"
#include "Framework/Process.h"
#include "Hcal/Event/HcalHit.h"
#include "Recon/Event/HgcrocDigiCollection.h"
#include "SimCore/Event/SimCalorimeterHit.h"
 
namespace hcal {
namespace test {
 
static const double PE_ENERGY = 4.66 / 68;  // 0.069 MeV
 
// static const double MAX_ENERGY_ERROR_DAQ = 4 * PE_ENERGY;
// static const double MAX_ENERGY_PERCENT_ERROR_DAQ = 0.2;
static const double MAX_PE_ERROR_DAQ = 40;
// large percentage error for now
static const double MAX_PE_PERCENT_ERROR_DAQ = 0.4;
 
static const int NUM_TEST_SIM_HITS = 1000;
 
class IsCloseEnough : public Catch::Matchers::MatcherBase<double> {
 private:
  double truth_;
 
  const double MAX_ABSOLUTE_DIFF;
 
  const double MAX_RELATIVE_DIFF;
 
 public:
  IsCloseEnough(double const &truth, double const &abs_diff,
                double const &rel_diff)
      : truth_{truth},
        MAX_ABSOLUTE_DIFF{abs_diff},
        MAX_RELATIVE_DIFF{rel_diff} {}
 
  bool match(const double &daq) const override {
    return (daq == Approx(truth_).epsilon(MAX_RELATIVE_DIFF) or
            daq == Approx(truth_).margin(MAX_ABSOLUTE_DIFF));
  }
 
  virtual std::string describe() const override {
    std::ostringstream ss;
    ss << "is within an absolute difference of " << MAX_ABSOLUTE_DIFF
       << " OR a relative difference of " << MAX_RELATIVE_DIFF << " with "
       << truth_;
    return ss.str();
  }
};
 
class HcalFakeSimHits : public framework::Producer {
  // Based on the current gain settings for the ADC readout mode
  // we will reach saturation ~ 20 MeV ~ 290 PEs
  const double MAX_ENERGY = 200 * PE_ENERGY;  // ~ 13 MeV
 
  const double MIN_ENERGY = 4 * PE_ENERGY;
  const double ENERGY_STEP = (MAX_ENERGY - MIN_ENERGY) / NUM_TEST_SIM_HITS;
 
  double curr_energy_ = MIN_ENERGY;
 
 public:
  HcalFakeSimHits(const std::string &name, framework::Process &p)
      : framework::Producer(name, p) {}
  ~HcalFakeSimHits() {}
 
  void beforeNewRun(ldmx::RunHeader &header) final override {
    header.setDetectorName("ldmx-det-v12");
  }
 
  void produce(framework::Event &event) final override {
    // put in a single sim hit
    std::vector<ldmx::SimCalorimeterHit> pretend_sim_hits(1);
 
    // We hard-code the position of one hit: back hcal, layer_ 1, strip 31
    // This real simHit position is obtained by looking at calorimeter
    // SimHits of a 4 GeV muon shoot through the beamline
    ldmx::HcalID id(0, 1, 31);
    pretend_sim_hits[0].setPosition(-6.70265, 3.70265, 879);  // mm
    pretend_sim_hits[0].setID(id.raw());
    pretend_sim_hits[0].addContrib(
        -1,            // incidentID
        -1,            // trackID
        0,             // pdg ID
        curr_energy_,  // edep
        2.96628  // time - 299mm is about 1ns from target and in middle of HCal
    );
 
    // needs to be correct collection name
    // REQUIRE_NOTHROW(event.add("HcalSimHits", pretendSimHits));
    REQUIRE_NOTHROW(event.add("HcalFakeSimHits", pretend_sim_hits));
    curr_energy_ += ENERGY_STEP;
 
    return;
  }
};  // HcalFakeSimHits
 
class HcalCheckReconstruction : public framework::Analyzer {
  // save ntuple? False by default because if ntuplizer is on, the HcalGeometry
  // test cannot be run
  const bool SAVE = false;
 
 private:
  std::string hcal_fake_sim_hits_passname_;
  std::string hcal_rec_hits_passname_;
  std::string hcal_digis_passname_;
 
 public:
  HcalCheckReconstruction(const std::string &name, framework::Process &p)
      : framework::Analyzer(name, p) {}
  ~HcalCheckReconstruction() {}
 
  void configure(framework::config::Parameters &ps) override {
    hcal_fake_sim_hits_passname_ =
        ps.getParameter("hcal_fake_sim_hits_passname", "");
    hcal_digis_passname_ = ps.getParameter("hcal_digis_passname", "");
    hcal_rec_hits_passname_ = ps.getParameter("hcal_rec_hits_passname", "");
  }
 
  void onProcessStart() final override {
    if (SAVE) {
      getHistoDirectory();
      ntuple_.create("HcalDigiTest");
      ntuple_.addVar<float>("HcalDigiTest", "SimEnergy");
      ntuple_.addVar<float>("HcalDigiTest", "RecEnergy");
      ntuple_.addVar<float>("HcalDigiTest", "SimX");
      ntuple_.addVar<float>("HcalDigiTest", "SimY");
      ntuple_.addVar<float>("HcalDigiTest", "SimZ");
      ntuple_.addVar<float>("HcalDigiTest", "SimTime");
      ntuple_.addVar<float>("HcalDigiTest", "RecX");
      ntuple_.addVar<float>("HcalDigiTest", "RecY");
      ntuple_.addVar<float>("HcalDigiTest", "RecZ");
      ntuple_.addVar<float>("HcalDigiTest", "RecTime");
      ntuple_.addVar<int>("HcalDigiTest", "DaqDigi");
      ntuple_.addVar<int>("HcalDigiTest", "DaqDigiIsADC");
      ntuple_.addVar<int>("HcalDigiTest", "DaqDigiADC");
      ntuple_.addVar<int>("HcalDigiTest", "DaqDigiTOT");
    }
  }
 
  void analyze(const framework::Event &event) final override {
    const auto sim_hits = event.getCollection<ldmx::SimCalorimeterHit>(
        "HcalFakeSimHits", hcal_fake_sim_hits_passname_);
 
    REQUIRE(sim_hits.size() == 1);
 
    float truth_energy = sim_hits.at(0).getEdep();
 
    if (SAVE) {
      ntuple_.setVar<float>("SimEnergy", truth_energy);
      ntuple_.setVar<float>("SimX", sim_hits.at(0).getPosition()[0]);
      ntuple_.setVar<float>("SimY", sim_hits.at(0).getPosition()[1]);
      ntuple_.setVar<float>("SimZ", sim_hits.at(0).getPosition()[2]);
      ntuple_.setVar<float>("SimTime", sim_hits.at(0).getContrib(0).time_);
    }
 
    const auto daq_digis{event.getObject<ldmx::HgcrocDigiCollection>(
        "HcalDigis", hcal_digis_passname_)};
    auto daq_digi = daq_digis.getDigi(0);
    bool is_in_adc_mode = daq_digi.isADC();
 
    if (SAVE) {
      ntuple_.setVar<int>("DaqDigi", daq_digi.soi().raw());
      ntuple_.setVar<int>("DaqDigiIsADC", is_in_adc_mode);
      ntuple_.setVar<int>("DaqDigiADC", daq_digi.soi().adcT());
      ntuple_.setVar<int>("DaqDigiTOT", daq_digi.tot());
    }
 
    const auto rec_hits = event.getCollection<ldmx::HcalHit>(
        "HcalRecHits", hcal_rec_hits_passname_);
    CHECK(rec_hits.size() == 1);
 
    auto hit = rec_hits.at(0);
    ldmx::HcalID id(hit.getID());
    CHECK_FALSE(hit.isNoise());
    CHECK(id.raw() == sim_hits.at(0).getID());
 
    if (SAVE) {
      ntuple_.setVar<float>("RecX", hit.getXPos());
      ntuple_.setVar<float>("RecY", hit.getYPos());
      ntuple_.setVar<float>("RecZ", hit.getZPos());
      ntuple_.setVar<float>("RecTime", hit.getTime());
      ntuple_.setVar<float>("RecPE", hit.getPE());
      ntuple_.setVar<float>("RecEnergy", hit.getEnergy());
    }
 
    // define target pe by using the settings at the top
    double daq_pe{hit.getPE()};
    CHECK_THAT(daq_pe, IsCloseEnough(truth_energy / PE_ENERGY, MAX_PE_ERROR_DAQ,
                                     MAX_PE_PERCENT_ERROR_DAQ));
 
    // std::cout << "rec energy " << hit.getEnergy() << " * approx sampl
    // fraction " << hit.getEnergy()*sampling_fraction << " truth " <<
    //   truth_energy
    //           << std::endl;
    // std::cout << "npes " << hit.getPE() << " approx PE " << int(truth_energy
    // / PE_ENERGY)  << std::endl;
    /*
        if (id.section() == 0) {
          double truth_pos, rec_pos;
          if ((id.layer() % 2) == 1) {
            truth_pos = simHits.at(0).getPosition()[0];
            rec_pos = hit.getXPos();
          } else {
            truth_pos = simHits.at(0).getPosition()[1];
            rec_pos = hit.getYPos();
          }
          // std::cout << "rec pos_ " << rec_pos << " truth " << truth_pos <<
          // std::endl;
          // comment position check for now
          // CHECK_THAT(rec_pos, isCloseEnough(truth_pos,
       MAX_POSITION_ERROR_DAQ,
          //                                 MAX_POSITION_PERCENT_ERROR_DAQ));
        }
    */
    return;
  }
};  // HcalCheckReconstruction
 
}  // namespace test
}  // namespace hcal
 
DECLARE_ANALYZER(hcal::test::HcalFakeSimHits);
DECLARE_PRODUCER(hcal::test::HcalCheckReconstruction);
 
TEST_CASE("Hcal Digi Pipeline test", "[Hcal][functionality]") {
  const std::string config_file{"hcal_digi_pipeline_test_config.py"};
  char **args{nullptr};
 
  auto cfg{framework::config::run("ldmxcfg.Process.lastProcess", config_file,
                                  args, 0)};
  auto p{std::make_unique<framework::Process>(cfg)};
  p->run();
}