TICL/plugins/SeedingRegionByTracks.cc

0001 // Author: Arabella Martelli, Felice Pantaleo, Marco Rovere
0002 // arabella.martelli@cern.ch, felice.pantaleo@cern.ch, marco.rovere@cern.ch
0003 // Date: 06/2019
0004 #include <algorithm>
0005 #include <set>
0006 #include <vector>
0007
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "SeedingRegionByTracks.h"
0010 #include "DataFormats/TrackReco/interface/Track.h"
0011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0012 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0013
0014 using namespace ticl;
0015
0016 SeedingRegionByTracks::SeedingRegionByTracks(const edm::ParameterSet &conf, edm::ConsumesCollector &sumes)
0017     : SeedingRegionAlgoBase(conf, sumes),
0018       tracks_token_(sumes.consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("tracks"))),
0019       cutTk_(conf.getParameter<std::string>("cutTk")),
0020       detector_(conf.getParameter<std::string>("detector")),
0021       propName_(conf.getParameter<std::string>("propagator")),
0022       bfield_token_(sumes.esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>()),
0023       propagator_token_(sumes.esConsumes<Propagator, TrackingComponentsRecord, edm::Transition::BeginRun>(
0024           edm::ESInputTag("", propName_))) {
0025   std::string detectorName_ = (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive";
0026   hdc_token_ = sumes.esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0027       edm::ESInputTag("", detectorName_));
0028 }
0029
0030 SeedingRegionByTracks::~SeedingRegionByTracks() {}
0031
0032 void SeedingRegionByTracks::initialize(const edm::EventSetup &es) {
0033   edm::ESHandle<HGCalDDDConstants> hdc = es.getHandle(hdc_token_);
0034   hgcons_ = hdc.product();
0035
0036   buildFirstLayers();
0037
0038   bfield_ = es.getHandle(bfield_token_);
0039   propagator_ = es.getHandle(propagator_token_);
0040 }
0041
0042 void SeedingRegionByTracks::makeRegions(const edm::Event &ev,
0043                                         const edm::EventSetup &es,
0044                                         std::vector<TICLSeedingRegion> &result) {
0045   edm::Handle<reco::TrackCollection> tracks_h;
0046   ev.getByToken(tracks_token_, tracks_h);
0047   edm::ProductID trkId = tracks_h.id();
0048   auto bFieldProd = bfield_.product();
0049   const Propagator &prop = (*propagator_);
0050
0051   int nTracks = tracks_h->size();
0052   for (int i = 0; i < nTracks; ++i) {
0053     const reco::Track &tk = (*tracks_h)[i];
0054     if (!cutTk_((tk))) {
0055       continue;
0056     }
0057
0058     FreeTrajectoryState fts = trajectoryStateTransform::outerFreeState((tk), bFieldProd);
0059     int iSide = int(tk.eta() > 0);
0060     TrajectoryStateOnSurface tsos = prop.propagate(fts, firstDisk_[iSide]->surface());
0061     if (tsos.isValid()) {
0062       result.emplace_back(tsos.globalPosition(), tsos.globalMomentum(), iSide, i, trkId);
0063     }
0064   }
0065   // sorting seeding region by descending momentum
0066   std::sort(result.begin(), result.end(), [](const TICLSeedingRegion &a, const TICLSeedingRegion &b) {
0067     return a.directionAtOrigin.perp2() > b.directionAtOrigin.perp2();
0068   });
0069 }
0070
0071 void SeedingRegionByTracks::fillPSetDescription(edm::ParameterSetDescription &desc) {
0072   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0073   desc.add<std::string>("cutTk",
0074                         "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && "
0075                         "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
0076   desc.add<std::string>("propagator", "PropagatorWithMaterial");
0077   desc.add<std::string>("detector", "HGCAL");
0078   SeedingRegionAlgoBase::fillPSetDescription(desc);
0079 }
0080
0081 void SeedingRegionByTracks::buildFirstLayers() {
0082   float zVal = hgcons_->waferZ(1, true);
0083   std::pair<double, double> rMinMax = hgcons_->rangeR(zVal, true);
0084
0085   for (int iSide = 0; iSide < 2; ++iSide) {
0086     float zSide = (iSide == 0) ? (-1. * zVal) : zVal;
0087     firstDisk_[iSide] =
0088         std::make_unique<GeomDet>(Disk::build(Disk::PositionType(0, 0, zSide),
0089                                               Disk::RotationType(),
0090                                               SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5))
0091                                       .get());
0092   }
0093 }