Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-29 11:58:11

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   if (conf.existsAs<edm::InputTag>("caloRecHits"))
0088     caloRecHitsToken_.emplace(consumes(conf.getUntrackedParameter<edm::InputTag>("caloRecHits")));
0089 
0090   const edm::InputTag input[2] = {conf.getUntrackedParameter<edm::InputTag>("pfRecHitsSource1"),
0091                                   conf.getUntrackedParameter<edm::InputTag>("pfRecHitsSource2")};
0092   const std::string type[2] = {conf.getUntrackedParameter<std::string>("pfRecHitsType1"),
0093                                conf.getUntrackedParameter<std::string>("pfRecHitsType2")};
0094   for (int i = 0; i < 2; i++) {
0095     if (type[i] == "legacy")
0096       pfRecHitsTokens_[i].emplace<LegacyToken>(consumes<LegacyHandle::element_type>(input[i]));
0097     else if (type[i] == "alpaka")
0098       pfRecHitsTokens_[i].emplace<AlpakaToken>(consumes<AlpakaHandle::element_type>(input[i]));
0099     else {
0100       fprintf(stderr, "Invalid value for PFRecHitProducerTest::pfRecHitsType%d: \"%s\"\n", i + 1, type[i].c_str());
0101       throw;
0102     }
0103   }
0104 }
0105 
0106 void PFRecHitProducerTest::endStream() {
0107   fprintf(stderr,
0108           "PFRecHitProducerTest%s%s%s has compared %u events and found %u problems: [%u, %u, %u, %u, %u]\n",
0109           title_.empty() ? "" : "[",
0110           title_.c_str(),
0111           title_.empty() ? "" : "]",
0112           num_events_,
0113           num_errors_,
0114           errors_[1],   // different number of PFRecHits
0115           errors_[2],   // detId not found
0116           errors_[3],   // depth,layer,time,energy or pos different
0117           errors_[4],   // different number of neighbours
0118           errors_[5]);  // neighbours different
0119 }
0120 
0121 void PFRecHitProducerTest::analyze(const edm::Event& event, const edm::EventSetup&) {
0122   GenericCollection pfRecHits[2];
0123   for (int i = 0; i < 2; i++)
0124     if (std::holds_alternative<LegacyToken>(pfRecHitsTokens_[i])) {
0125       pfRecHits[i].emplace<LegacyCollection>(&event.get(std::get<LegacyToken>(pfRecHitsTokens_[i])));
0126     } else {
0127       pfRecHits[i].emplace<AlpakaCollection>(&event.get(std::get<AlpakaToken>(pfRecHitsTokens_[i])).const_view());
0128     }
0129 
0130   int error = 0;
0131   const size_t n = GenericCollectionSize(pfRecHits[0]);
0132   if (n != GenericCollectionSize(pfRecHits[1]))
0133     error = 1;  // different number of PFRecHits
0134   else {
0135     std::vector<GenericPFRecHit> first, second;
0136     std::unordered_map<uint32_t, size_t> detId2Idx;  // for second vector
0137     first.reserve(n);
0138     second.reserve(n);
0139     for (size_t i = 0; i < n; i++) {
0140       first.emplace_back(GenericPFRecHit::Construct(pfRecHits[0], i));
0141       second.emplace_back(GenericPFRecHit::Construct(pfRecHits[1], i));
0142       detId2Idx[second.at(i).detId] = i;
0143     }
0144     for (size_t i = 0; i < n && error == 0; i++) {
0145       error = [&]() {
0146         const GenericPFRecHit& rh1 = first.at(i);
0147         if (detId2Idx.find(rh1.detId) == detId2Idx.end())
0148           return 2;  // detId not found
0149 
0150         const GenericPFRecHit& rh2 = second.at(detId2Idx.at(rh1.detId));
0151         assert(rh1.detId == rh2.detId);
0152         if (rh1.depth != rh2.depth || rh1.layer != rh2.layer || rh1.x != rh2.x || rh1.y != rh2.y || rh1.z != rh2.z)
0153           return 3;  // depth, layer or pos different
0154         if (strictCompare_ && (rh1.time != rh2.time || rh1.energy != rh2.energy))
0155           return 3;  // time or energy different
0156         hist_energy_->Fill(rh1.energy, rh2.energy);
0157         hist_time_->Fill(rh1.time, rh2.time);
0158 
0159         if (rh1.neighbours4.size() != rh2.neighbours4.size() || rh1.neighbours8.size() != rh2.neighbours8.size())
0160           return 4;  // different number of neighbours
0161 
0162         for (size_t i = 0; i < rh1.neighbours4.size(); i++)
0163           if (first.at(rh1.neighbours4[i]).detId != second.at(rh2.neighbours4[i]).detId)
0164             return 5;  // neighbours4 different
0165         for (size_t i = 0; i < rh1.neighbours8.size(); i++)
0166           if (first.at(rh1.neighbours8[i]).detId != second.at(rh2.neighbours8[i]).detId)
0167             return 5;  // neighbours8 different
0168 
0169         return 0;  // no error
0170       }();
0171     }
0172   }
0173 
0174   if (dumpFirstEvent_ && num_events_ == 0)
0175     dumpEvent(event, pfRecHits[0], pfRecHits[1]);
0176 
0177   if (error) {
0178     if (dumpFirstError_ && num_errors_ == 0) {
0179       printf("Error: %d\n", error);
0180       dumpEvent(event, pfRecHits[0], pfRecHits[1]);
0181     }
0182     num_errors_++;
0183     errors_[error]++;
0184   }
0185   num_events_++;
0186 }
0187 
0188 void PFRecHitProducerTest::dumpEvent(const edm::Event& event,
0189                                      const GenericCollection& pfRecHits1,
0190                                      const GenericCollection& pfRecHits2) {
0191   if (caloRecHitsToken_) {
0192     edm::Handle<reco::CaloRecHitHostCollection> caloRecHits;
0193     event.getByToken(*caloRecHitsToken_, caloRecHits);
0194     const reco::CaloRecHitHostCollection::ConstView view = caloRecHits->view();
0195     printf("Found %d recHits\n", view.metadata().size());
0196     for (int i = 0; i < view.metadata().size(); i++)
0197       printf("recHit %4d detId:%u energy:%f time:%f flags:%d\n",
0198              i,
0199              view.detId(i),
0200              view.energy(i),
0201              view.time(i),
0202              view.flags(i));
0203   }
0204 
0205   printf("Found %zd/%zd pfRecHits from first/second origin\n",
0206          GenericCollectionSize(pfRecHits1),
0207          GenericCollectionSize(pfRecHits2));
0208   for (size_t i = 0; i < GenericCollectionSize(pfRecHits1); i++)
0209     GenericPFRecHit::Construct(pfRecHits1, i).print("First", i);
0210   for (size_t i = 0; i < GenericCollectionSize(pfRecHits2); i++)
0211     GenericPFRecHit::Construct(pfRecHits2, i).print("Second", i);
0212 }
0213 
0214 void PFRecHitProducerTest::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0215   edm::ParameterSetDescription desc;
0216   desc.addOptionalUntracked<edm::InputTag>("caloRecHits")
0217       ->setComment("CaloRecHitSoA, if supplied, it is dumped alongside the PFRecHits");
0218   desc.addUntracked<edm::InputTag>("pfRecHitsSource1")->setComment("First PFRecHit list for comparison");
0219   desc.addUntracked<edm::InputTag>("pfRecHitsSource2")->setComment("Second PFRecHit list for comparison");
0220   desc.ifValue(edm::ParameterDescription<std::string>("pfRecHitsType1", "legacy", false),
0221                edm::allowedValues<std::string>("legacy", "alpaka"))
0222       ->setComment("Format of first PFRecHit list (legacy or alpaka)");
0223   desc.ifValue(edm::ParameterDescription<std::string>("pfRecHitsType2", "alpaka", false),
0224                edm::allowedValues<std::string>("legacy", "alpaka"))
0225       ->setComment("Format of second PFRecHit list (legacy or alpaka)");
0226   desc.addUntracked<std::string>("title", "")->setComment("Module name for printout");
0227   desc.addUntracked<bool>("dumpFirstEvent", false)
0228       ->setComment("Dump PFRecHits of first event, regardless of result of comparison");
0229   desc.addUntracked<bool>("dumpFirstError", false)->setComment("Dump PFRecHits upon first encountered error");
0230   desc.addUntracked<bool>("strictCompare", false)->setComment("Compare all floats for equality");
0231   descriptions.addDefault(desc);
0232 }
0233 
0234 void PFRecHitProducerTest::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0235   ibooker.setCurrentFolder("ParticleFlow/PFRecHitV");
0236   hist_energy_ = ibooker.book2D("energy", "energy;Input1;Input2;Entries", 100, 0, 100, 100, 0, 100);
0237   hist_time_ = ibooker.book2D("time", "time;Input1;Input2;Entries", 100, 0, 100, 100, 0, 100);
0238 }
0239 
0240 PFRecHitProducerTest::GenericPFRecHit::GenericPFRecHit(const reco::PFRecHit& pfRecHit)
0241     : detId(pfRecHit.detId()),
0242       depth(pfRecHit.depth()),
0243       layer(pfRecHit.layer()),
0244       time(pfRecHit.time()),
0245       energy(pfRecHit.energy()),
0246       x(pfRecHit.position().x()),
0247       y(pfRecHit.position().y()),
0248       z(pfRecHit.position().z()) {
0249   // Fill neighbours4 and neighbours8, then remove elements of neighbours4 from neighbours8
0250   // This is necessary, because there can be duplicates in the neighbour lists
0251   // This procedure correctly accounts for these multiplicities
0252   reco::PFRecHit::Neighbours pfRecHitNeighbours4 = pfRecHit.neighbours4();
0253   reco::PFRecHit::Neighbours pfRecHitNeighbours8 = pfRecHit.neighbours8();
0254   neighbours4.reserve(4);
0255   neighbours8.reserve(8);
0256   for (auto p = pfRecHitNeighbours8.begin(); p < pfRecHitNeighbours8.end(); p++)
0257     neighbours8.emplace_back(*p);
0258   for (auto p = pfRecHitNeighbours4.begin(); p < pfRecHitNeighbours4.end(); p++) {
0259     neighbours4.emplace_back(*p);
0260     auto idx = std::find(neighbours8.begin(), neighbours8.end(), *p);
0261     std::copy(idx + 1, neighbours8.end(), idx);
0262   }
0263   neighbours8.resize(pfRecHitNeighbours8.size() - pfRecHitNeighbours4.size());
0264 }
0265 
0266 PFRecHitProducerTest::GenericPFRecHit::GenericPFRecHit(
0267     const reco::PFRecHitHostCollection::ConstView::const_element pfRecHit)
0268     : detId(pfRecHit.detId()),
0269       depth(pfRecHit.depth()),
0270       layer(pfRecHit.layer()),
0271       time(pfRecHit.time()),
0272       energy(pfRecHit.energy()),
0273       x(pfRecHit.x()),
0274       y(pfRecHit.y()),
0275       z(pfRecHit.z()) {
0276   // Copy first four elements into neighbours4 and last four into neighbours8
0277   neighbours4.reserve(4);
0278   neighbours8.reserve(4);
0279   for (size_t k = 0; k < 4; k++)
0280     if (pfRecHit.neighbours()(k) != -1)
0281       neighbours4.emplace_back((uint32_t)pfRecHit.neighbours()(k));
0282   for (size_t k = 4; k < 8; k++)
0283     if (pfRecHit.neighbours()(k) != -1)
0284       neighbours8.emplace_back((uint32_t)pfRecHit.neighbours()(k));
0285 }
0286 
0287 void PFRecHitProducerTest::GenericPFRecHit::print(const char* prefix, size_t idx) {
0288   printf("%s %4lu detId:%u depth:%d layer:%d energy:%f time:%f pos:%f,%f,%f neighbours:%lu+%lu(",
0289          prefix,
0290          idx,
0291          detId,
0292          depth,
0293          layer,
0294          energy,
0295          time,
0296          x,
0297          y,
0298          z,
0299          neighbours4.size(),
0300          neighbours8.size());
0301   for (uint32_t j = 0; j < neighbours4.size(); j++)
0302     printf("%s%u", j ? "," : "", neighbours4[j]);
0303   printf(";");
0304   for (uint32_t j = 0; j < neighbours8.size(); j++)
0305     printf("%s%u", j ? "," : "", neighbours8[j]);
0306   printf(")\n");
0307 }
0308 
0309 #include "FWCore/Framework/interface/MakerMacros.h"
0310 DEFINE_FWK_MODULE(PFRecHitProducerTest);