Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:20:10 

0001 // -*- C++ -*-
0002 //
0003 // Package: L1CaloTrigger
0004 // Class: Phase1L1TJetProducer
0005 //
0006 /**\class Phase1L1TJetProducer Phase1L1TJetProducer.cc L1Trigger/L1CaloTrigger/plugin/Phase1L1TJetProducer.cc
0007 
0008 Description: Produces jets with a phase-1 like sliding window algorithm using a collection of reco::Candidates in input.  Also calculates MET from the histogram used to find the jets.
0009 
0010 *** INPUT PARAMETERS ***
0011   * etaBinning: vdouble with eta binning (allows non-homogeneous binning in eta)
0012   * nBinsPhi: uint32, number of bins in phi
0013   * phiLow: double, min phi (typically -pi)
0014   * phiUp: double, max phi (typically +pi)
0015   * jetIEtaSize: uint32, jet cluster size in ieta
0016   * jetIPhiSize: uint32, jet cluster size in iphi
0017   * trimmedGrid: Flag (bool) to remove three bins in each corner of grid in jet finding
0018   * seedPtThreshold: double, threshold of the seed tower
0019   * pt/eta/philsb : lsb of quantities used in firmware implementation
0020   * puSubtraction: bool, runs chunky doughnut pile-up subtraction, 9x9 jet only
0021   * eta/phiRegionEdges: Boundaries of the input (PF) regions
0022   * maxInputsPerRegion: Truncate number of inputes per input (PF) region
0023   * sin/cosPhi: Value of sin/cos phi in the middle of each bin of the grid.
0024   * met{HF}AbsETaCut: Eta selection of input candidates for calculation of MET
0025   * outputCollectionName: string, tag for the output collection
0026   * vetoZeroPt: bool, controls whether jets with 0 pt should be save. 
0027     It matters if PU is ON, as you can get negative or zero pt jets after it.
0028   * inputCollectionTag: inputtag, collection of reco::candidates used as input to the algo
0029 
0030 */
0031 //
0032 // Original Simone Bologna
0033 // Created: Mon Jul 02 2018
0034 // Modified 2020 Emyr Clement
0035 //
0036 
0037 #include "FWCore/Framework/interface/Frameworkfwd.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039 #include "FWCore/Framework/interface/one/EDProducer.h"
0040 #include "DataFormats/JetReco/interface/CaloJet.h"
0041 #include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
0042 #include "DataFormats/L1TParticleFlow/interface/PFCluster.h"
0043 #include "DataFormats/L1Trigger/interface/L1Candidate.h"
0044 #include "DataFormats/Common/interface/View.h"
0045 #include "DataFormats/Candidate/interface/Candidate.h"
0046 #include "FWCore/Framework/interface/MakerMacros.h"
0047 #include "FWCore/Framework/interface/Event.h"
0048 #include "DataFormats/Math/interface/LorentzVector.h"
0049 #include "FWCore/ServiceRegistry/interface/Service.h"
0050 #include "DataFormats/L1Trigger/interface/EtSum.h"
0051 
0052 #include "TH2F.h"
0053 
0054 #include <cmath>
0055 
0056 #include <algorithm>
0057 constexpr int x_scroll_min = -4;
0058 constexpr int x_scroll_max = 4;
0059 constexpr int y_scroll_min = 0;
0060 constexpr int y_scroll_max = 3;
0061 
0062 class Phase1L1TJetProducer : public edm::one::EDProducer<> {
0063 public:
0064   explicit Phase1L1TJetProducer(const edm::ParameterSet&);
0065   ~Phase1L1TJetProducer() override;
0066 
0067   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0068 
0069 private:
0070   void produce(edm::Event&, const edm::EventSetup&) override;
0071 
0072   /// Finds the seeds in the caloGrid, seeds are saved in a vector that contain the index in the TH2F of each seed
0073   std::vector<std::tuple<int, int>> findSeeds(float seedThreshold) const;
0074 
0075   // Calculates the pt sum of the bins around the seeds
0076   // And applies some PU mitigation
0077   // Warning - not used for some time, so user beware
0078   std::vector<reco::CaloJet> buildJetsFromSeedsWithPUSubtraction(const std::vector<std::tuple<int, int>>& seeds,
0079                                                                  bool killZeroPt) const;
0080 
0081   // Calculates the pt sum of the bins around the seeds
0082   std::vector<reco::CaloJet> buildJetsFromSeeds(const std::vector<std::tuple<int, int>>& seeds) const;
0083 
0084   // Used by buildJetsFromSeedsWithPUSubtraction
0085   // Implementation of pileup mitigation
0086   // Warning - not used for some time, so user beware
0087   void subtract9x9Pileup(reco::CaloJet& jet) const;
0088 
0089   /// Get the energy of a certain tower while correctly handling phi periodicity in case of overflow
0090   float getTowerEnergy(int iEta, int iPhi) const;
0091 
0092   // Implementation of pt sum of bins around one seed
0093   reco::CaloJet buildJetFromSeed(const std::tuple<int, int>& seed) const;
0094 
0095   // <3 handy method to fill the calogrid with whatever type
0096   template <class Container>
0097   void fillCaloGrid(TH2F& caloGrid, const Container& triggerPrimitives, const unsigned int regionIndex);
0098 
0099   // Digitise the eta and phi coordinates of input candidates
0100   // This converts the quantities to integers to reduce precision
0101   // And takes account of bin edge effects i.e. makes sure the
0102   // candidate ends up in the correct (i.e. same behaviour as the firmware) bin of caloGrid_
0103   std::pair<float, float> getCandidateDigiEtaPhi(const float eta,
0104                                                  const float phi,
0105                                                  const unsigned int regionIndex) const;
0106 
0107   // Sorts the input candidates into the PF regions they arrive in
0108   // Truncates the inputs.  Takes the first N candidates as they are provided, without any sorting (this may be needed in the future and/or provided in this way from emulation of layer 1)
0109   template <class Handle>
0110   std::vector<std::vector<reco::CandidatePtr>> prepareInputsIntoRegions(const Handle triggerPrimitives);
0111 
0112   // Converts phi and eta (PF) region indices to a single index
0113   unsigned int getRegionIndex(const unsigned int phiRegion, const unsigned int etaRegion) const;
0114   // From the single index, calculated by getRegionIndex, provides the lower eta and phi boundaries of the input (PF) region index
0115   std::pair<double, double> regionEtaPhiLowEdges(const unsigned int regionIndex) const;
0116   // From the single index, calculated by getRegionIndex, provides the upper eta and phi boundaries of the input (PF) region index
0117   std::pair<double, double> regionEtaPhiUpEdges(const unsigned int regionIndex) const;
0118 
0119   // computes MET
0120   // Takes grid used by jet finder and projects to 1D histogram of phi, bin contents are total pt in that phi bin
0121   // the phi bin index is used to retrieve the sin-cos value from the LUT emulator
0122   // the pt of the input is multiplied by that sin cos value to obtain px and py that is added to the total event px & py
0123   // after all the inputs have been processed we compute the total pt of the event, and set that as MET
0124   l1t::EtSum computeMET(const double etaCut, l1t::EtSum::EtSumType sumType) const;
0125 
0126   // Determine if this tower should be trimmed or not
0127   // Used only when trimmedGrid_ option is set to true
0128   // Trim means removing 3 towers in each corner of the square grid
0129   // giving a cross shaped grid, which is a bit more circular in shape than a square
0130   bool trimTower(const int etaIndex, const int phiIndex) const;
0131 
0132   edm::EDGetTokenT<edm::View<reco::Candidate>> inputCollectionTag_;
0133   // histogram containing our clustered inputs
0134   std::unique_ptr<TH2F> caloGrid_;
0135 
0136   std::vector<double> etaBinning_;
0137   size_t nBinsEta_;
0138   unsigned int nBinsPhi_;
0139   double phiLow_;
0140   double phiUp_;
0141   unsigned int jetIEtaSize_;
0142   unsigned int jetIPhiSize_;
0143   bool trimmedGrid_;
0144   double seedPtThreshold_;
0145   double ptlsb_;
0146   double philsb_;
0147   double etalsb_;
0148   bool puSubtraction_;
0149   bool vetoZeroPt_;
0150   // Eta and phi edges of input PF regions
0151   std::vector<double> etaRegionEdges_;
0152   std::vector<double> phiRegionEdges_;
0153   // Maximum number of candidates per input PF region
0154   unsigned int maxInputsPerRegion_;
0155   // LUT for sin and cos phi, to match that used in firmware
0156   std::vector<double> sinPhi_;
0157   std::vector<double> cosPhi_;
0158   // input eta cut for met calculation
0159   double metAbsEtaCut_;
0160   // input eta cut for metHF calculation
0161   double metHFAbsEtaCut_;
0162   std::string outputCollectionName_;
0163 };
0164 
0165 Phase1L1TJetProducer::Phase1L1TJetProducer(const edm::ParameterSet& iConfig)
0166     :  // getting configuration settings
0167       inputCollectionTag_{
0168           consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("inputCollectionTag"))},
0169       etaBinning_(iConfig.getParameter<std::vector<double>>("etaBinning")),
0170       nBinsEta_(etaBinning_.size() - 1),
0171       nBinsPhi_(iConfig.getParameter<unsigned int>("nBinsPhi")),
0172       phiLow_(iConfig.getParameter<double>("phiLow")),
0173       phiUp_(iConfig.getParameter<double>("phiUp")),
0174       jetIEtaSize_(iConfig.getParameter<unsigned int>("jetIEtaSize")),
0175       jetIPhiSize_(iConfig.getParameter<unsigned int>("jetIPhiSize")),
0176       trimmedGrid_(iConfig.getParameter<bool>("trimmedGrid")),
0177       seedPtThreshold_(iConfig.getParameter<double>("seedPtThreshold")),
0178       ptlsb_(iConfig.getParameter<double>("ptlsb")),
0179       philsb_(iConfig.getParameter<double>("philsb")),
0180       etalsb_(iConfig.getParameter<double>("etalsb")),
0181       puSubtraction_(iConfig.getParameter<bool>("puSubtraction")),
0182       vetoZeroPt_(iConfig.getParameter<bool>("vetoZeroPt")),
0183       etaRegionEdges_(iConfig.getParameter<std::vector<double>>("etaRegions")),
0184       phiRegionEdges_(iConfig.getParameter<std::vector<double>>("phiRegions")),
0185       maxInputsPerRegion_(iConfig.getParameter<unsigned int>("maxInputsPerRegion")),
0186       sinPhi_(iConfig.getParameter<std::vector<double>>("sinPhi")),
0187       cosPhi_(iConfig.getParameter<std::vector<double>>("cosPhi")),
0188       metAbsEtaCut_(iConfig.getParameter<double>("metAbsEtaCut")),
0189       metHFAbsEtaCut_(iConfig.getParameter<double>("metHFAbsEtaCut")),
0190       outputCollectionName_(iConfig.getParameter<std::string>("outputCollectionName")) {
0191   caloGrid_ =
0192       std::make_unique<TH2F>("caloGrid", "Calorimeter grid", nBinsEta_, etaBinning_.data(), nBinsPhi_, phiLow_, phiUp_);
0193   caloGrid_->GetXaxis()->SetTitle("#eta");
0194   caloGrid_->GetYaxis()->SetTitle("#phi");
0195   produces<std::vector<reco::CaloJet>>(outputCollectionName_).setBranchAlias(outputCollectionName_);
0196   produces<std::vector<l1t::EtSum>>(outputCollectionName_ + "MET").setBranchAlias(outputCollectionName_ + "MET");
0197 }
0198 
0199 Phase1L1TJetProducer::~Phase1L1TJetProducer() {}
0200 
0201 float Phase1L1TJetProducer::getTowerEnergy(int iEta, int iPhi) const {
0202   // We return the pt of a certain bin in the calo grid, taking account of the phi periodicity when overflowing (e.g. phi > phiSize), and returning 0 for the eta out of bounds
0203 
0204   int nBinsEta = caloGrid_->GetNbinsX();
0205   int nBinsPhi = caloGrid_->GetNbinsY();
0206   while (iPhi < 1) {
0207     iPhi += nBinsPhi;
0208   }
0209   while (iPhi > nBinsPhi) {
0210     iPhi -= nBinsPhi;
0211   }
0212   if (iEta < 1) {
0213     return 0;
0214   }
0215   if (iEta > nBinsEta) {
0216     return 0;
0217   }
0218   return caloGrid_->GetBinContent(iEta, iPhi);
0219 }
0220 
0221 void Phase1L1TJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0222   edm::Handle<edm::View<reco::Candidate>> inputCollectionHandle;
0223   iEvent.getByToken(inputCollectionTag_, inputCollectionHandle);
0224 
0225   // sort inputs into PF regions
0226   std::vector<std::vector<reco::CandidatePtr>> inputsInRegions = prepareInputsIntoRegions<>(inputCollectionHandle);
0227 
0228   // histogramming the data
0229   caloGrid_->Reset();
0230   for (unsigned int iInputRegion = 0; iInputRegion < inputsInRegions.size(); ++iInputRegion) {
0231     fillCaloGrid<>(*(caloGrid_), inputsInRegions[iInputRegion], iInputRegion);
0232   }
0233 
0234   // find the seeds
0235   const auto& seedsVector = findSeeds(seedPtThreshold_);  // seedPtThreshold = 5
0236   // build jets from the seeds
0237   auto l1jetVector =
0238       puSubtraction_ ? buildJetsFromSeedsWithPUSubtraction(seedsVector, vetoZeroPt_) : buildJetsFromSeeds(seedsVector);
0239 
0240   // sort by pt
0241   std::sort(l1jetVector.begin(), l1jetVector.end(), [](const reco::CaloJet& jet1, const reco::CaloJet& jet2) {
0242     return jet1.pt() > jet2.pt();
0243   });
0244 
0245   auto l1jetVectorPtr = std::make_unique<std::vector<reco::CaloJet>>(l1jetVector);
0246   iEvent.put(std::move(l1jetVectorPtr), outputCollectionName_);
0247 
0248   // calculate METs
0249   l1t::EtSum lMET = computeMET(metAbsEtaCut_, l1t::EtSum::EtSumType::kMissingEt);
0250   l1t::EtSum lMETHF = computeMET(metHFAbsEtaCut_, l1t::EtSum::EtSumType::kMissingEtHF);
0251   std::unique_ptr<std::vector<l1t::EtSum>> lSumVectorPtr(new std::vector<l1t::EtSum>(0));
0252   lSumVectorPtr->push_back(lMET);
0253   lSumVectorPtr->push_back(lMETHF);
0254   iEvent.put(std::move(lSumVectorPtr), this->outputCollectionName_ + "MET");
0255 
0256   return;
0257 }
0258 
0259 void Phase1L1TJetProducer::subtract9x9Pileup(reco::CaloJet& jet) const {
0260   // these variables host the total pt in each sideband and the total pileup contribution
0261   float topBandPt = 0;
0262   float leftBandPt = 0;
0263   float rightBandPt = 0;
0264   float bottomBandPt = 0;
0265   float pileUpEnergy;
0266 
0267   // hold the jet's x-y (and z, as I have to use it, even if 2D) location in the histo
0268   int xCenter, yCenter, zCenter;
0269   // Retrieving histo-coords for seed
0270   caloGrid_->GetBinXYZ(caloGrid_->FindFixBin(jet.eta(), jet.phi()), xCenter, yCenter, zCenter);
0271 
0272   // Computing pileup
0273   for (int x = x_scroll_min; x <= x_scroll_max; x++) {
0274     for (int y = y_scroll_min; y < y_scroll_max; y++) {
0275       // top band, I go up 5 squares to reach the bottom of the top band
0276       // +x scrolls from left to right, +y scrolls up
0277       topBandPt += getTowerEnergy(xCenter + x, yCenter + (5 + y));
0278       // left band, I go left 5 squares (-5) to reach the bottom of the top band
0279       // +x scrolls from bottom to top, +y scrolls left
0280       leftBandPt += getTowerEnergy(xCenter - (5 + y), yCenter + x);
0281       // right band, I go right 5 squares (+5) to reach the bottom of the top band
0282       // +x scrolls from bottom to top, +y scrolls right
0283       rightBandPt += getTowerEnergy(xCenter + (5 + y), yCenter + x);
0284       // right band, I go right 5 squares (+5) to reach the bottom of the top band
0285       // +x scrolls from bottom to top, +y scrolls right
0286       bottomBandPt += getTowerEnergy(xCenter + x, yCenter - (5 + y));
0287     }
0288   }
0289   // adding bands and removing the maximum band (equivalent to adding the three minimum bands)
0290   pileUpEnergy = topBandPt + leftBandPt + rightBandPt + bottomBandPt -
0291                  std::max(topBandPt, std::max(leftBandPt, std::max(rightBandPt, bottomBandPt)));
0292 
0293   //preparing the new 4-momentum vector
0294   reco::Candidate::PolarLorentzVector ptVector;
0295   // removing pu contribution
0296   float ptAfterPUSubtraction = jet.pt() - pileUpEnergy;
0297   ptVector.SetPt((ptAfterPUSubtraction > 0) ? ptAfterPUSubtraction : 0);
0298   ptVector.SetEta(jet.eta());
0299   ptVector.SetPhi(jet.phi());
0300   //updating the jet
0301   jet.setP4(ptVector);
0302   jet.setPileup(pileUpEnergy);
0303   return;
0304 }
0305 
0306 std::vector<std::tuple<int, int>> Phase1L1TJetProducer::findSeeds(float seedThreshold) const {
0307   int nBinsX = caloGrid_->GetNbinsX();
0308   int nBinsY = caloGrid_->GetNbinsY();
0309 
0310   std::vector<std::tuple<int, int>> seeds;
0311 
0312   int etaHalfSize = (int)jetIEtaSize_ / 2;
0313   int phiHalfSize = (int)jetIPhiSize_ / 2;
0314 
0315   // for each point of the grid check if it is a local maximum
0316   // to do so I take a point, and look if is greater than the points around it (in the 9x9 neighborhood)
0317   // to prevent mutual exclusion, I check greater or equal for points above and right to the one I am considering (including the top-left point)
0318   // to prevent mutual exclusion, I check greater for points below and left to the one I am considering (including the bottom-right point)
0319 
0320   for (int iPhi = 1; iPhi <= nBinsY; iPhi++) {
0321     for (int iEta = 1; iEta <= nBinsX; iEta++) {
0322       float centralPt = caloGrid_->GetBinContent(iEta, iPhi);
0323       if (centralPt < seedThreshold)
0324         continue;
0325       bool isLocalMaximum = true;
0326 
0327       // Scanning through the grid centered on the seed
0328       for (int etaIndex = -etaHalfSize; etaIndex <= etaHalfSize; etaIndex++) {
0329         for (int phiIndex = -phiHalfSize; phiIndex <= phiHalfSize; phiIndex++) {
0330           if (trimmedGrid_) {
0331             if (trimTower(etaIndex, phiIndex))
0332               continue;
0333           }
0334 
0335           if ((etaIndex == 0) && (phiIndex == 0))
0336             continue;
0337           if (phiIndex > 0) {
0338             if (phiIndex > -etaIndex) {
0339               isLocalMaximum = ((isLocalMaximum) && (centralPt > getTowerEnergy(iEta + etaIndex, iPhi + phiIndex)));
0340             } else {
0341               isLocalMaximum = ((isLocalMaximum) && (centralPt >= getTowerEnergy(iEta + etaIndex, iPhi + phiIndex)));
0342             }
0343           } else {
0344             if (phiIndex >= -etaIndex) {
0345               isLocalMaximum = ((isLocalMaximum) && (centralPt > getTowerEnergy(iEta + etaIndex, iPhi + phiIndex)));
0346             } else {
0347               isLocalMaximum = ((isLocalMaximum) && (centralPt >= getTowerEnergy(iEta + etaIndex, iPhi + phiIndex)));
0348             }
0349           }
0350         }
0351       }
0352       if (isLocalMaximum) {
0353         seeds.emplace_back(iEta, iPhi);
0354       }
0355     }
0356   }
0357 
0358   return seeds;
0359 }
0360 
0361 reco::CaloJet Phase1L1TJetProducer::buildJetFromSeed(const std::tuple<int, int>& seed) const {
0362   int iEta = std::get<0>(seed);
0363   int iPhi = std::get<1>(seed);
0364 
0365   int etaHalfSize = (int)jetIEtaSize_ / 2;
0366   int phiHalfSize = (int)jetIPhiSize_ / 2;
0367 
0368   float ptSum = 0;
0369   // Scanning through the grid centered on the seed
0370   for (int etaIndex = -etaHalfSize; etaIndex <= etaHalfSize; etaIndex++) {
0371     for (int phiIndex = -phiHalfSize; phiIndex <= phiHalfSize; phiIndex++) {
0372       if (trimmedGrid_) {
0373         if (trimTower(etaIndex, phiIndex))
0374           continue;
0375       }
0376       ptSum += getTowerEnergy(iEta + etaIndex, iPhi + phiIndex);
0377     }
0378   }
0379 
0380   // Creating a jet with eta phi centered on the seed and momentum equal to the sum of the pt of the components
0381   reco::Candidate::PolarLorentzVector ptVector;
0382   ptVector.SetPt(ptSum);
0383   ptVector.SetEta(caloGrid_->GetXaxis()->GetBinCenter(iEta));
0384   ptVector.SetPhi(caloGrid_->GetYaxis()->GetBinCenter(iPhi));
0385   reco::CaloJet jet;
0386   jet.setP4(ptVector);
0387   return jet;
0388 }
0389 
0390 std::vector<reco::CaloJet> Phase1L1TJetProducer::buildJetsFromSeedsWithPUSubtraction(
0391     const std::vector<std::tuple<int, int>>& seeds, bool killZeroPt) const {
0392   // For each seed take a grid centered on the seed of the size specified by the user
0393   // Sum the pf in the grid, that will be the pt of the l1t jet. Eta and phi of the jet is taken from the seed.
0394   std::vector<reco::CaloJet> jets;
0395   for (const auto& seed : seeds) {
0396     reco::CaloJet jet = buildJetFromSeed(seed);
0397     subtract9x9Pileup(jet);
0398     //killing jets with 0 pt
0399     if ((vetoZeroPt_) && (jet.pt() <= 0))
0400       continue;
0401     jets.push_back(jet);
0402   }
0403   return jets;
0404 }
0405 
0406 std::vector<reco::CaloJet> Phase1L1TJetProducer::buildJetsFromSeeds(
0407     const std::vector<std::tuple<int, int>>& seeds) const {
0408   // For each seed take a grid centered on the seed of the size specified by the user
0409   // Sum the pf in the grid, that will be the pt of the l1t jet. Eta and phi of the jet is taken from the seed.
0410   std::vector<reco::CaloJet> jets;
0411   for (const auto& seed : seeds) {
0412     reco::CaloJet jet = buildJetFromSeed(seed);
0413     jets.push_back(jet);
0414   }
0415   return jets;
0416 }
0417 
0418 template <class Container>
0419 void Phase1L1TJetProducer::fillCaloGrid(TH2F& caloGrid,
0420                                         const Container& triggerPrimitives,
0421                                         const unsigned int regionIndex) {
0422   //Filling the calo grid with the primitives
0423   for (const auto& primitiveIterator : triggerPrimitives) {
0424     // Get digitised (floating point with reduced precision) eta and phi
0425     std::pair<float, float> digi_EtaPhi =
0426         getCandidateDigiEtaPhi(primitiveIterator->eta(), primitiveIterator->phi(), regionIndex);
0427 
0428     caloGrid.Fill(static_cast<float>(digi_EtaPhi.first),
0429                   static_cast<float>(digi_EtaPhi.second),
0430                   static_cast<float>(primitiveIterator->pt()));
0431   }
0432 }
0433 
0434 std::pair<float, float> Phase1L1TJetProducer::getCandidateDigiEtaPhi(const float eta,
0435                                                                      const float phi,
0436                                                                      const unsigned int regionIndex) const {
0437   std::pair<double, double> regionLowEdges = regionEtaPhiLowEdges(regionIndex);
0438 
0439   int digitisedEta = floor((eta - regionLowEdges.second) / etalsb_);
0440   int digitisedPhi = floor((phi - regionLowEdges.first) / philsb_);
0441 
0442   // If eta or phi is on a bin edge
0443   // Put in bin above, to match behaviour of HLS
0444   // Unless it's on the last bin of this pf region
0445   // Then it is placed in the last bin, not the overflow
0446   TAxis* etaAxis = caloGrid_->GetXaxis();
0447   std::pair<double, double> regionUpEdges = regionEtaPhiUpEdges(regionIndex);
0448   int digiEtaEdgeLastBinUp = floor((regionUpEdges.second - regionLowEdges.second) / etalsb_);
0449   // If the digi eta is outside the last bin of this pf region
0450   // Set the digitised quantity so it would be in the last bin
0451   // These cases could be avoided by sorting input candidates based on digitised eta/phi
0452   if (digitisedEta >= digiEtaEdgeLastBinUp) {
0453     digitisedEta = digiEtaEdgeLastBinUp - 1;
0454   } else {
0455     for (int i = 0; i < etaAxis->GetNbins(); ++i) {
0456       if (etaAxis->GetBinUpEdge(i) < regionLowEdges.second)
0457         continue;
0458       int digiEdgeBinUp = floor((etaAxis->GetBinUpEdge(i) - regionLowEdges.second) / etalsb_);
0459       if (digiEdgeBinUp == digitisedEta) {
0460         digitisedEta += 1;
0461       }
0462     }
0463   }
0464 
0465   // Similar for phi
0466   TAxis* phiAxis = caloGrid_->GetYaxis();
0467   int digiPhiEdgeLastBinUp = floor((regionUpEdges.first - regionLowEdges.first) / philsb_);
0468   if (digitisedPhi >= digiPhiEdgeLastBinUp) {
0469     digitisedPhi = digiPhiEdgeLastBinUp - 1;
0470   } else {
0471     for (int i = 0; i < phiAxis->GetNbins(); ++i) {
0472       if (phiAxis->GetBinUpEdge(i) < regionLowEdges.first)
0473         continue;
0474       int digiEdgeBinUp = floor((phiAxis->GetBinUpEdge(i) - regionLowEdges.first) / philsb_);
0475       if (digiEdgeBinUp == digitisedPhi) {
0476         digitisedPhi += 1;
0477       }
0478     }
0479   }
0480 
0481   // Convert digitised eta and phi back to floating point quantities with reduced precision
0482   float floatDigitisedEta = (digitisedEta + 0.5) * etalsb_ + regionLowEdges.second;
0483   float floatDigitisedPhi = (digitisedPhi + 0.5) * philsb_ + regionLowEdges.first;
0484 
0485   return std::pair<float, float>{floatDigitisedEta, floatDigitisedPhi};
0486 }
0487 
0488 void Phase1L1TJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0489   edm::ParameterSetDescription desc;
0490   desc.add<edm::InputTag>("inputCollectionTag", edm::InputTag("l1pfCandidates", "Puppi"));
0491   desc.add<std::vector<double>>("etaBinning");
0492   desc.add<unsigned int>("nBinsPhi", 72);
0493   desc.add<double>("phiLow", -M_PI);
0494   desc.add<double>("phiUp", M_PI);
0495   desc.add<unsigned int>("jetIEtaSize", 7);
0496   desc.add<unsigned int>("jetIPhiSize", 7);
0497   desc.add<bool>("trimmedGrid", false);
0498   desc.add<double>("seedPtThreshold", 5);
0499   desc.add<double>("ptlsb", 0.25), desc.add<double>("philsb", 0.0043633231), desc.add<double>("etalsb", 0.0043633231),
0500       desc.add<bool>("puSubtraction", false);
0501   desc.add<string>("outputCollectionName", "UncalibratedPhase1L1TJetFromPfCandidates");
0502   desc.add<bool>("vetoZeroPt", true);
0503   desc.add<std::vector<double>>("etaRegions");
0504   desc.add<std::vector<double>>("phiRegions");
0505   desc.add<unsigned int>("maxInputsPerRegion", 18);
0506   desc.add<std::vector<double>>("sinPhi");
0507   desc.add<std::vector<double>>("cosPhi");
0508   desc.add<double>("metAbsEtaCut", 3);
0509   desc.add<double>("metHFAbsEtaCut", 5);
0510   descriptions.add("Phase1L1TJetProducer", desc);
0511 }
0512 
0513 template <class Handle>
0514 std::vector<std::vector<edm::Ptr<reco::Candidate>>> Phase1L1TJetProducer::prepareInputsIntoRegions(
0515     const Handle triggerPrimitives) {
0516   std::vector<std::vector<reco::CandidatePtr>> inputsInRegions{etaRegionEdges_.size() * (phiRegionEdges_.size() - 1)};
0517 
0518   for (unsigned int i = 0; i < triggerPrimitives->size(); ++i) {
0519     reco::CandidatePtr tp(triggerPrimitives, i);
0520 
0521     if (tp->phi() < phiRegionEdges_.front() || tp->phi() >= phiRegionEdges_.back() ||
0522         tp->eta() < etaRegionEdges_.front() || tp->eta() >= etaRegionEdges_.back())
0523       continue;
0524 
0525     // Which phi region does this tp belong to
0526     auto it_phi = phiRegionEdges_.begin();
0527     it_phi = std::upper_bound(phiRegionEdges_.begin(), phiRegionEdges_.end(), tp->phi()) - 1;
0528 
0529     // Which eta region does this tp belong to
0530     auto it_eta = etaRegionEdges_.begin();
0531     it_eta = std::upper_bound(etaRegionEdges_.begin(), etaRegionEdges_.end(), tp->eta()) - 1;
0532 
0533     if (it_phi != phiRegionEdges_.end() && it_eta != etaRegionEdges_.end()) {
0534       auto phiRegion = it_phi - phiRegionEdges_.begin();
0535       auto etaRegion = it_eta - etaRegionEdges_.begin();
0536       inputsInRegions[getRegionIndex(phiRegion, etaRegion)].emplace_back(tp);
0537     }
0538   }
0539 
0540   // Truncate number of inputs in each pf region
0541   for (auto& inputs : inputsInRegions) {
0542     if (inputs.size() > maxInputsPerRegion_) {
0543       inputs.resize(maxInputsPerRegion_);
0544     }
0545   }
0546 
0547   return inputsInRegions;
0548 }
0549 
0550 unsigned int Phase1L1TJetProducer::getRegionIndex(const unsigned int phiRegion, const unsigned int etaRegion) const {
0551   return etaRegion * (phiRegionEdges_.size() - 1) + phiRegion;
0552 }
0553 
0554 std::pair<double, double> Phase1L1TJetProducer::regionEtaPhiLowEdges(const unsigned int regionIndex) const {
0555   unsigned int phiRegion = regionIndex % (phiRegionEdges_.size() - 1);
0556   unsigned int etaRegion = (regionIndex - phiRegion) / (phiRegionEdges_.size() - 1);
0557   return std::pair<double, double>{phiRegionEdges_.at(phiRegion), etaRegionEdges_.at(etaRegion)};
0558 }
0559 
0560 std::pair<double, double> Phase1L1TJetProducer::regionEtaPhiUpEdges(const unsigned int regionIndex) const {
0561   unsigned int phiRegion = regionIndex % (phiRegionEdges_.size() - 1);
0562   unsigned int etaRegion = (regionIndex - phiRegion) / (phiRegionEdges_.size() - 1);
0563   if (phiRegion == phiRegionEdges_.size() - 1) {
0564     return std::pair<double, double>{phiRegionEdges_.at(phiRegion), etaRegionEdges_.at(etaRegion + 1)};
0565   } else if (etaRegion == etaRegionEdges_.size() - 1) {
0566     return std::pair<double, double>{phiRegionEdges_.at(phiRegion + 1), etaRegionEdges_.at(etaRegion)};
0567   }
0568 
0569   return std::pair<double, double>{phiRegionEdges_.at(phiRegion + 1), etaRegionEdges_.at(etaRegion + 1)};
0570 }
0571 
0572 l1t::EtSum Phase1L1TJetProducer::computeMET(const double etaCut, l1t::EtSum::EtSumType sumType) const {
0573   const auto lowEtaBin = caloGrid_->GetXaxis()->FindBin(-1.0 * etaCut);
0574   const auto highEtaBin = caloGrid_->GetXaxis()->FindBin(etaCut) - 1;
0575   const auto phiProjection = caloGrid_->ProjectionY("temp", lowEtaBin, highEtaBin);
0576 
0577   // Use digitised quantities when summing to improve agreement with firmware
0578   int totalDigiPx{0};
0579   int totalDigiPy{0};
0580 
0581   for (int i = 1; i < phiProjection->GetNbinsX() + 1; ++i) {
0582     double pt = phiProjection->GetBinContent(i);
0583     totalDigiPx += trunc(floor(pt / ptlsb_) * cosPhi_[i - 1]);
0584     totalDigiPy += trunc(floor(pt / ptlsb_) * sinPhi_[i - 1]);
0585   }
0586 
0587   double lMET = floor(sqrt(totalDigiPx * totalDigiPx + totalDigiPy * totalDigiPy)) * ptlsb_;
0588 
0589   math::PtEtaPhiMLorentzVector lMETVector(lMET, 0, acos(totalDigiPx / (lMET / ptlsb_)), 0);
0590   l1t::EtSum lMETSum(lMETVector, sumType, 0, 0, 0, 0);
0591   return lMETSum;
0592 }
0593 
0594 bool Phase1L1TJetProducer::trimTower(const int etaIndex, const int phiIndex) const {
0595   int etaHalfSize = jetIEtaSize_ / 2;
0596   int phiHalfSize = jetIPhiSize_ / 2;
0597 
0598   if (etaIndex == -etaHalfSize || etaIndex == etaHalfSize) {
0599     if (phiIndex <= -phiHalfSize + 1 || phiIndex >= phiHalfSize - 1) {
0600       return true;
0601     }
0602   } else if (etaIndex == -etaHalfSize + 1 || etaIndex == etaHalfSize - 1) {
0603     if (phiIndex == -phiHalfSize || phiIndex == phiHalfSize) {
0604       return true;
0605     }
0606   }
0607 
0608   return false;
0609 }
0610 DEFINE_FWK_MODULE(Phase1L1TJetProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!