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/ConfigurePython.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 {
  const double maxEnergy_ = 10000. * MeV_per_fC;
 
  const double minEnergy_ = MIP_SI_ENERGY;
 
  const double energyStep_ = (maxEnergy_ - minEnergy_) / NUM_TEST_SIM_HITS;
 
  double currEnergy_ = minEnergy_;
 
 public:
  EcalFakeSimHits(const std::string &name, framework::Process &p)
      : framework::Producer(name, p) {}
  ~EcalFakeSimHits() {}
 
  void beforeNewRun(ldmx::RunHeader &header) final override {
    header.setDetectorName("ldmx-det-v14-8gev");
  }
 
  void produce(framework::Event &event) final override {
    // put in a single sim hit
    std::vector<ldmx::SimCalorimeterHit> pretendSimHits(1);
 
    ldmx::EcalID id(0, 0, 0);
    pretendSimHits[0].setID(id.raw());
    // incidentID, trackID, pdg ID, edep, time - 299mm is about 1ns from target
    // and in middle of ECal
    pretendSimHits[0].addContrib(-1, -1, 0, currEnergy_, 1.);
    // sim position in middle of ECal
    pretendSimHits[0].setPosition(0., 0., 299.);
 
    // needs to be correct collection name
    REQUIRE_NOTHROW(event.add("EcalSimHits", pretendSimHits));
 
    currEnergy_ += energyStep_;
 
    return;
  }
};  // EcalFakeSimHits
 
class EcalCheckEnergyReconstruction : public framework::Analyzer {
 public:
  EcalCheckEnergyReconstruction(const std::string &name, framework::Process &p)
      : framework::Analyzer(name, p) {}
  ~EcalCheckEnergyReconstruction() {}
 
  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 simHits =
        event.getCollection<ldmx::SimCalorimeterHit>("EcalSimHits");
 
    REQUIRE(simHits.size() == 1);
 
    float truth_energy = simHits.at(0).getEdep();
    ntuple_.setVar<float>("SimEnergy", truth_energy);
 
    const auto daqDigis{
        event.getObject<ldmx::HgcrocDigiCollection>("EcalDigis")};
 
    CHECK(daqDigis.getNumDigis() == 1);
    auto daqDigi = daqDigis.getDigi(0);
    ntuple_.setVar<int>("DaqDigi", daqDigi.soi().raw());
    bool is_in_adc_mode = daqDigi.isADC();
    ntuple_.setVar<int>("DaqDigiIsADC", is_in_adc_mode);
    ntuple_.setVar<int>("DaqDigiADC", daqDigi.soi().adc_t());
    ntuple_.setVar<int>("DaqDigiTOT", daqDigi.tot());
 
    const auto recHits = event.getCollection<ldmx::EcalHit>("EcalRecHits");
    CHECK(recHits.size() == 1);
 
    auto hit = recHits.at(0);
    ldmx::EcalID id(hit.getID());
    CHECK_FALSE(hit.isNoise());
    CHECK(id.raw() == simHits.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());
 
    const auto trigDigis{
        event.getObject<ldmx::HgcrocTrigDigiCollection>("ecalTrigDigis")};
    CHECK(trigDigis.size() == 1);
 
    auto trigDigi = trigDigis.at(0);
    float tp_energy = 8 * trigDigi.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", trigDigi.getPrimitive());
    ntuple_.setVar<int>("TrigPrimDigiLinear", trigDigi.linearPrimitive());
 
    return;
  }
};  // EcalCheckEnergyReconstruction
 
}  // namespace test
}  // namespace ecal
 
DECLARE_PRODUCER_NS(ecal::test, EcalFakeSimHits)
DECLARE_ANALYZER_NS(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};
  framework::ProcessHandle p;
 
  framework::ConfigurePython cfg(config_file, args, 0);
  REQUIRE_NOTHROW(p = cfg.makeProcess());
  p->run();
}