File indexing completed on 2023-06-02 00:50:51
0001 #include <memory> // unique_ptr
0002 #include "FWCore/Framework/interface/stream/EDProducer.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0005 #include "FWCore/ParameterSet/interface/PluginDescription.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Utilities/interface/ESGetToken.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0012 #include "FWCore/Framework/interface/ConsumesCollector.h"
0013
0014 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0015 #include "DataFormats/HGCalReco/interface/Common.h"
0016 #include "DataFormats/HGCalReco/interface/TICLLayerTile.h"
0017 #include "DataFormats/HGCalReco/interface/Trackster.h"
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019 #include "DataFormats/MuonReco/interface/Muon.h"
0020 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0021 #include "DataFormats/HGCalReco/interface/TICLCandidate.h"
0022 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0023 #include "DataFormats/Math/interface/Vector3D.h"
0024
0025 #include "RecoHGCal/TICL/interface/GlobalCache.h"
0026
0027 #include "PhysicsTools/TensorFlow/interface/TfGraphRecord.h"
0028 #include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
0029 #include "PhysicsTools/TensorFlow/interface/TfGraphDefWrapper.h"
0030
0031 #include "RecoHGCal/TICL/plugins/LinkingAlgoBase.h"
0032 #include "RecoHGCal/TICL/plugins/LinkingAlgoFactory.h"
0033 #include "RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h"
0034
0035 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0036 #include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
0037
0038 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0039 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0040 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0041
0042 #include "MagneticField/Engine/interface/MagneticField.h"
0043 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0044
0045 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0046 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0047 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0048
0049 #include "TrackstersPCA.h"
0050
0051 using namespace ticl;
0052
0053 class TrackstersMergeProducer : public edm::stream::EDProducer<> {
0054 public:
0055 explicit TrackstersMergeProducer(const edm::ParameterSet &ps);
0056 ~TrackstersMergeProducer() override{};
0057 void produce(edm::Event &, const edm::EventSetup &) override;
0058 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0059
0060
0061 static std::unique_ptr<TrackstersCache> initializeGlobalCache(const edm::ParameterSet &);
0062 static void globalEndJob(TrackstersCache *);
0063
0064 void beginJob();
0065 void endJob();
0066
0067 void beginRun(edm::Run const &iEvent, edm::EventSetup const &es) override;
0068
0069 private:
0070 typedef ticl::Trackster::IterationIndex TracksterIterIndex;
0071 typedef math::XYZVector Vector;
0072
0073 void fillTile(TICLTracksterTiles &, const std::vector<Trackster> &, TracksterIterIndex);
0074
0075 void energyRegressionAndID(const std::vector<reco::CaloCluster> &layerClusters,
0076 const tensorflow::Session *,
0077 std::vector<Trackster> &result) const;
0078 void printTrackstersDebug(const std::vector<Trackster> &, const char *label) const;
0079 void assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates) const;
0080 void dumpTrackster(const Trackster &) const;
0081
0082 std::unique_ptr<LinkingAlgoBase> linkingAlgo_;
0083
0084 const edm::EDGetTokenT<std::vector<Trackster>> tracksters_clue3d_token_;
0085 const edm::EDGetTokenT<std::vector<reco::CaloCluster>> clusters_token_;
0086 const edm::EDGetTokenT<edm::ValueMap<std::pair<float, float>>> clustersTime_token_;
0087 const edm::EDGetTokenT<std::vector<reco::Track>> tracks_token_;
0088 const edm::EDGetTokenT<edm::ValueMap<float>> tracks_time_token_;
0089 const edm::EDGetTokenT<edm::ValueMap<float>> tracks_time_quality_token_;
0090 const edm::EDGetTokenT<edm::ValueMap<float>> tracks_time_err_token_;
0091 const edm::EDGetTokenT<std::vector<reco::Muon>> muons_token_;
0092 const std::string tfDnnLabel_;
0093 const edm::ESGetToken<TfGraphDefWrapper, TfGraphRecord> tfDnnToken_;
0094 const tensorflow::Session *tfSession_;
0095
0096 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometry_token_;
0097 const std::string detector_;
0098 const std::string propName_;
0099
0100 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> bfield_token_;
0101 const edm::ESGetToken<Propagator, TrackingComponentsRecord> propagator_token_;
0102 const bool optimiseAcrossTracksters_;
0103 const int eta_bin_window_;
0104 const int phi_bin_window_;
0105 const double pt_sigma_high_;
0106 const double pt_sigma_low_;
0107 const double halo_max_distance2_;
0108 const double track_min_pt_;
0109 const double track_min_eta_;
0110 const double track_max_eta_;
0111 const int track_max_missing_outerhits_;
0112 const double cosangle_align_;
0113 const double e_over_h_threshold_;
0114 const double pt_neutral_threshold_;
0115 const double resol_calo_offset_had_;
0116 const double resol_calo_scale_had_;
0117 const double resol_calo_offset_em_;
0118 const double resol_calo_scale_em_;
0119 const std::string eidInputName_;
0120 const std::string eidOutputNameEnergy_;
0121 const std::string eidOutputNameId_;
0122 const float eidMinClusterEnergy_;
0123 const int eidNLayers_;
0124 const int eidNClusters_;
0125 std::once_flag initializeGeometry_;
0126
0127 const HGCalDDDConstants *hgcons_;
0128
0129 std::unique_ptr<GeomDet> firstDisk_[2];
0130
0131 tensorflow::Session *eidSession_;
0132 hgcal::RecHitTools rhtools_;
0133
0134 static constexpr int eidNFeatures_ = 3;
0135
0136 edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> hdc_token_;
0137 };
0138
0139 TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps)
0140 : tracksters_clue3d_token_(consumes<std::vector<Trackster>>(ps.getParameter<edm::InputTag>("trackstersclue3d"))),
0141 clusters_token_(consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("layer_clusters"))),
0142 clustersTime_token_(
0143 consumes<edm::ValueMap<std::pair<float, float>>>(ps.getParameter<edm::InputTag>("layer_clustersTime"))),
0144 tracks_token_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("tracks"))),
0145 tracks_time_token_(consumes<edm::ValueMap<float>>(ps.getParameter<edm::InputTag>("tracksTime"))),
0146 tracks_time_quality_token_(consumes<edm::ValueMap<float>>(ps.getParameter<edm::InputTag>("tracksTimeQual"))),
0147 tracks_time_err_token_(consumes<edm::ValueMap<float>>(ps.getParameter<edm::InputTag>("tracksTimeErr"))),
0148 muons_token_(consumes<std::vector<reco::Muon>>(ps.getParameter<edm::InputTag>("muons"))),
0149 tfDnnLabel_(ps.getParameter<std::string>("tfDnnLabel")),
0150 tfDnnToken_(esConsumes(edm::ESInputTag("", tfDnnLabel_))),
0151 tfSession_(nullptr),
0152 geometry_token_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0153 detector_(ps.getParameter<std::string>("detector")),
0154 propName_(ps.getParameter<std::string>("propagator")),
0155 bfield_token_(esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>()),
0156 propagator_token_(
0157 esConsumes<Propagator, TrackingComponentsRecord, edm::Transition::BeginRun>(edm::ESInputTag("", propName_))),
0158 optimiseAcrossTracksters_(ps.getParameter<bool>("optimiseAcrossTracksters")),
0159 eta_bin_window_(ps.getParameter<int>("eta_bin_window")),
0160 phi_bin_window_(ps.getParameter<int>("phi_bin_window")),
0161 pt_sigma_high_(ps.getParameter<double>("pt_sigma_high")),
0162 pt_sigma_low_(ps.getParameter<double>("pt_sigma_low")),
0163 halo_max_distance2_(ps.getParameter<double>("halo_max_distance2")),
0164 track_min_pt_(ps.getParameter<double>("track_min_pt")),
0165 track_min_eta_(ps.getParameter<double>("track_min_eta")),
0166 track_max_eta_(ps.getParameter<double>("track_max_eta")),
0167 track_max_missing_outerhits_(ps.getParameter<int>("track_max_missing_outerhits")),
0168 cosangle_align_(ps.getParameter<double>("cosangle_align")),
0169 e_over_h_threshold_(ps.getParameter<double>("e_over_h_threshold")),
0170 pt_neutral_threshold_(ps.getParameter<double>("pt_neutral_threshold")),
0171 resol_calo_offset_had_(ps.getParameter<double>("resol_calo_offset_had")),
0172 resol_calo_scale_had_(ps.getParameter<double>("resol_calo_scale_had")),
0173 resol_calo_offset_em_(ps.getParameter<double>("resol_calo_offset_em")),
0174 resol_calo_scale_em_(ps.getParameter<double>("resol_calo_scale_em")),
0175 eidInputName_(ps.getParameter<std::string>("eid_input_name")),
0176 eidOutputNameEnergy_(ps.getParameter<std::string>("eid_output_name_energy")),
0177 eidOutputNameId_(ps.getParameter<std::string>("eid_output_name_id")),
0178 eidMinClusterEnergy_(ps.getParameter<double>("eid_min_cluster_energy")),
0179 eidNLayers_(ps.getParameter<int>("eid_n_layers")),
0180 eidNClusters_(ps.getParameter<int>("eid_n_clusters")),
0181 eidSession_(nullptr) {
0182 produces<std::vector<Trackster>>();
0183 produces<std::vector<TICLCandidate>>();
0184
0185 std::string detectorName_ = (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive";
0186 hdc_token_ =
0187 esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag("", detectorName_));
0188
0189 auto linkingPSet = ps.getParameter<edm::ParameterSet>("linkingPSet");
0190 auto algoType = linkingPSet.getParameter<std::string>("type");
0191 linkingAlgo_ = LinkingAlgoFactory::get()->create(algoType, linkingPSet);
0192 }
0193
0194 void TrackstersMergeProducer::beginJob() {}
0195
0196 void TrackstersMergeProducer::endJob(){};
0197
0198 void TrackstersMergeProducer::beginRun(edm::Run const &iEvent, edm::EventSetup const &es) {
0199 edm::ESHandle<HGCalDDDConstants> hdc = es.getHandle(hdc_token_);
0200 hgcons_ = hdc.product();
0201
0202 edm::ESHandle<CaloGeometry> geom = es.getHandle(geometry_token_);
0203 rhtools_.setGeometry(*geom);
0204
0205 edm::ESHandle<MagneticField> bfield = es.getHandle(bfield_token_);
0206 edm::ESHandle<Propagator> propagator = es.getHandle(propagator_token_);
0207
0208 linkingAlgo_->initialize(hgcons_, rhtools_, bfield, propagator);
0209 };
0210
0211 void TrackstersMergeProducer::fillTile(TICLTracksterTiles &tracksterTile,
0212 const std::vector<Trackster> &tracksters,
0213 TracksterIterIndex tracksterIteration) {
0214 int tracksterId = 0;
0215 for (auto const &t : tracksters) {
0216 tracksterTile.fill(tracksterIteration, t.barycenter().eta(), t.barycenter().phi(), tracksterId);
0217 LogDebug("TrackstersMergeProducer") << "Adding tracksterId: " << tracksterId << " into bin [eta,phi]: [ "
0218 << tracksterTile[tracksterIteration].etaBin(t.barycenter().eta()) << ", "
0219 << tracksterTile[tracksterIteration].phiBin(t.barycenter().phi())
0220 << "] for iteration: " << tracksterIteration << std::endl;
0221
0222 tracksterId++;
0223 }
0224 }
0225
0226 void TrackstersMergeProducer::dumpTrackster(const Trackster &t) const {
0227 auto e_over_h = (t.raw_em_pt() / ((t.raw_pt() - t.raw_em_pt()) != 0. ? (t.raw_pt() - t.raw_em_pt()) : 1.));
0228 LogDebug("TrackstersMergeProducer")
0229 << "\nTrackster raw_pt: " << t.raw_pt() << " raw_em_pt: " << t.raw_em_pt() << " eoh: " << e_over_h
0230 << " barycenter: " << t.barycenter() << " eta,phi (baricenter): " << t.barycenter().eta() << ", "
0231 << t.barycenter().phi() << " eta,phi (eigen): " << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi()
0232 << " pt(eigen): " << std::sqrt(t.eigenvectors(0).Unit().perp2()) * t.raw_energy() << " seedID: " << t.seedID()
0233 << " seedIndex: " << t.seedIndex() << " size: " << t.vertices().size() << " average usage: "
0234 << (std::accumulate(std::begin(t.vertex_multiplicity()), std::end(t.vertex_multiplicity()), 0.) /
0235 (float)t.vertex_multiplicity().size())
0236 << " raw_energy: " << t.raw_energy() << " regressed energy: " << t.regressed_energy()
0237 << " probs(ga/e/mu/np/cp/nh/am/unk): ";
0238 for (auto const &p : t.id_probabilities()) {
0239 LogDebug("TrackstersMergeProducer") << std::fixed << p << " ";
0240 }
0241 LogDebug("TrackstersMergeProducer") << " sigmas: ";
0242 for (auto const &s : t.sigmas()) {
0243 LogDebug("TrackstersMergeProducer") << s << " ";
0244 }
0245 LogDebug("TrackstersMergeProducer") << std::endl;
0246 }
0247
0248 void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es) {
0249 auto resultTrackstersMerged = std::make_unique<std::vector<Trackster>>();
0250 auto resultCandidates = std::make_unique<std::vector<TICLCandidate>>();
0251 auto resultFromTracks = std::make_unique<std::vector<TICLCandidate>>();
0252 tfSession_ = es.getData(tfDnnToken_).getSession();
0253
0254 edm::Handle<std::vector<Trackster>> trackstersclue3d_h;
0255 evt.getByToken(tracksters_clue3d_token_, trackstersclue3d_h);
0256
0257 edm::Handle<std::vector<reco::Track>> track_h;
0258 evt.getByToken(tracks_token_, track_h);
0259 const auto &tracks = *track_h;
0260
0261 const auto &layerClusters = evt.get(clusters_token_);
0262 const auto &layerClustersTimes = evt.get(clustersTime_token_);
0263 const auto &muons = evt.get(muons_token_);
0264 const auto &trackTime = evt.get(tracks_time_token_);
0265 const auto &trackTimeErr = evt.get(tracks_time_err_token_);
0266 const auto &trackTimeQual = evt.get(tracks_time_quality_token_);
0267
0268
0269 linkingAlgo_->linkTracksters(
0270 track_h, trackTime, trackTimeErr, trackTimeQual, muons, trackstersclue3d_h, *resultCandidates, *resultFromTracks);
0271
0272
0273 LogDebug("TrackstersMergeProducer") << "Results from the linking step : " << std::endl
0274 << "No. of Tracks : " << tracks.size()
0275 << " No. of Tracksters : " << (*trackstersclue3d_h).size() << std::endl
0276 << "(neutral candidates have track id -1)" << std::endl;
0277
0278 std::vector<TICLCandidate> &candidates = *resultCandidates;
0279 for (const auto &cand : candidates) {
0280 auto track_ptr = cand.trackPtr();
0281 auto trackster_ptrs = cand.tracksters();
0282
0283 #ifdef EDM_ML_DEBUG
0284 auto track_idx = track_ptr.get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
0285 track_idx = (track_ptr.isNull()) ? -1 : track_idx;
0286 LogDebug("TrackstersMergeProducer") << "PDG ID " << cand.pdgId() << " charge " << cand.charge() << " p " << cand.p()
0287 << std::endl;
0288 LogDebug("TrackstersMergeProducer") << "track id (p) : " << track_idx << " ("
0289 << (track_ptr.isNull() ? -1 : track_ptr->p()) << ") "
0290 << " trackster ids (E) : ";
0291 #endif
0292
0293
0294 ticl::Trackster outTrackster;
0295 auto updated_size = 0;
0296 for (const auto &ts_ptr : trackster_ptrs) {
0297 #ifdef EDM_ML_DEBUG
0298 auto ts_idx = ts_ptr.get() - (edm::Ptr<ticl::Trackster>(trackstersclue3d_h, 0)).get();
0299 LogDebug("TrackstersMergeProducer") << ts_idx << " (" << ts_ptr->raw_energy() << ") ";
0300 #endif
0301
0302 auto &thisTrackster = *ts_ptr;
0303 updated_size += thisTrackster.vertices().size();
0304 outTrackster.vertices().reserve(updated_size);
0305 outTrackster.vertex_multiplicity().reserve(updated_size);
0306 std::copy(std::begin(thisTrackster.vertices()),
0307 std::end(thisTrackster.vertices()),
0308 std::back_inserter(outTrackster.vertices()));
0309 std::copy(std::begin(thisTrackster.vertex_multiplicity()),
0310 std::end(thisTrackster.vertex_multiplicity()),
0311 std::back_inserter(outTrackster.vertex_multiplicity()));
0312 }
0313
0314 LogDebug("TrackstersMergeProducer") << std::endl;
0315
0316
0317 auto &orig_vtx = outTrackster.vertices();
0318 auto vtx_sorted{orig_vtx};
0319 std::sort(std::begin(vtx_sorted), std::end(vtx_sorted));
0320 for (unsigned int iLC = 1; iLC < vtx_sorted.size(); ++iLC) {
0321 if (vtx_sorted[iLC] == vtx_sorted[iLC - 1]) {
0322
0323 const auto lcIdx = vtx_sorted[iLC];
0324 const auto firstEl = std::find(orig_vtx.begin(), orig_vtx.end(), lcIdx);
0325 const auto firstPos = std::distance(std::begin(orig_vtx), firstEl);
0326 auto iDup = std::find(std::next(firstEl), orig_vtx.end(), lcIdx);
0327 while (iDup != orig_vtx.end()) {
0328 orig_vtx.erase(iDup);
0329 outTrackster.vertex_multiplicity().erase(outTrackster.vertex_multiplicity().begin() +
0330 std::distance(std::begin(orig_vtx), iDup));
0331 outTrackster.vertex_multiplicity()[firstPos] -= 1;
0332 iDup = std::find(std::next(firstEl), orig_vtx.end(), lcIdx);
0333 };
0334 }
0335 }
0336
0337 outTrackster.zeroProbabilities();
0338 if (!track_ptr.isNull())
0339 outTrackster.setSeed(track_h.id(), track_ptr.get() - (edm::Ptr<reco::Track>(track_h, 0)).get());
0340 if (!outTrackster.vertices().empty()) {
0341 resultTrackstersMerged->push_back(outTrackster);
0342 }
0343 }
0344
0345 assignPCAtoTracksters(*resultTrackstersMerged,
0346 layerClusters,
0347 layerClustersTimes,
0348 rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z());
0349 energyRegressionAndID(layerClusters, tfSession_, *resultTrackstersMerged);
0350
0351
0352 assert(resultTrackstersMerged->size() == resultCandidates->size());
0353
0354 auto isHad = [](const Trackster &tracksterMerge) {
0355 return tracksterMerge.id_probability(Trackster::ParticleType::photon) +
0356 tracksterMerge.id_probability(Trackster::ParticleType::electron) <
0357 0.5;
0358 };
0359 for (size_t i = 0; i < resultTrackstersMerged->size(); i++) {
0360 auto const &tm = (*resultTrackstersMerged)[i];
0361 auto &cand = (*resultCandidates)[i];
0362
0363 cand.setIdProbabilities(tm.id_probabilities());
0364
0365 if (!cand.trackPtr().isNull()) {
0366 auto pdgId = isHad(tm) ? 211 : 11;
0367 auto const &tk = cand.trackPtr().get();
0368 cand.setPdgId(pdgId * tk->charge());
0369 cand.setCharge(tk->charge());
0370 cand.setRawEnergy(tm.raw_energy());
0371 auto const ®rE = tm.regressed_energy();
0372 math::XYZTLorentzVector p4(regrE * tk->momentum().unit().x(),
0373 regrE * tk->momentum().unit().y(),
0374 regrE * tk->momentum().unit().z(),
0375 regrE);
0376 cand.setP4(p4);
0377 } else {
0378 auto pdgId = isHad(tm) ? 130 : 22;
0379 cand.setPdgId(pdgId);
0380 cand.setCharge(0);
0381 cand.setRawEnergy(tm.raw_energy());
0382 const float ®rE = tm.regressed_energy();
0383 math::XYZTLorentzVector p4(regrE * tm.barycenter().unit().x(),
0384 regrE * tm.barycenter().unit().y(),
0385 regrE * tm.barycenter().unit().z(),
0386 regrE);
0387 cand.setP4(p4);
0388 }
0389 }
0390 for (auto &cand : *resultFromTracks) {
0391 auto const &tk = cand.trackPtr().get();
0392 cand.setPdgId(211 * tk->charge());
0393 cand.setCharge(tk->charge());
0394 const float energy = std::sqrt(tk->p() * tk->p() + ticl::mpion2);
0395 cand.setRawEnergy(energy);
0396 math::PtEtaPhiMLorentzVector p4Polar(tk->pt(), tk->eta(), tk->phi(), ticl::mpion);
0397 cand.setP4(p4Polar);
0398 }
0399
0400 resultCandidates->insert(resultCandidates->end(), resultFromTracks->begin(), resultFromTracks->end());
0401 assignTimeToCandidates(*resultCandidates);
0402
0403 evt.put(std::move(resultTrackstersMerged));
0404 evt.put(std::move(resultCandidates));
0405 }
0406
0407 void TrackstersMergeProducer::energyRegressionAndID(const std::vector<reco::CaloCluster> &layerClusters,
0408 const tensorflow::Session *eidSession,
0409 std::vector<Trackster> &tracksters) const {
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432 int batchSize = (int)tracksters.size();
0433 if (batchSize == 0) {
0434 return;
0435 }
0436
0437 for (auto &t : tracksters) {
0438 t.setRegressedEnergy(0.f);
0439 t.zeroProbabilities();
0440 }
0441
0442
0443 tensorflow::TensorShape shape({batchSize, eidNLayers_, eidNClusters_, eidNFeatures_});
0444 tensorflow::Tensor input(tensorflow::DT_FLOAT, shape);
0445 tensorflow::NamedTensorList inputList = {{eidInputName_, input}};
0446 static constexpr int inputDimension = 4;
0447
0448 std::vector<tensorflow::Tensor> outputs;
0449 std::vector<std::string> outputNames;
0450 if (!eidOutputNameEnergy_.empty()) {
0451 outputNames.push_back(eidOutputNameEnergy_);
0452 }
0453 if (!eidOutputNameId_.empty()) {
0454 outputNames.push_back(eidOutputNameId_);
0455 }
0456
0457
0458 for (int i = 0; i < batchSize; i++) {
0459 const Trackster &trackster = tracksters[i];
0460
0461
0462
0463
0464
0465 std::vector<int> clusterIndices(trackster.vertices().size());
0466 for (int k = 0; k < (int)trackster.vertices().size(); k++) {
0467 clusterIndices[k] = k;
0468 }
0469 sort(clusterIndices.begin(), clusterIndices.end(), [&layerClusters, &trackster](const int &a, const int &b) {
0470 return layerClusters[trackster.vertices(a)].energy() > layerClusters[trackster.vertices(b)].energy();
0471 });
0472
0473
0474 std::vector<int> seenClusters(eidNLayers_);
0475
0476
0477 for (const int &k : clusterIndices) {
0478
0479 const reco::CaloCluster &cluster = layerClusters[trackster.vertices(k)];
0480 int j = rhtools_.getLayerWithOffset(cluster.hitsAndFractions()[0].first) - 1;
0481 if (j < eidNLayers_ && seenClusters[j] < eidNClusters_) {
0482
0483 float *features = &input.tensor<float, inputDimension>()(i, j, seenClusters[j], 0);
0484
0485
0486 *(features++) = float(cluster.energy() / float(trackster.vertex_multiplicity(k)));
0487 *(features++) = float(std::abs(cluster.eta()));
0488 *(features) = float(cluster.phi());
0489
0490
0491 seenClusters[j]++;
0492 }
0493 }
0494
0495
0496 for (int j = 0; j < eidNLayers_; j++) {
0497 for (int k = seenClusters[j]; k < eidNClusters_; k++) {
0498 float *features = &input.tensor<float, inputDimension>()(i, j, k, 0);
0499 for (int l = 0; l < eidNFeatures_; l++) {
0500 *(features++) = 0.f;
0501 }
0502 }
0503 }
0504 }
0505
0506
0507 tensorflow::run(eidSession, inputList, outputNames, &outputs);
0508
0509
0510 if (!eidOutputNameEnergy_.empty()) {
0511
0512 float *energy = outputs[0].flat<float>().data();
0513
0514 for (int i = 0; i < batchSize; ++i) {
0515 float regressedEnergy =
0516 tracksters[i].raw_energy() > eidMinClusterEnergy_ ? energy[i] : tracksters[i].raw_energy();
0517 tracksters[i].setRegressedEnergy(regressedEnergy);
0518 }
0519 }
0520
0521
0522 if (!eidOutputNameId_.empty()) {
0523
0524 int probsIdx = !eidOutputNameEnergy_.empty();
0525 float *probs = outputs[probsIdx].flat<float>().data();
0526 int probsNumber = tracksters[0].id_probabilities().size();
0527 for (int i = 0; i < batchSize; ++i) {
0528 tracksters[i].setProbabilities(&probs[i * probsNumber]);
0529 }
0530 }
0531 }
0532
0533 void TrackstersMergeProducer::assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates) const {
0534 for (auto &cand : resultCandidates) {
0535 if (cand.tracksters().size() > 1) {
0536 float time = 0.f;
0537 float timeErr = 0.f;
0538 for (const auto &tr : cand.tracksters()) {
0539 if (tr->timeError() > 0) {
0540 auto invTimeESq = pow(tr->timeError(), -2);
0541 time += tr->time() * invTimeESq;
0542 timeErr += invTimeESq;
0543 }
0544 }
0545 if (timeErr > 0) {
0546 timeErr = 1. / timeErr;
0547
0548 cand.setTime(time * timeErr);
0549 cand.setTimeError(sqrt(timeErr));
0550 }
0551 }
0552 }
0553 }
0554
0555 void TrackstersMergeProducer::printTrackstersDebug(const std::vector<Trackster> &tracksters, const char *label) const {
0556 #ifdef EDM_ML_DEBUG
0557 int counter = 0;
0558 for (auto const &t : tracksters) {
0559 LogDebug("TrackstersMergeProducer")
0560 << counter++ << " TrackstersMergeProducer (" << label << ") obj barycenter: " << t.barycenter()
0561 << " eta,phi (baricenter): " << t.barycenter().eta() << ", " << t.barycenter().phi()
0562 << " eta,phi (eigen): " << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi()
0563 << " pt(eigen): " << std::sqrt(t.eigenvectors(0).Unit().perp2()) * t.raw_energy() << " seedID: " << t.seedID()
0564 << " seedIndex: " << t.seedIndex() << " size: " << t.vertices().size() << " average usage: "
0565 << (std::accumulate(std::begin(t.vertex_multiplicity()), std::end(t.vertex_multiplicity()), 0.) /
0566 (float)t.vertex_multiplicity().size())
0567 << " raw_energy: " << t.raw_energy() << " regressed energy: " << t.regressed_energy()
0568 << " probs(ga/e/mu/np/cp/nh/am/unk): ";
0569 for (auto const &p : t.id_probabilities()) {
0570 LogDebug("TrackstersMergeProducer") << std::fixed << p << " ";
0571 }
0572 LogDebug("TrackstersMergeProducer") << " sigmas: ";
0573 for (auto const &s : t.sigmas()) {
0574 LogDebug("TrackstersMergeProducer") << s << " ";
0575 }
0576 LogDebug("TrackstersMergeProducer") << std::endl;
0577 }
0578 #endif
0579 }
0580
0581 void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0582 edm::ParameterSetDescription desc;
0583
0584 edm::ParameterSetDescription linkingDesc;
0585 linkingDesc.addNode(edm::PluginDescription<LinkingAlgoFactory>("type", "LinkingAlgoByDirectionGeometric", true));
0586 desc.add<edm::ParameterSetDescription>("linkingPSet", linkingDesc);
0587
0588 desc.add<edm::InputTag>("trackstersclue3d", edm::InputTag("ticlTrackstersCLUE3DHigh"));
0589 desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalMergeLayerClusters"));
0590 desc.add<edm::InputTag>("layer_clustersTime", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster"));
0591 desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0592 desc.add<edm::InputTag>("tracksTime", edm::InputTag("tofPID:t0"));
0593 desc.add<edm::InputTag>("tracksTimeQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
0594 desc.add<edm::InputTag>("tracksTimeErr", edm::InputTag("tofPID:sigmat0"));
0595 desc.add<edm::InputTag>("muons", edm::InputTag("muons1stStep"));
0596 desc.add<std::string>("detector", "HGCAL");
0597 desc.add<std::string>("propagator", "PropagatorWithMaterial");
0598 desc.add<bool>("optimiseAcrossTracksters", true);
0599 desc.add<int>("eta_bin_window", 1);
0600 desc.add<int>("phi_bin_window", 1);
0601 desc.add<double>("pt_sigma_high", 2.);
0602 desc.add<double>("pt_sigma_low", 2.);
0603 desc.add<double>("halo_max_distance2", 4.);
0604 desc.add<double>("track_min_pt", 1.);
0605 desc.add<double>("track_min_eta", 1.48);
0606 desc.add<double>("track_max_eta", 3.);
0607 desc.add<int>("track_max_missing_outerhits", 5);
0608 desc.add<double>("cosangle_align", 0.9945);
0609 desc.add<double>("e_over_h_threshold", 1.);
0610 desc.add<double>("pt_neutral_threshold", 2.);
0611 desc.add<double>("resol_calo_offset_had", 1.5);
0612 desc.add<double>("resol_calo_scale_had", 0.15);
0613 desc.add<double>("resol_calo_offset_em", 1.5);
0614 desc.add<double>("resol_calo_scale_em", 0.15);
0615 desc.add<std::string>("tfDnnLabel", "tracksterSelectionTf");
0616 desc.add<std::string>("eid_input_name", "input");
0617 desc.add<std::string>("eid_output_name_energy", "output/regressed_energy");
0618 desc.add<std::string>("eid_output_name_id", "output/id_probabilities");
0619 desc.add<double>("eid_min_cluster_energy", 2.5);
0620 desc.add<int>("eid_n_layers", 50);
0621 desc.add<int>("eid_n_clusters", 10);
0622 descriptions.add("trackstersMergeProducer", desc);
0623 }
0624
0625 DEFINE_FWK_MODULE(TrackstersMergeProducer);