PFRecHitProducer/test/PFRecHitProducerTest.cc

0001 #include <algorithm>
0002 #include <cstdio>
0003 #include <string>
0004 #include <optional>
0005 #include <map>
0006 #include <unordered_map>
0007 #include <utility>
0008 #include <variant>
0009
0010 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0011 #include "DQMServices/Core/interface/DQMStore.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "DataFormats/DetId/interface/DetId.h"
0014 #include "DataFormats/ParticleFlowReco/interface/CaloRecHitHostCollection.h"
0015 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0016 #include "DataFormats/ParticleFlowReco/interface/PFRecHitHostCollection.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/ParameterSet/interface/allowedValues.h"
0022
0023 class PFRecHitProducerTest : public DQMEDAnalyzer {
0024 public:
0025   PFRecHitProducerTest(edm::ParameterSet const& conf);
0026   void analyze(edm::Event const& e, edm::EventSetup const& c) override;
0027   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0028   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0029   void endStream() override;
0030
0031 private:
0032   // Define generic types for token, handle and collection, such that origin of PFRecHits can be selected at runtime.
0033   // The idea is to read desired format from config (pfRecHitsType1/2) and construct token accordingly. Then use handle
0034   // and collection corresponding to token type. Finally, construct GenericPFRecHits from each source.
0035   // This way, this module can be used to validate any combination of legacy and Alpaka formats.
0036
0037   using LegacyToken = edm::EDGetTokenT<reco::PFRecHitCollection>;
0038   using AlpakaToken = edm::EDGetTokenT<reco::PFRecHitHostCollection>;
0039   using GenericPFRecHitToken = std::variant<LegacyToken, AlpakaToken>;
0040   using LegacyHandle = edm::Handle<reco::PFRecHitCollection>;
0041   using AlpakaHandle = edm::Handle<reco::PFRecHitHostCollection>;
0042   using LegacyCollection = const reco::PFRecHitCollection*;
0043   using AlpakaCollection = const reco::PFRecHitHostCollection::ConstView*;
0044   using GenericCollection = std::variant<LegacyCollection, AlpakaCollection>;
0045
0046   static size_t GenericCollectionSize(const GenericCollection& collection) {
0047     return std::visit([](auto& c) -> size_t { return c->size(); }, collection);
0048   };
0049
0050   std::optional<edm::EDGetTokenT<reco::CaloRecHitHostCollection>> caloRecHitsToken_{};
0051   GenericPFRecHitToken pfRecHitsTokens_[2];
0052
0053   void dumpEvent(const edm::Event&, const GenericCollection&, const GenericCollection&);
0054   int32_t num_events_ = 0, num_errors_ = 0;
0055   std::map<int32_t, uint32_t> errors_;
0056   const std::string title_;
0057   const bool strictCompare_, dumpFirstEvent_, dumpFirstError_;
0058   MonitorElement *hist_energy_, *hist_time_;
0059
0060   // Container for PFRecHit, independent of how it was constructed
0061   struct GenericPFRecHit {
0062     uint32_t detId;
0063     int depth;
0064     PFLayer::Layer layer;
0065     float time;
0066     float energy;
0067     float x, y, z;
0068     std::vector<uint32_t> neighbours4;  // nearest neighbours
0069     std::vector<uint32_t> neighbours8;  // non-nearest neighbours
0070
0071     static GenericPFRecHit Construct(const GenericCollection& collection, size_t i) {  // Select appropriate constructor
0072       return std::visit([i](auto& c) { return GenericPFRecHit{(*c)[i]}; }, collection);
0073     };
0074     GenericPFRecHit(const reco::PFRecHit& pfRecHit);  // Constructor from legacy format
0075     GenericPFRecHit(
0076         const reco::PFRecHitHostCollection::ConstView::const_element pfRecHitsAlpaka);  // Constructor from Alpaka SoA
0077
0078     void print(const char* prefix, size_t idx);
0079   };
0080 };
0081
0082 PFRecHitProducerTest::PFRecHitProducerTest(const edm::ParameterSet& conf)
0083     : title_(conf.getUntrackedParameter<std::string>("title")),
0084       strictCompare_(conf.getUntrackedParameter<bool>("strictCompare")),
0085       dumpFirstEvent_(conf.getUntrackedParameter<bool>("dumpFirstEvent")),
0086       dumpFirstError_(conf.getUntrackedParameter<bool>("dumpFirstError")) {
0087   const auto& caloRecHits = conf.getUntrackedParameter<edm::InputTag>("caloRecHits", {});
0088   if (!caloRecHits.label().empty())
0089     caloRecHitsToken_.emplace(consumes(caloRecHits));
0090
0091   const edm::InputTag input[2] = {conf.getUntrackedParameter<edm::InputTag>("pfRecHitsSource1"),
0092                                   conf.getUntrackedParameter<edm::InputTag>("pfRecHitsSource2")};
0093   const std::string type[2] = {conf.getUntrackedParameter<std::string>("pfRecHitsType1"),
0094                                conf.getUntrackedParameter<std::string>("pfRecHitsType2")};
0095   for (int i = 0; i < 2; i++) {
0096     if (type[i] == "legacy")
0097       pfRecHitsTokens_[i].emplace<LegacyToken>(consumes<LegacyHandle::element_type>(input[i]));
0098     else if (type[i] == "alpaka")
0099       pfRecHitsTokens_[i].emplace<AlpakaToken>(consumes<AlpakaHandle::element_type>(input[i]));
0100     else {
0101       fprintf(stderr, "Invalid value for PFRecHitProducerTest::pfRecHitsType%d: \"%s\"\n", i + 1, type[i].c_str());
0102       throw;
0103     }
0104   }
0105 }
0106
0107 void PFRecHitProducerTest::endStream() {
0108   fprintf(stderr,
0109           "PFRecHitProducerTest%s%s%s has compared %u events and found %u problems: [%u, %u, %u, %u, %u]\n",
0110           title_.empty() ? "" : "[",
0111           title_.c_str(),
0112           title_.empty() ? "" : "]",
0113           num_events_,
0114           num_errors_,
0115           errors_[1],   // different number of PFRecHits
0116           errors_[2],   // detId not found
0117           errors_[3],   // depth,layer,time,energy or pos different
0118           errors_[4],   // different number of neighbours
0119           errors_[5]);  // neighbours different
0120 }
0121
0122 void PFRecHitProducerTest::analyze(const edm::Event& event, const edm::EventSetup&) {
0123   GenericCollection pfRecHits[2];
0124   for (int i = 0; i < 2; i++)
0125     if (std::holds_alternative<LegacyToken>(pfRecHitsTokens_[i])) {
0126       pfRecHits[i].emplace<LegacyCollection>(&event.get(std::get<LegacyToken>(pfRecHitsTokens_[i])));
0127     } else {
0128       pfRecHits[i].emplace<AlpakaCollection>(&event.get(std::get<AlpakaToken>(pfRecHitsTokens_[i])).const_view());
0129     }
0130
0131   int error = 0;
0132   const size_t n = GenericCollectionSize(pfRecHits[0]);
0133   if (n != GenericCollectionSize(pfRecHits[1]))
0134     error = 1;  // different number of PFRecHits
0135   else {
0136     std::vector<GenericPFRecHit> first, second;
0137     std::unordered_map<uint32_t, size_t> detId2Idx;  // for second vector
0138     first.reserve(n);
0139     second.reserve(n);
0140     for (size_t i = 0; i < n; i++) {
0141       first.emplace_back(GenericPFRecHit::Construct(pfRecHits[0], i));
0142       second.emplace_back(GenericPFRecHit::Construct(pfRecHits[1], i));
0143       detId2Idx[second.at(i).detId] = i;
0144     }
0145     for (size_t i = 0; i < n && error == 0; i++) {
0146       error = [&]() {
0147         const GenericPFRecHit& rh1 = first.at(i);
0148         if (detId2Idx.find(rh1.detId) == detId2Idx.end())
0149           return 2;  // detId not found
0150
0151         const GenericPFRecHit& rh2 = second.at(detId2Idx.at(rh1.detId));
0152         assert(rh1.detId == rh2.detId);
0153         if (rh1.depth != rh2.depth || rh1.layer != rh2.layer || rh1.x != rh2.x || rh1.y != rh2.y || rh1.z != rh2.z)
0154           return 3;  // depth, layer or pos different
0155         if (strictCompare_ && (rh1.time != rh2.time || rh1.energy != rh2.energy))
0156           return 3;  // time or energy different
0157         hist_energy_->Fill(rh1.energy, rh2.energy);
0158         hist_time_->Fill(rh1.time, rh2.time);
0159
0160         if (rh1.neighbours4.size() != rh2.neighbours4.size() || rh1.neighbours8.size() != rh2.neighbours8.size())
0161           return 4;  // different number of neighbours
0162
0163         for (size_t i = 0; i < rh1.neighbours4.size(); i++)
0164           if (first.at(rh1.neighbours4[i]).detId != second.at(rh2.neighbours4[i]).detId)
0165             return 5;  // neighbours4 different
0166         for (size_t i = 0; i < rh1.neighbours8.size(); i++)
0167           if (first.at(rh1.neighbours8[i]).detId != second.at(rh2.neighbours8[i]).detId)
0168             return 5;  // neighbours8 different
0169
0170         return 0;  // no error
0171       }();
0172     }
0173   }
0174
0175   if (dumpFirstEvent_ && num_events_ == 0)
0176     dumpEvent(event, pfRecHits[0], pfRecHits[1]);
0177
0178   if (error) {
0179     if (dumpFirstError_ && num_errors_ == 0) {
0180       printf("Error: %d\n", error);
0181       dumpEvent(event, pfRecHits[0], pfRecHits[1]);
0182     }
0183     num_errors_++;
0184     errors_[error]++;
0185   }
0186   num_events_++;
0187 }
0188
0189 void PFRecHitProducerTest::dumpEvent(const edm::Event& event,
0190                                      const GenericCollection& pfRecHits1,
0191                                      const GenericCollection& pfRecHits2) {
0192   if (caloRecHitsToken_) {
0193     edm::Handle<reco::CaloRecHitHostCollection> caloRecHits;
0194     event.getByToken(*caloRecHitsToken_, caloRecHits);
0195     const reco::CaloRecHitHostCollection::ConstView view = caloRecHits->view();
0196     printf("Found %d recHits\n", view.metadata().size());
0197     for (int i = 0; i < view.metadata().size(); i++)
0198       printf("recHit %4d detId:%u energy:%f time:%f flags:%d\n",
0199              i,
0200              view.detId(i),
0201              view.energy(i),
0202              view.time(i),
0203              view.flags(i));
0204   }
0205
0206   printf("Found %zd/%zd pfRecHits from first/second origin\n",
0207          GenericCollectionSize(pfRecHits1),
0208          GenericCollectionSize(pfRecHits2));
0209   for (size_t i = 0; i < GenericCollectionSize(pfRecHits1); i++)
0210     GenericPFRecHit::Construct(pfRecHits1, i).print("First", i);
0211   for (size_t i = 0; i < GenericCollectionSize(pfRecHits2); i++)
0212     GenericPFRecHit::Construct(pfRecHits2, i).print("Second", i);
0213 }
0214
0215 void PFRecHitProducerTest::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0216   edm::ParameterSetDescription desc;
0217   desc.addUntracked<edm::InputTag>("caloRecHits", {})
0218       ->setComment("CaloRecHitSoA, if supplied, it is dumped alongside the PFRecHits");
0219   desc.addUntracked<edm::InputTag>("pfRecHitsSource1")->setComment("First PFRecHit list for comparison");
0220   desc.addUntracked<edm::InputTag>("pfRecHitsSource2")->setComment("Second PFRecHit list for comparison");
0221   desc.ifValue(edm::ParameterDescription<std::string>("pfRecHitsType1", "legacy", false),
0222                edm::allowedValues<std::string>("legacy", "alpaka"))
0223       ->setComment("Format of first PFRecHit list (legacy or alpaka)");
0224   desc.ifValue(edm::ParameterDescription<std::string>("pfRecHitsType2", "alpaka", false),
0225                edm::allowedValues<std::string>("legacy", "alpaka"))
0226       ->setComment("Format of second PFRecHit list (legacy or alpaka)");
0227   desc.addUntracked<std::string>("title", "")->setComment("Module name for printout");
0228   desc.addUntracked<bool>("dumpFirstEvent", false)
0229       ->setComment("Dump PFRecHits of first event, regardless of result of comparison");
0230   desc.addUntracked<bool>("dumpFirstError", false)->setComment("Dump PFRecHits upon first encountered error");
0231   desc.addUntracked<bool>("strictCompare", false)->setComment("Compare all floats for equality");
0232   descriptions.addDefault(desc);
0233 }
0234
0235 void PFRecHitProducerTest::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0236   ibooker.setCurrentFolder("ParticleFlow/PFRecHitV");
0237   hist_energy_ = ibooker.book2D("energy", "energy;Input1;Input2;Entries", 100, 0, 100, 100, 0, 100);
0238   hist_time_ = ibooker.book2D("time", "time;Input1;Input2;Entries", 100, 0, 100, 100, 0, 100);
0239 }
0240
0241 PFRecHitProducerTest::GenericPFRecHit::GenericPFRecHit(const reco::PFRecHit& pfRecHit)
0242     : detId(pfRecHit.detId()),
0243       depth(pfRecHit.depth()),
0244       layer(pfRecHit.layer()),
0245       time(pfRecHit.time()),
0246       energy(pfRecHit.energy()),
0247       x(pfRecHit.position().x()),
0248       y(pfRecHit.position().y()),
0249       z(pfRecHit.position().z()) {
0250   // Fill neighbours4 and neighbours8, then remove elements of neighbours4 from neighbours8
0251   // This is necessary, because there can be duplicates in the neighbour lists
0252   // This procedure correctly accounts for these multiplicities
0253   reco::PFRecHit::Neighbours pfRecHitNeighbours4 = pfRecHit.neighbours4();
0254   reco::PFRecHit::Neighbours pfRecHitNeighbours8 = pfRecHit.neighbours8();
0255   neighbours4.reserve(4);
0256   neighbours8.reserve(8);
0257   for (auto p = pfRecHitNeighbours8.begin(); p < pfRecHitNeighbours8.end(); p++)
0258     neighbours8.emplace_back(*p);
0259   for (auto p = pfRecHitNeighbours4.begin(); p < pfRecHitNeighbours4.end(); p++) {
0260     neighbours4.emplace_back(*p);
0261     auto idx = std::find(neighbours8.begin(), neighbours8.end(), *p);
0262     std::copy(idx + 1, neighbours8.end(), idx);
0263   }
0264   neighbours8.resize(pfRecHitNeighbours8.size() - pfRecHitNeighbours4.size());
0265 }
0266
0267 PFRecHitProducerTest::GenericPFRecHit::GenericPFRecHit(
0268     const reco::PFRecHitHostCollection::ConstView::const_element pfRecHit)
0269     : detId(pfRecHit.detId()),
0270       depth(pfRecHit.depth()),
0271       layer(pfRecHit.layer()),
0272       time(pfRecHit.time()),
0273       energy(pfRecHit.energy()),
0274       x(pfRecHit.x()),
0275       y(pfRecHit.y()),
0276       z(pfRecHit.z()) {
0277   // Copy first four elements into neighbours4 and last four into neighbours8
0278   neighbours4.reserve(4);
0279   neighbours8.reserve(4);
0280   for (size_t k = 0; k < 4; k++)
0281     if (pfRecHit.neighbours()(k) != -1)
0282       neighbours4.emplace_back((uint32_t)pfRecHit.neighbours()(k));
0283   for (size_t k = 4; k < 8; k++)
0284     if (pfRecHit.neighbours()(k) != -1)
0285       neighbours8.emplace_back((uint32_t)pfRecHit.neighbours()(k));
0286 }
0287
0288 void PFRecHitProducerTest::GenericPFRecHit::print(const char* prefix, size_t idx) {
0289   printf("%s %4lu detId:%u depth:%d layer:%d energy:%f time:%f pos:%f,%f,%f neighbours:%lu+%lu(",
0290          prefix,
0291          idx,
0292          detId,
0293          depth,
0294          layer,
0295          energy,
0296          time,
0297          x,
0298          y,
0299          z,
0300          neighbours4.size(),
0301          neighbours8.size());
0302   for (uint32_t j = 0; j < neighbours4.size(); j++)
0303     printf("%s%u", j ? "," : "", neighbours4[j]);
0304   printf(";");
0305   for (uint32_t j = 0; j < neighbours8.size(); j++)
0306     printf("%s%u", j ? "," : "", neighbours8[j]);
0307   printf(")\n");
0308 }
0309
0310 #include "FWCore/Framework/interface/MakerMacros.h"
0311 DEFINE_FWK_MODULE(PFRecHitProducerTest);