TICL/plugins/PatternRecognitionbyFastJet.cc

0001 // Author: Marco Rovere - marco.rovere@cern.ch
0002 // Date: 10/2021
0003 #include <algorithm>
0004 #include <set>
0005 #include <vector>
0006
0007 #include "tbb/task_arena.h"
0008 #include "tbb/tbb.h"
0009
0010 #include "DataFormats/Math/interface/deltaR.h"
0011 #include "DataFormats/Math/interface/LorentzVector.h"
0012 #include "DataFormats/Candidate/interface/Candidate.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 #include "PatternRecognitionbyFastJet.h"
0017
0018 #include "TrackstersPCA.h"
0019 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0020 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022
0023 #include "fastjet/ClusterSequence.hh"
0024
0025 using namespace ticl;
0026 using namespace fastjet;
0027
0028 template <typename TILES>
0029 PatternRecognitionbyFastJet<TILES>::PatternRecognitionbyFastJet(const edm::ParameterSet &conf,
0030                                                                 edm::ConsumesCollector iC)
0031     : PatternRecognitionAlgoBaseT<TILES>(conf, iC),
0032       caloGeomToken_(iC.esConsumes<CaloGeometry, CaloGeometryRecord>()),
0033       antikt_radius_(conf.getParameter<double>("antikt_radius")),
0034       minNumLayerCluster_(conf.getParameter<int>("minNumLayerCluster")),
0035       computeLocalTime_(conf.getParameter<bool>("computeLocalTime")){};
0036
0037 template <typename TILES>
0038 void PatternRecognitionbyFastJet<TILES>::buildJetAndTracksters(std::vector<PseudoJet> &fjInputs,
0039                                                                std::vector<ticl::Trackster> &result) {
0040   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Basic) {
0041     edm::LogVerbatim("PatternRecogntionbyFastJet")
0042         << "Creating FastJet with " << fjInputs.size() << " LayerClusters in input";
0043   }
0044   fastjet::ClusterSequence sequence(fjInputs, JetDefinition(antikt_algorithm, antikt_radius_));
0045   auto jets = fastjet::sorted_by_pt(sequence.inclusive_jets(0));
0046   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Basic) {
0047     edm::LogVerbatim("PatternRecogntionbyFastJet") << "FastJet produced " << jets.size() << " jets/trackster";
0048   }
0049
0050   auto trackster_idx = result.size();
0051   auto jetsSize = std::count_if(jets.begin(), jets.end(), [this](fastjet::PseudoJet jet) {
0052     return jet.constituents().size() > static_cast<unsigned int>(minNumLayerCluster_);
0053   });
0054   result.resize(trackster_idx + jetsSize);
0055
0056   for (const auto &pj : jets) {
0057     if (pj.constituents().size() > static_cast<unsigned int>(minNumLayerCluster_)) {
0058       for (const auto &component : pj.constituents()) {
0059         result[trackster_idx].vertices().push_back(component.user_index());
0060         result[trackster_idx].vertex_multiplicity().push_back(1);
0061         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Basic) {
0062           edm::LogVerbatim("PatternRecogntionbyFastJet")
0063               << "Jet has " << pj.constituents().size() << " components that are stored in trackster " << trackster_idx;
0064         }
0065       }
0066       trackster_idx++;
0067     } else {
0068       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0069         edm::LogVerbatim("PatternRecogntionbyFastJet")
0070             << "Jet with " << pj.constituents().size() << " constituents discarded since too small wrt "
0071             << minNumLayerCluster_;
0072       }
0073     }
0074   }
0075   fjInputs.clear();
0076 }
0077
0078 template <typename TILES>
0079 void PatternRecognitionbyFastJet<TILES>::makeTracksters(
0080     const typename PatternRecognitionAlgoBaseT<TILES>::Inputs &input,
0081     std::vector<Trackster> &result,
0082     std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) {
0083   // Protect from events with no seeding regions
0084   if (input.regions.empty())
0085     return;
0086
0087   edm::EventSetup const &es = input.es;
0088   const CaloGeometry &geom = es.getData(caloGeomToken_);
0089   rhtools_.setGeometry(geom);
0090
0091   constexpr auto isHFnose = std::is_same<TILES, TICLLayerTilesHFNose>::value;
0092   constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
0093   constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
0094
0095   // We need to partition the two sides of the HGCAL detector
0096   auto lastLayerPerSide = static_cast<unsigned int>(rhtools_.lastLayer(isHFnose)) - 1;
0097   unsigned int maxLayer = 2 * lastLayerPerSide - 1;
0098   std::vector<fastjet::PseudoJet> fjInputs;
0099   fjInputs.clear();
0100   for (unsigned int currentLayer = 0; currentLayer <= maxLayer; ++currentLayer) {
0101     if (currentLayer == lastLayerPerSide) {
0102       buildJetAndTracksters(fjInputs, result);
0103     }
0104     const auto &tileOnLayer = input.tiles[currentLayer];
0105     for (int ieta = 0; ieta <= nEtaBin; ++ieta) {
0106       auto offset = ieta * nPhiBin;
0107       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0108         edm::LogVerbatim("PatternRecogntionbyFastJet") << "offset: " << offset;
0109       }
0110       for (int iphi = 0; iphi <= nPhiBin; ++iphi) {
0111         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0112           edm::LogVerbatim("PatternRecogntionbyFastJet") << "iphi: " << iphi;
0113           edm::LogVerbatim("PatternRecogntionbyFastJet") << "Entries in tileBin: " << tileOnLayer[offset + iphi].size();
0114         }
0115         for (auto clusterIdx : tileOnLayer[offset + iphi]) {
0116           // Skip masked layer clusters
0117           if (input.mask[clusterIdx] == 0.) {
0118             if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0119               edm::LogVerbatim("PatternRecogntionbyFastJet") << "Skipping masked layerIdx " << clusterIdx;
0120             }
0121             continue;
0122           }
0123           // Should we correct for the position of the PV?
0124           auto const &cl = input.layerClusters[clusterIdx];
0125           math::XYZVector direction(cl.x(), cl.y(), cl.z());
0126           direction = direction.Unit();
0127           direction *= cl.energy();
0128           auto fpj = fastjet::PseudoJet(direction.X(), direction.Y(), direction.Z(), cl.energy());
0129           fpj.set_user_index(clusterIdx);
0130           fjInputs.push_back(fpj);
0131         }  // End of loop on the clusters on currentLayer
0132       }  // End of loop over phi-bin region
0133     }  // End of loop over eta-bin region
0134   }  // End of loop over layers
0135
0136   // Collect the jet from the other side wrt to the one taken care of inside the main loop above.
0137   buildJetAndTracksters(fjInputs, result);
0138
0139   ticl::assignPCAtoTracksters(result,
0140                               input.layerClusters,
0141                               input.layerClustersTime,
0142                               rhtools_.getPositionLayer(rhtools_.lastLayerEE(isHFnose), isHFnose).z(),
0143                               rhtools_,
0144                               computeLocalTime_);
0145
0146   // run energy regression and ID
0147   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Basic) {
0148     for (auto const &t : result) {
0149       edm::LogVerbatim("PatternRecogntionbyFastJet") << "Barycenter: " << t.barycenter();
0150       edm::LogVerbatim("PatternRecogntionbyFastJet") << "LCs: " << t.vertices().size();
0151       edm::LogVerbatim("PatternRecogntionbyFastJet") << "Energy: " << t.raw_energy();
0152       edm::LogVerbatim("PatternRecogntionbyFastJet") << "Regressed: " << t.regressed_energy();
0153     }
0154   }
0155 }
0156
0157 template <typename TILES>
0158 void PatternRecognitionbyFastJet<TILES>::filter(std::vector<Trackster> &output,
0159                                                 const std::vector<Trackster> &inTracksters,
0160                                                 const typename PatternRecognitionAlgoBaseT<TILES>::Inputs &input,
0161                                                 std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) {
0162   output = inTracksters;
0163 }
0164
0165 template <typename TILES>
0166 void PatternRecognitionbyFastJet<TILES>::fillPSetDescription(edm::ParameterSetDescription &iDesc) {
0167   iDesc.add<int>("algo_verbosity", 0);
0168   iDesc.add<double>("antikt_radius", 0.09)->setComment("Radius to be used while running the Anti-kt clustering");
0169   iDesc.add<int>("minNumLayerCluster", 5)->setComment("Not Inclusive");
0170   iDesc.add<bool>("computeLocalTime", false);
0171 }
0172
0173 template class ticl::PatternRecognitionbyFastJet<TICLLayerTiles>;