Local/plugins/TotemRPUVPatternFinder.cc

0001 /****************************************************************************
0002 *
0003 * This is a part of TOTEM offline software.
0004 * Authors:
0005 *   Jan Kašpar (jan.kaspar@gmail.com)
0006 *
0007 ****************************************************************************/
0008
0009 #include "FWCore/Framework/interface/stream/EDProducer.h"
0010 #include "FWCore/Framework/interface/ESWatcher.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/ESGetToken.h"
0017
0018 #include "DataFormats/Common/interface/DetSetVector.h"
0019 #include "DataFormats/Common/interface/DetSet.h"
0020 #include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h"
0021 #include "DataFormats/CTPPSReco/interface/TotemRPRecHit.h"
0022 #include "DataFormats/CTPPSReco/interface/TotemRPUVPattern.h"
0023
0024 #include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
0025 #include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
0026
0027 #include "RecoPPS/Local/interface/FastLineRecognition.h"
0028
0029 //----------------------------------------------------------------------------------------------------
0030
0031 /**
0032  * \brief Class to recognize straight line tracks, based on optimized Hough trasform.
0033  *
0034  * The search is perfomed in global U,V coordinates (wrt. beam). In this way (some of)
0035  * the alignment corrections can be taken into account.
0036 **/
0037 class TotemRPUVPatternFinder : public edm::stream::EDProducer<> {
0038 public:
0039   TotemRPUVPatternFinder(const edm::ParameterSet &conf);
0040
0041   ~TotemRPUVPatternFinder() override;
0042
0043   void produce(edm::Event &e, const edm::EventSetup &c) override;
0044   static void fillDescriptions(edm::ConfigurationDescriptions &);
0045
0046 private:
0047   edm::InputTag tagRecHit;
0048   edm::EDGetTokenT<edm::DetSetVector<TotemRPRecHit>> detSetVectorTotemRPRecHitToken;
0049
0050   edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> ctppsGeometryToken;
0051   edm::ESWatcher<VeryForwardRealGeometryRecord> geometryWatcher;
0052
0053   unsigned int verbosity;
0054
0055   /// minimal required number of active planes per projection to even start track recognition
0056   unsigned char minPlanesPerProjectionToSearch;
0057
0058   /// minimal required number of active planes per projection to mark track candidate as fittable
0059   unsigned char minPlanesPerProjectionToFit;
0060
0061   /// above this limit, planes are considered noisy
0062   unsigned int maxHitsPerPlaneToSearch;
0063
0064   /// the line recognition algorithm
0065   FastLineRecognition *lrcgn;
0066
0067   /// minimal weight of (Hough) cluster to accept it as candidate
0068   double threshold;
0069
0070   /// maximal angle (in any projection) to mark candidate as fittable - controls track parallelity
0071   double max_a_toFit;
0072
0073   /// block of (exceptional) settings for 1 RP
0074   struct RPSettings {
0075     unsigned char minPlanesPerProjectionToFit_U, minPlanesPerProjectionToFit_V;
0076     double threshold_U, threshold_V;
0077   };
0078
0079   /// exceptional settings: RP Id --> settings
0080   std::map<unsigned int, RPSettings> exceptionalSettings;
0081
0082   /// executes line recognition in a projection
0083   void recognizeAndSelect(TotemRPUVPattern::ProjectionType proj,
0084                           double z0,
0085                           double threshold,
0086                           unsigned int planes_required,
0087                           const edm::DetSetVector<TotemRPRecHit> &hits,
0088                           edm::DetSet<TotemRPUVPattern> &patterns);
0089 };
0090
0091 //----------------------------------------------------------------------------------------------------
0092
0093 using namespace std;
0094 using namespace edm;
0095
0096 //----------------------------------------------------------------------------------------------------
0097
0098 TotemRPUVPatternFinder::TotemRPUVPatternFinder(const edm::ParameterSet &conf)
0099     : tagRecHit(conf.getParameter<edm::InputTag>("tagRecHit")),
0100       ctppsGeometryToken(esConsumes()),
0101       verbosity(conf.getUntrackedParameter<unsigned int>("verbosity", 0)),
0102       minPlanesPerProjectionToSearch(conf.getParameter<unsigned int>("minPlanesPerProjectionToSearch")),
0103       minPlanesPerProjectionToFit(conf.getParameter<unsigned int>("minPlanesPerProjectionToFit")),
0104       maxHitsPerPlaneToSearch(conf.getParameter<unsigned int>("maxHitsPerPlaneToSearch")),
0105       lrcgn(new FastLineRecognition(conf.getParameter<double>("clusterSize_a"),
0106                                     conf.getParameter<double>("clusterSize_b"))),
0107       threshold(conf.getParameter<double>("threshold")),
0108       max_a_toFit(conf.getParameter<double>("max_a_toFit")) {
0109   for (const auto &ps : conf.getParameter<vector<ParameterSet>>("exceptionalSettings")) {
0110     unsigned int rpId = ps.getParameter<unsigned int>("rpId");
0111
0112     RPSettings settings;
0113     settings.minPlanesPerProjectionToFit_U = ps.getParameter<unsigned int>("minPlanesPerProjectionToFit_U");
0114     settings.minPlanesPerProjectionToFit_V = ps.getParameter<unsigned int>("minPlanesPerProjectionToFit_V");
0115     settings.threshold_U = ps.getParameter<double>("threshold_U");
0116     settings.threshold_V = ps.getParameter<double>("threshold_V");
0117
0118     exceptionalSettings[rpId] = settings;
0119   }
0120
0121   detSetVectorTotemRPRecHitToken = consumes<edm::DetSetVector<TotemRPRecHit>>(tagRecHit);
0122
0123   produces<DetSetVector<TotemRPUVPattern>>();
0124 }
0125
0126 //----------------------------------------------------------------------------------------------------
0127
0128 TotemRPUVPatternFinder::~TotemRPUVPatternFinder() { delete lrcgn; }
0129
0130 //----------------------------------------------------------------------------------------------------
0131
0132 void TotemRPUVPatternFinder::recognizeAndSelect(TotemRPUVPattern::ProjectionType proj,
0133                                                 double z0,
0134                                                 double threshold_loc,
0135                                                 unsigned int planes_required,
0136                                                 const DetSetVector<TotemRPRecHit> &hits,
0137                                                 DetSet<TotemRPUVPattern> &patterns) {
0138   // run recognition
0139   DetSet<TotemRPUVPattern> newPatterns;
0140   lrcgn->getPatterns(hits, z0, threshold_loc, newPatterns);
0141
0142   // set pattern properties and copy to the global pattern collection
0143   for (auto &p : newPatterns) {
0144     p.setProjection(proj);
0145
0146     p.setFittable(true);
0147
0148     set<unsigned int> planes;
0149     for (const auto &ds : p.hits())
0150       planes.insert(TotemRPDetId(ds.detId()).plane());
0151
0152     if (planes.size() < planes_required)
0153       p.setFittable(false);
0154
0155     if (fabs(p.a()) > max_a_toFit)
0156       p.setFittable(false);
0157
0158     patterns.push_back(p);
0159   }
0160 }
0161
0162 //----------------------------------------------------------------------------------------------------
0163
0164 void TotemRPUVPatternFinder::produce(edm::Event &event, const edm::EventSetup &es) {
0165   if (verbosity > 5)
0166     LogVerbatim("TotemRPUVPatternFinder")
0167         << ">> TotemRPUVPatternFinder::produce " << event.id().run() << ":" << event.id().event();
0168
0169   // geometry
0170   const auto &geometry = es.getData(ctppsGeometryToken);
0171   if (geometryWatcher.check(es))
0172     lrcgn->resetGeometry(&geometry);
0173
0174   // get input
0175   edm::Handle<edm::DetSetVector<TotemRPRecHit>> input;
0176   event.getByToken(detSetVectorTotemRPRecHitToken, input);
0177
0178   // prepare output
0179   DetSetVector<TotemRPUVPattern> patternsVector;
0180
0181   // split input per RP and per U/V projection
0182   struct RPData {
0183     DetSetVector<TotemRPRecHit> hits_U, hits_V;
0184     map<uint8_t, uint16_t> planeOccupancy_U, planeOccupancy_V;
0185   };
0186   map<unsigned int, RPData> rpData;
0187
0188   for (auto &ids : *input) {
0189     TotemRPDetId detId(ids.detId());
0190     unsigned int plane = detId.plane();
0191     bool uDir = detId.isStripsCoordinateUDirection();
0192
0193     CTPPSDetId rpId = detId.rpId();
0194
0195     RPData &data = rpData[rpId];
0196
0197     for (auto &h : ids) {
0198       if (uDir) {
0199         auto &ods = data.hits_U.find_or_insert(ids.detId());
0200         ods.push_back(h);
0201         data.planeOccupancy_U[plane]++;
0202       } else {
0203         auto &ods = data.hits_V.find_or_insert(ids.detId());
0204         ods.push_back(h);
0205         data.planeOccupancy_V[plane]++;
0206       }
0207     }
0208   }
0209
0210   // track recognition pot by pot
0211   for (auto const &it : rpData) {
0212     CTPPSDetId rpId(it.first);
0213     RPData const &data = it.second;
0214
0215     // merge default and exceptional settings (if available)
0216     unsigned int minPlanesPerProjectionToFit_U = minPlanesPerProjectionToFit;
0217     unsigned int minPlanesPerProjectionToFit_V = minPlanesPerProjectionToFit;
0218     double threshold_U = threshold;
0219     double threshold_V = threshold;
0220
0221     auto setIt = exceptionalSettings.find(rpId);
0222     if (setIt != exceptionalSettings.end()) {
0223       minPlanesPerProjectionToFit_U = setIt->second.minPlanesPerProjectionToFit_U;
0224       minPlanesPerProjectionToFit_V = setIt->second.minPlanesPerProjectionToFit_V;
0225       threshold_U = setIt->second.threshold_U;
0226       threshold_V = setIt->second.threshold_V;
0227     }
0228
0229     auto &uColl = data.planeOccupancy_U;
0230     auto &vColl = data.planeOccupancy_V;
0231
0232     if (verbosity > 5) {
0233       LogVerbatim("TotemRPUVPatternFinder")
0234           << "\tRP " << rpId << "\n\t\tall planes: u = " << uColl.size() << ", v = " << vColl.size();
0235     }
0236
0237     // count planes with clean data (no showers, noise, ...)
0238     unsigned int uPlanes = 0, vPlanes = 0;
0239     for (auto pit : uColl)
0240       if (pit.second <= maxHitsPerPlaneToSearch)
0241         uPlanes++;
0242
0243     for (auto pit : vColl)
0244       if (pit.second <= maxHitsPerPlaneToSearch)
0245         vPlanes++;
0246
0247     if (verbosity > 5)
0248       LogVerbatim("TotemRPUVPatternFinder") << "\t\tplanes with clean data: u = " << uPlanes << ", v = " << vPlanes;
0249
0250     // discard RPs with too few reasonable planes
0251     if (uPlanes < minPlanesPerProjectionToSearch || vPlanes < minPlanesPerProjectionToSearch)
0252       continue;
0253
0254     // prepare data containers
0255     DetSet<TotemRPUVPattern> &patterns = patternsVector.find_or_insert(rpId);
0256
0257     // "typical" z0 for the RP
0258     double z0 = geometry.rp(rpId)->translation().z();
0259
0260     // u then v recognition
0261     recognizeAndSelect(TotemRPUVPattern::projU, z0, threshold_U, minPlanesPerProjectionToFit_U, data.hits_U, patterns);
0262
0263     recognizeAndSelect(TotemRPUVPattern::projV, z0, threshold_V, minPlanesPerProjectionToFit_V, data.hits_V, patterns);
0264
0265     if (verbosity > 5) {
0266       LogVerbatim("TotemRPUVPatternFinder") << "\t\tpatterns:";
0267       for (const auto &p : patterns) {
0268         unsigned int n_hits = 0;
0269         for (auto &hds : p.hits())
0270           n_hits += hds.size();
0271
0272         LogVerbatim("TotemRPUVPatternFinder")
0273             << "\t\t\tproj = " << ((p.projection() == TotemRPUVPattern::projU) ? "U" : "V") << ", a = " << p.a()
0274             << ", b = " << p.b() << ", w = " << p.w() << ", fittable = " << p.fittable() << ", hits = " << n_hits;
0275       }
0276     }
0277   }
0278
0279   // save output
0280   event.put(make_unique<DetSetVector<TotemRPUVPattern>>(patternsVector));
0281 }
0282
0283 //----------------------------------------------------------------------------------------------------
0284
0285 void TotemRPUVPatternFinder::fillDescriptions(edm::ConfigurationDescriptions &descr) {
0286   edm::ParameterSetDescription desc;
0287
0288   desc.add<edm::InputTag>("tagRecHit", edm::InputTag("totemRPRecHitProducer"))
0289       ->setComment("input rechits collection to retrieve");
0290   desc.addUntracked<unsigned int>("verbosity", 0);
0291   desc.add<unsigned int>("maxHitsPerPlaneToSearch", 5)
0292       ->setComment("minimum threshold of hits multiplicity to flag the pattern as dirty");
0293   desc.add<unsigned int>("minPlanesPerProjectionToSearch", 3)
0294       ->setComment(
0295           "minimal number of reasonable (= not empty and not dirty) planes per projection and per RP, to start the "
0296           "pattern search");
0297   desc.add<double>("clusterSize_a", 0.02 /* rad */)
0298       ->setComment("(full) cluster size (in rad) in slope-intercept space");
0299   desc.add<double>("clusterSize_b", 0.3 /* mm */);
0300
0301   desc.add<double>("threshold", 2.99)
0302       ->setComment(
0303           "minimal weight of (Hough) cluster to accept it as candidate\n"
0304           "  weight of cluster = sum of weights of contributing points\n"
0305           "  weight of point = sigma0 / sigma_of_point\n"
0306           "most often: weight of point ~ 1, thus cluster weight is roughly number of contributing points");
0307
0308   desc.add<unsigned int>("minPlanesPerProjectionToFit", 3)
0309       ->setComment(
0310           "minimal number of planes (in the recognised patterns) per projection and per RP, to tag the candidate as "
0311           "fittable");
0312
0313   desc.add<bool>("allowAmbiguousCombination", false)
0314       ->setComment(
0315           "whether to allow combination of most significant U and V pattern, in case there several of them.\n"
0316           "don't set it to True, unless you have reason");
0317
0318   desc.add<double>("max_a_toFit", 10.0)
0319       ->setComment(
0320           "maximal angle (in any projection) to mark the candidate as fittable -> controls track parallelity with "
0321           "beam\n"
0322           "huge value -> no constraint");
0323
0324   edm::ParameterSetDescription exceptions_validator;
0325   exceptions_validator.add<unsigned int>("rpId")->setComment("RP id according to CTPPSDetId");
0326   exceptions_validator.add<unsigned int>("minPlanesPerProjectionToFit_U");
0327   exceptions_validator.add<unsigned int>("minPlanesPerProjectionToFit_V");
0328   exceptions_validator.add<double>("threshold_U");
0329   exceptions_validator.add<double>("threshold_V");
0330
0331   std::vector<edm::ParameterSet> exceptions_default;
0332   desc.addVPSet("exceptionalSettings", exceptions_validator, exceptions_default);
0333
0334   descr.add("totemRPUVPatternFinder", desc);
0335 }
0336
0337 //----------------------------------------------------------------------------------------------------
0338
0339 DEFINE_FWK_MODULE(TotemRPUVPatternFinder);