EcalDigiPipelineTest.cxx

#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers.hpp>
 
#include "DetDescr/EcalID.h"  //creating unique cell IDs
#include "Ecal/Event/EcalHit.h"
#include "Framework/Configure/Python.h"
#include "Framework/EventProcessor.h"
#include "Framework/Process.h"
#include "Recon/Event/HgcrocDigiCollection.h"
#include "Recon/Event/HgcrocTrigDigi.h"
#include "SimCore/Event/SimCalorimeterHit.h"
 
using Catch::Approx;
 
namespace ecal {
namespace test {
 
static const double MIP_SI_ENERGY = 0.130;
 
static const double MEV_PER_FC = MIP_SI_ENERGY / (37 * 0.1602);
 
static const double MAX_ENERGY_PERCENT_ERROR_DAQ = 0.025;
 
static const double MAX_ENERGY_PERCENT_ERROR_TP = 0.15;
 
static const double MAX_ENERGY_ERROR_DAQ = MIP_SI_ENERGY / 2;
 
static const double MAX_ENERGY_ERROR_TP = 2 * MIP_SI_ENERGY;
 
static const int NUM_TEST_SIM_HITS = 2000;
 
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_energy) const override {
    return (daq_energy == Approx(truth_).epsilon(MAX_RELATIVE_DIFF) or
            daq_energy == 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
       << "MeV OR a relative difference of " << MAX_RELATIVE_DIFF << " with "
       << truth_ << " MeV.";
    return ss.str();
  }
};
 
class EcalFakeSimHits : public framework::Producer {
  double max_energy_;
 
  double min_energy_;
 
  double max_time_;
 
  double min_time_;
 
  double energy_step_;
  double time_step_;
 
  double curr_energy_;
 
  double curr_time_;
 
 public:
  EcalFakeSimHits(const std::string& name, framework::Process& p)
      : framework::Producer(name, p) {}
  ~EcalFakeSimHits() {}
 
  void configure(framework::config::Parameters& ps) final override {
    min_energy_ = ps.get<double>("min_energy", MIP_SI_ENERGY);
    max_energy_ = ps.get<double>("max_energy", 10000. * MEV_PER_FC);
    // 299mm is about 1ns from target and in middle of ECal,
    // so the default arrival time of 1ns is an in-time hit
    min_time_ = ps.get<double>("min_time", 1.);
    max_time_ = ps.get<double>("max_time", 1.);
    // should match the number of events the process runs
    int n_steps = ps.get<int>("n_steps", NUM_TEST_SIM_HITS);
 
    energy_step_ = (max_energy_ - min_energy_) / n_steps;
    time_step_ = (max_time_ - min_time_) / n_steps;
 
    curr_energy_ = min_energy_;
    curr_time_ = min_time_;
  }
 
  void beforeNewRun(ldmx::RunHeader& header) final override {
    header.setDetectorName("ldmx-det-v15-8gev");
  }
 
  void produce(framework::Event& event) final override {
    // put in a single sim hit
    std::vector<ldmx::SimCalorimeterHit> pretend_sim_hits(1);
 
    ldmx::EcalID id(0, 0, 0);
    pretend_sim_hits[0].setID(id.raw());
    // incidentID, trackID, pdg ID, edep, time
    pretend_sim_hits[0].addContrib(-1, -1, 0, curr_energy_, curr_time_);
    // sim position in middle of ECal
    pretend_sim_hits[0].setPosition(0., 0., 299.);
 
    // needs to be correct collection name
    REQUIRE_NOTHROW(event.add("EcalSimHits", pretend_sim_hits));
 
    curr_energy_ += energy_step_;
    curr_time_ += time_step_;
 
    return;
  }
};  // EcalFakeSimHits
 
class EcalCheckEnergyReconstruction : public framework::Analyzer {
 private:
  std::string ecal_simhits_passname_;
  std::string ecal_digis_passname_;
  std::string ecal_rechits_passname_;
  std::string ecal_trig_digis_passname_;
  bool check_trig_prim_;
 
 public:
  EcalCheckEnergyReconstruction(const std::string& name, framework::Process& p)
      : framework::Analyzer(name, p) {}
  ~EcalCheckEnergyReconstruction() {}
 
  void configure(framework::config::Parameters& parameters) final override {
    // the trigger primitives are an in-time, single-sample estimate,
    // so they are not expected to see out-of-time hits
    check_trig_prim_ = parameters.get<bool>("check_trig_prim", true);
    ecal_simhits_passname_ =
        parameters.get<std::string>("ecal_simhits_passname", "");
    ecal_digis_passname_ =
        parameters.get<std::string>("ecal_digis_passname", "");
    ecal_rechits_passname_ =
        parameters.get<std::string>("ecal_rechits_passname", "");
    ecal_trig_digis_passname_ =
        parameters.get<std::string>("ecal_trig_digis_passname", "");
  }
 
  void onProcessStart() final override {
    getHistoDirectory();
    ntuple_.create("EcalDigiTest");
    ntuple_.addVar<float>("EcalDigiTest", "SimEnergy");
    ntuple_.addVar<float>("EcalDigiTest", "RecEnergy");
    ntuple_.addVar<float>("EcalDigiTest", "TrigPrimEnergy");
 
    ntuple_.addVar<int>("EcalDigiTest", "DaqDigi");
    ntuple_.addVar<int>("EcalDigiTest", "DaqDigiIsADC");
    ntuple_.addVar<int>("EcalDigiTest", "DaqDigiADC");
    ntuple_.addVar<int>("EcalDigiTest", "DaqDigiTOT");
    ntuple_.addVar<int>("EcalDigiTest", "TrigPrimDigiEncoded");
    ntuple_.addVar<int>("EcalDigiTest", "TrigPrimDigiLinear");
  }
 
  void analyze(const framework::Event& event) final override {
    const auto sim_hits = event.getCollection<ldmx::SimCalorimeterHit>(
        "EcalSimHits", ecal_simhits_passname_);
 
    REQUIRE(sim_hits.size() == 1);
 
    float truth_energy = sim_hits.at(0).getEdep();
    ntuple_.setVar<float>("SimEnergy", truth_energy);
 
    const auto daq_digis{event.getObject<ldmx::HgcrocDigiCollection>(
        "EcalDigis", ecal_digis_passname_)};
 
    if (daq_digis.getNumDigis() == 1) {
      auto daq_digi = daq_digis.getDigi(0);
      ntuple_.setVar<int>("DaqDigi", daq_digi.soi().raw());
      bool is_in_adc_mode = daq_digi.isADC();
      ntuple_.setVar<int>("DaqDigiIsADC", is_in_adc_mode);
      ntuple_.setVar<int>("DaqDigiADC", daq_digi.soi().adcT());
      ntuple_.setVar<int>("DaqDigiTOT", daq_digi.tot());
 
      // arrival time of the hit at the chip, helpful when the failure
      // depends on where in the readout window the pulse lands
      INFO("sim hit arrival time = " << sim_hits.at(0).getContrib(0).time_
                                     << " ns");
      INFO("digi is " << (is_in_adc_mode ? "ADC" : "TOT") << " mode");
 
      const auto rec_hits = event.getCollection<ldmx::EcalHit>(
          "EcalRecHits", ecal_rechits_passname_);
      CHECK(rec_hits.size() == 1);
      // a hit that was read out but not reconstructed has nothing left
      // for us to check, keep going so that all events are checked
      if (rec_hits.size() != 1) return;
 
      auto hit = rec_hits.at(0);
      ldmx::EcalID id(hit.getID());
      CHECK_FALSE(hit.isNoise());
      CHECK(id.raw() == sim_hits.at(0).getID());
 
      double daq_energy{hit.getAmplitude()};
      CHECK_THAT(daq_energy, IsCloseEnough(truth_energy, MAX_ENERGY_ERROR_DAQ,
                                           MAX_ENERGY_PERCENT_ERROR_DAQ));
      ntuple_.setVar<float>("RecEnergy", hit.getAmplitude());
 
      if (not check_trig_prim_) return;
 
      const auto trig_digis{event.getObject<ldmx::HgcrocTrigDigiCollection>(
          "ecalTrigDigis", ecal_trig_digis_passname_)};
      CHECK(trig_digis.size() == 1);
 
      auto trig_digi = trig_digis.at(0);
      float tp_energy =
          8 * trig_digi.linearPrimitive() * 320. / 1024 * MEV_PER_FC;
 
      CHECK_THAT(tp_energy, IsCloseEnough(truth_energy, MAX_ENERGY_ERROR_TP,
                                          MAX_ENERGY_PERCENT_ERROR_TP));
      ntuple_.setVar<float>("TrigPrimEnergy", tp_energy);
      ntuple_.setVar<int>("TrigPrimDigiEncoded", trig_digi.getPrimitive());
      ntuple_.setVar<int>("TrigPrimDigiLinear", trig_digi.linearPrimitive());
    }
 
    return;
  }
};  // EcalCheckEnergyReconstruction
 
}  // namespace test
}  // namespace ecal
 
DECLARE_PRODUCER(ecal::test::EcalFakeSimHits)
DECLARE_ANALYZER(ecal::test::EcalCheckEnergyReconstruction)
 
 
TEST_CASE("Ecal Digi Pipeline test", "[Ecal][functionality]") {
  const std::string config_file{"ecal_digi_pipeline_test_config.py"};
  char** args{nullptr};
 
  auto cfg{framework::config::run("ldmxcfg.Process.last_process", config_file,
                                  args, 0)};
  auto p{std::make_unique<framework::Process>(cfg)};
  p->run();
}
 
TEST_CASE("Ecal Digi Pipeline out-of-time test", "[Ecal][functionality]") {
  const std::string config_file{"ecal_digi_out_of_time_test_config.py"};
  char** args{nullptr};
 
  auto cfg{framework::config::run("ldmxcfg.Process.last_process", config_file,
                                  args, 0)};
  auto p{std::make_unique<framework::Process>(cfg)};
  p->run();
}