CaloChecker_module.cc

Go to the documentation of this file.

#include <cmath>
#include <limits> // std::numeric_limits<>
#include <optional>
#include <string>

#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "larreco/Calorimetry/CalorimetryAlg.h"

#include "larcorealg/CoreUtils/NumericUtils.h" // util::absDiff()
#include "larcorealg/Geometry/PlaneGeo.h"
#include "larcorealg/Geometry/WireGeo.h"
#include "lardata/ArtDataHelper/TrackUtils.h" // lar::util::TrackPitchInView()
#include "lardata/Utilities/AssociationUtil.h"
#include "lardataobj/AnalysisBase/Calorimetry.h"
#include "lardataobj/AnalysisBase/T0.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/TrackHitMeta.h"
#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"

#include "larevt/SpaceCharge/SpaceCharge.h"
#include "larevt/SpaceChargeServices/SpaceChargeService.h"

// ROOT includes
#include <TF1.h>
#include <TGraph.h>
#include <TMath.h>

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "cetlib/pow.h" // cet::sum_of_squares()
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

namespace calo {

  class CaloChecker : public art::EDAnalyzer {
  public:
    explicit CaloChecker(fhicl::ParameterSet const& config);
    void analyze(const art::Event& evt) override;

  private:
    std::vector<std::string> fCaloLabels;
    std::string fTrackLabel;
  };

} // end namespace calo

calo::CaloChecker::CaloChecker(fhicl::ParameterSet const& config)
  : EDAnalyzer{config}
  , fCaloLabels(config.get<std::vector<std::string>>("CaloLabels"))
  , fTrackLabel(config.get<std::string>("TrackLabel"))
{
  assert(fCaloLabels.size() >= 2);
}

void calo::CaloChecker::analyze(const art::Event& evt)
{
  art::Handle<std::vector<recob::Track>> trackListHandle;
  std::vector<art::Ptr<recob::Track>> tracklist;
  if (evt.getByLabel(fTrackLabel, trackListHandle)) {
    art::fill_ptr_vector(tracklist, trackListHandle);
  }

  std::vector<art::FindManyP<anab::Calorimetry>> calos;
  for (const std::string& l : fCaloLabels) {
    calos.emplace_back(tracklist, evt, l);
  }

  float EPS = 1e-3;

  for (unsigned trk_i = 0; trk_i < tracklist.size(); trk_i++) {
    std::vector<art::Ptr<anab::Calorimetry>> base = calos[0].at(trk_i);

    std::cout << "New Track!\n";
    std::cout << "Base calo (" << fCaloLabels[0] << ") n calo: " << base.size() << std::endl;
    for (unsigned plane = 0; plane < base.size(); plane++) {
      std::cout << "N points on plane (" << plane << ") ID (" << base[plane]->PlaneID()
                << ") n points: " << base[plane]->dEdx().size() << std::endl;
    }

    for (unsigned i = 1; i < fCaloLabels.size(); i++) {
      bool equal = true;
      std::vector<art::Ptr<anab::Calorimetry>> othr = calos[i].at(trk_i);
      if (base.size() != othr.size()) {
        equal = false;
        std::cout << "Track: " << trk_i << " calo (" << fCaloLabels[0] << ") has (" << base.size()
                  << "). Calo (" << fCaloLabels[i] << ") has size (" << othr.size() << ") mismatch."
                  << std::endl;
      }

      for (unsigned plane = 0; plane < base.size(); plane++) {
        // check plane index
        if (base[plane]->PlaneID() != othr[plane]->PlaneID()) {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": ";
          std::cout << "Plane ID mismatch (" << base[plane]->PlaneID() << ") ("
                    << othr[plane]->PlaneID() << ")" << std::endl;
        }

        // check the range
        if (abs(base[plane]->Range() - othr[plane]->Range()) > EPS) {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": ";
          std::cout << "Range mismatch (" << base[plane]->Range() << ") (" << othr[plane]->Range()
                    << ")" << std::endl;
        }

        // check the kinetic energy
        if (abs(base[plane]->KineticEnergy() - othr[plane]->KineticEnergy()) > EPS) {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": ";
          std::cout << "KineticEnergy mismatch (" << base[plane]->KineticEnergy() << ") ("
                    << othr[plane]->KineticEnergy() << ")" << std::endl;
        }

        // check dE dx
        const std::vector<float>& base_dedx = base[plane]->dEdx();
        const std::vector<float>& othr_dedx = othr[plane]->dEdx();

        if (base_dedx.size() == othr_dedx.size()) {
          for (unsigned i_dedx = 0; i_dedx < base_dedx.size(); i_dedx++) {
            if (abs(base_dedx[i_dedx] - othr_dedx[i_dedx]) > EPS) {
              equal = false;
              std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") ("
                        << fCaloLabels[i] << ") plane " << plane << ": ";
              std::cout << "dEdx mismatch at index: " << i_dedx << " (" << base_dedx[i_dedx]
                        << ") (" << othr_dedx[i_dedx] << ")" << std::endl;
            }
          }
        }
        else {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": "
                    << "dEdx size mismatch (" << base_dedx.size() << ") (" << othr_dedx.size()
                    << ")" << std::endl;
        }

        // check dQdx
        const std::vector<float>& base_dqdx = base[plane]->dQdx();
        const std::vector<float>& othr_dqdx = othr[plane]->dQdx();

        if (base_dqdx.size() == othr_dqdx.size()) {
          for (unsigned i_dqdx = 0; i_dqdx < base_dqdx.size(); i_dqdx++) {
            if (abs(base_dqdx[i_dqdx] - othr_dqdx[i_dqdx]) > EPS) {
              equal = false;
              std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") ("
                        << fCaloLabels[i] << ") plane " << plane << ": ";
              std::cout << "dQdx mismatch at index: " << i_dqdx << " (" << base_dqdx[i_dqdx]
                        << ") (" << othr_dqdx[i_dqdx] << ")" << std::endl;
            }
          }
        }
        else {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": "
                    << "dQdx size mismatch (" << base_dqdx.size() << ") (" << othr_dqdx.size()
                    << ")" << std::endl;
        }

        // check Residual Range
        const std::vector<float>& base_rr = base[plane]->ResidualRange();
        const std::vector<float>& othr_rr = othr[plane]->ResidualRange();

        if (base_rr.size() == othr_rr.size()) {
          for (unsigned i_rr = 0; i_rr < base_rr.size(); i_rr++) {
            if (abs(base_rr[i_rr] - othr_rr[i_rr]) > EPS) {
              equal = false;
              std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") ("
                        << fCaloLabels[i] << ") plane " << plane << ": ";
              std::cout << "ResidualRange mismatch at index: " << i_rr << " (" << base_rr[i_rr]
                        << ") (" << othr_rr[i_rr] << ")" << std::endl;
            }
          }
        }
        else {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": "
                    << "ResidualRange size mismatch (" << base_rr.size() << ") (" << othr_rr.size()
                    << ")" << std::endl;
        }

        // check DeadWireRC
        const std::vector<float>& base_dwrr = base[plane]->DeadWireResRC();
        const std::vector<float>& othr_dwrr = othr[plane]->DeadWireResRC();

        if (base_dwrr.size() == othr_dwrr.size()) {
          for (unsigned i_dwrr = 0; i_dwrr < base_dwrr.size(); i_dwrr++) {
            if (abs(base_dwrr[i_dwrr] - othr_dwrr[i_dwrr]) > EPS) {
              equal = false;
              std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") ("
                        << fCaloLabels[i] << ") plane " << plane << ": ";
              std::cout << "DeadWireResRC mismatch at index: " << i_dwrr << " ("
                        << base_dwrr[i_dwrr] << ") (" << othr_dwrr[i_dwrr] << ")" << std::endl;
            }
          }
        }
        else {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": "
                    << "DeadWireResRC size mismatch (" << base_dwrr.size() << ") ("
                    << othr_dwrr.size() << ")" << std::endl;
        }

        // check Track Pitch
        const std::vector<float>& base_tpv = base[plane]->TrkPitchVec();
        const std::vector<float>& othr_tpv = othr[plane]->TrkPitchVec();

        if (base_tpv.size() == othr_tpv.size()) {
          for (unsigned i_tpv = 0; i_tpv < base_tpv.size(); i_tpv++) {
            if (abs(base_tpv[i_tpv] - othr_tpv[i_tpv]) > EPS) {
              equal = false;
              std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") ("
                        << fCaloLabels[i] << ") plane " << plane << ": ";
              std::cout << "TrkPitchVec mismatch at index: " << i_tpv << " (" << base_tpv[i_tpv]
                        << ") (" << othr_tpv[i_tpv] << ")" << std::endl;
            }
          }
        }
        else {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": "
                    << "TrkPitchVec size mismatch (" << base_tpv.size() << ") (" << othr_tpv.size()
                    << ")" << std::endl;
        }

        // check TP Indices
        const std::vector<size_t>& base_tpi = base[plane]->TpIndices();
        const std::vector<size_t>& othr_tpi = othr[plane]->TpIndices();

        if (base_tpi.size() == othr_tpi.size()) {
          for (unsigned i_tpi = 0; i_tpi < base_tpi.size(); i_tpi++) {
            if (base_tpi[i_tpi] != othr_tpi[i_tpi]) {
              equal = false;
              std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") ("
                        << fCaloLabels[i] << ") plane " << plane << ": ";
              std::cout << "TpIndices mismatch at index: " << i_tpi << " (" << base_tpi[i_tpi]
                        << ") (" << othr_tpi[i_tpi] << ")" << std::endl;
            }
          }
        }
        else {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": "
                    << "TpIndices size mismatch (" << base_tpi.size() << ") (" << othr_tpi.size()
                    << ")" << std::endl;
        }

        // check XYZ
        const std::vector<geo::Point_t>& base_xyz = base[plane]->XYZ();
        const std::vector<geo::Point_t>& othr_xyz = othr[plane]->XYZ();

        if (base_xyz.size() == othr_xyz.size()) {
          for (unsigned i_xyz = 0; i_xyz < base_xyz.size(); i_xyz++) {
            if (abs(base_xyz[i_xyz].X() - othr_xyz[i_xyz].X()) > EPS ||
                abs(base_xyz[i_xyz].Y() - othr_xyz[i_xyz].Y()) > EPS ||
                abs(base_xyz[i_xyz].Z() - othr_xyz[i_xyz].Z()) > EPS) {
              equal = false;
              std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") ("
                        << fCaloLabels[i] << ") plane " << plane << ": ";
              std::cout << "XYZ mismatch at index: " << i_xyz;
              std::cout << "X (" << base_xyz[i_xyz].X() << ") (" << othr_xyz[i_xyz].X() << ") ";
              std::cout << "Y (" << base_xyz[i_xyz].Y() << ") (" << othr_xyz[i_xyz].Y() << ") ";
              std::cout << "Z (" << base_xyz[i_xyz].Z() << ") (" << othr_xyz[i_xyz].Z() << ") "
                        << std::endl;
            }
          }
        }
        else {
          equal = false;
          std::cout << "Track: " << trk_i << " calos (" << fCaloLabels[0] << ") (" << fCaloLabels[i]
                    << ") plane " << plane << ": "
                    << "XYZ size mismatch (" << base_xyz.size() << ") (" << othr_xyz.size() << ")"
                    << std::endl;
        }
      }

      if (equal) {
        std::cout << "Track: " << trk_i << " calo (" << fCaloLabels[0] << ") is equal to calo ("
                  << fCaloLabels[i] << ")" << std::endl;
      }
    }
  }
}

DEFINE_ART_MODULE(calo::CaloChecker)