Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-16 00:33:18

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoHI/HiJetAlgos/plugins/HiPuRhoProducer
0004 // Class:      HiPuRhoProducer
0005 //
0006 /**\class HiPuRhoProducer HiPuRhoProducer.cc RecoHI/HiJetAlgos/plugins/HiPuRhoProducer.cc
0007  Description: Producer to dump Pu-jet style rho into event content
0008  Implementation:
0009  Just see MultipleAlgoIterator - re-implenting for use in CS jet with sigma subtraction and zeroing
0010 */
0011 //
0012 // Original Author:  Chris McGinn - in Style of HiFJRhoProducer
0013 //         Created:  Mon, 29 May 2017
0014 //
0015 //
0016 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0017 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0018 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0019 #include "DataFormats/Candidate/interface/Candidate.h"
0020 #include "DataFormats/Common/interface/Handle.h"
0021 #include "DataFormats/Common/interface/View.h"
0022 #include "DataFormats/DetId/interface/DetId.h"
0023 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0024 #include "DataFormats/Math/interface/LorentzVector.h"
0025 #include "DataFormats/Math/interface/deltaR.h"
0026 #include "FWCore/Framework/interface/ESHandle.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
0034 #include "FWCore/Utilities/interface/Exception.h"
0035 #include "FWCore/Utilities/interface/isFinite.h"
0036 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0037 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0038 #include "DataFormats/Common/interface/Ptr.h"
0039 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0040 #include "FWCore/Framework/interface/stream/EDProducer.h"
0041 #include "FWCore/Framework/interface/EventSetup.h"
0042 #include "FWCore/Utilities/interface/EDGetToken.h"
0043 #include "FWCore/Utilities/interface/InputTag.h"
0044 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0045 #include "RecoHI/HiJetAlgos/plugins/HITowerHelper.h"
0046 
0047 #include "fastjet/ClusterSequence.hh"
0048 #include "fastjet/JetDefinition.hh"
0049 #include "fastjet/PseudoJet.hh"
0050 
0051 #include <cmath>
0052 #include <cstdint>
0053 #include <cstdlib>
0054 #include <algorithm>
0055 #include <limits>
0056 #include <memory>
0057 #include <sstream>
0058 #include <string>
0059 #include <utility>
0060 #include <map>
0061 #include <vector>
0062 
0063 struct EtaPhiTower {
0064   int ieta, iphi;
0065   float eta, phi;
0066 };
0067 
0068 class HiPuRhoProducer : public edm::stream::EDProducer<> {
0069 public:
0070   explicit HiPuRhoProducer(const edm::ParameterSet&);
0071 
0072   using ClusterSequencePtr = std::shared_ptr<fastjet::ClusterSequence>;
0073   using JetDefPtr = std::shared_ptr<fastjet::JetDefinition>;
0074 
0075   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0076   virtual void setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup);
0077   virtual void calculatePedestal(std::vector<fastjet::PseudoJet> const& coll);
0078   virtual void subtractPedestal(std::vector<fastjet::PseudoJet>& coll);
0079   virtual void calculateOrphanInput(std::vector<fastjet::PseudoJet>& orphanInput);
0080   virtual void putRho(edm::Event& iEvent, const edm::EventSetup& iSetup);
0081 
0082 private:
0083   void produce(edm::Event&, const edm::EventSetup&) override;
0084 
0085   // This checks if the tower is anomalous (if a calo tower).
0086   virtual void inputTowers();
0087 
0088   bool postOrphan_ = false;
0089   bool setInitialValue_;
0090 
0091   const bool dropZeroTowers_;
0092   const int medianWindowWidth_;
0093   const double minimumTowersFraction_;
0094   const double nSigmaPU_;  // number of sigma for pileup
0095   const double puPtMin_;
0096   const double radiusPU_;  // pileup radius
0097   const double rParam_;    // the R parameter to use
0098   const double towSigmaCut_;
0099   const edm::EDGetTokenT<CaloTowerCollection> caloTowerToken_;
0100   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0101   const int initialValue_ = -99;
0102   constexpr static int nEtaTow_ = 82;
0103 
0104   std::array<int, nEtaTow_> vngeom_;
0105   std::array<int, nEtaTow_> vntow_;
0106   std::array<float, nEtaTow_> vmean0_;
0107   std::array<float, nEtaTow_> vrms0_;
0108   std::array<float, nEtaTow_> vrho0_;
0109   std::array<float, nEtaTow_> vmean1_;
0110   std::array<float, nEtaTow_> vrms1_;
0111   std::array<float, nEtaTow_> vrho1_;
0112 
0113   std::vector<double> etaEdgeLow_;
0114   std::vector<double> etaEdgeHi_;
0115   std::vector<double> etaEdges_;
0116 
0117   std::vector<double> rho_;
0118   std::vector<double> rhoExtra_;
0119   std::vector<double> rhoM_;
0120   std::vector<int> nTow_;
0121 
0122   std::vector<double> towExcludePt_;
0123   std::vector<double> towExcludePhi_;
0124   std::vector<double> towExcludeEta_;
0125 
0126   std::vector<const CaloTower*> inputs_;
0127   ClusterSequencePtr fjClusterSeq_;  // fastjet cluster sequence
0128   JetDefPtr fjJetDefinition_;        // fastjet jet definition
0129 
0130   std::vector<fastjet::PseudoJet> fjInputs_;          // fastjet inputs
0131   std::vector<fastjet::PseudoJet> fjJets_;            // fastjet jets
0132   std::vector<fastjet::PseudoJet> fjOriginalInputs_;  // to back-up unsubtracted fastjet inputs
0133 
0134   CaloGeometry const* geo_ = nullptr;  // geometry
0135   std::vector<HcalDetId> allgeomid_;   // all det ids in the geometry
0136 
0137   int ietamax_;                                 // maximum eta in geometry
0138   int ietamin_;                                 // minimum eta in geometry
0139   std::map<int, int> ntowersWithJets_;          // number of towers with jets
0140   std::map<int, int> geomtowers_;               // map of geometry towers to det id
0141   std::map<int, double> esigma_;                // energy sigma
0142   std::map<int, double> emean_;                 // energy mean
0143   std::map<int, std::array<double, 4>> eTop4_;  // energy mean
0144 
0145   typedef std::pair<double, double> EtaPhi;
0146   std::vector<EtaPhiTower> towermap_;
0147 };
0148 
0149 HiPuRhoProducer::HiPuRhoProducer(const edm::ParameterSet& iConfig)
0150     : dropZeroTowers_(iConfig.getParameter<bool>("dropZeroTowers")),
0151       medianWindowWidth_(iConfig.getParameter<int>("medianWindowWidth")),
0152       minimumTowersFraction_(iConfig.getParameter<double>("minimumTowersFraction")),
0153       nSigmaPU_(iConfig.getParameter<double>("nSigmaPU")),
0154       puPtMin_(iConfig.getParameter<double>("puPtMin")),
0155       radiusPU_(iConfig.getParameter<double>("radiusPU")),
0156       rParam_(iConfig.getParameter<double>("rParam")),
0157       towSigmaCut_(iConfig.getParameter<double>("towSigmaCut")),
0158       caloTowerToken_(consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("src"))),
0159       caloGeometryToken_(esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})) {
0160   //register your products
0161   produces<std::vector<double>>("mapEtaEdges");
0162   produces<std::vector<double>>("mapToRho");
0163   produces<std::vector<double>>("mapToRhoMedian");
0164   produces<std::vector<double>>("mapToRhoExtra");
0165   produces<std::vector<double>>("mapToRhoM");
0166   produces<std::vector<int>>("mapToNTow");
0167   produces<std::vector<double>>("mapToTowExcludePt");
0168   produces<std::vector<double>>("mapToTowExcludePhi");
0169   produces<std::vector<double>>("mapToTowExcludeEta");
0170 }
0171 
0172 // ------------ method called to produce the data  ------------
0173 void HiPuRhoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0174   setupGeometryMap(iEvent, iSetup);
0175 
0176   for (int i = ietamin_; i < ietamax_ + 1; i++) {
0177     ntowersWithJets_[i] = 0;
0178   }
0179 
0180   auto const& inputView = iEvent.get(caloTowerToken_);
0181   inputs_.reserve(inputView.size());
0182   for (auto const& input : inputView)
0183     inputs_.push_back(&input);
0184 
0185   fjInputs_.reserve(inputs_.size());
0186   inputTowers();
0187   fjOriginalInputs_ = fjInputs_;
0188   setInitialValue_ = true;
0189   calculatePedestal(fjInputs_);
0190   subtractPedestal(fjInputs_);
0191 
0192   fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, rParam_);
0193   fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, *fjJetDefinition_);
0194   fjJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(puPtMin_));
0195 
0196   etaEdgeLow_.clear();
0197   etaEdgeHi_.clear();
0198   etaEdges_.clear();
0199 
0200   rho_.clear();
0201   rhoExtra_.clear();
0202   rhoM_.clear();
0203   nTow_.clear();
0204 
0205   towExcludePt_.clear();
0206   towExcludePhi_.clear();
0207   towExcludeEta_.clear();
0208 
0209   setInitialValue_ = false;
0210   std::vector<fastjet::PseudoJet> orphanInput;
0211   calculateOrphanInput(orphanInput);
0212   calculatePedestal(orphanInput);
0213   putRho(iEvent, iSetup);
0214 
0215   inputs_.clear();
0216   fjInputs_.clear();
0217   fjJets_.clear();
0218 }
0219 
0220 void HiPuRhoProducer::inputTowers() {
0221   int index = -1;
0222   for (auto const& input : inputs_) {
0223     index++;
0224 
0225     if (edm::isNotFinite(input->pt()))
0226       continue;
0227     if (input->pt() < 100 * std::numeric_limits<double>::epsilon())
0228       continue;
0229 
0230     fjInputs_.push_back(fastjet::PseudoJet(input->px(), input->py(), input->pz(), input->energy()));
0231     fjInputs_.back().set_user_index(index);
0232   }
0233 }
0234 
0235 void HiPuRhoProducer::setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0236   LogDebug("PileUpSubtractor") << "The subtractor setting up geometry...\n";
0237   const auto& pG = iSetup.getData(caloGeometryToken_);
0238   geo_ = &pG;
0239   std::vector<DetId> alldid = geo_->getValidDetIds();
0240   int ietaold = -10000;
0241   ietamax_ = -10000;
0242   ietamin_ = 10000;
0243   towermap_.clear();
0244 
0245   for (auto const& did : alldid) {
0246     if (did.det() == DetId::Hcal) {
0247       HcalDetId hid = HcalDetId(did);
0248       allgeomid_.push_back(did);
0249       towermap_.push_back({hid.ieta(), hid.iphi(), geo_->getPosition(did).eta(), geo_->getPosition(did).phi()});
0250       if (hid.ieta() != ietaold) {
0251         ietaold = hid.ieta();
0252         geomtowers_[hid.ieta()] = 1;
0253         if (hid.ieta() > ietamax_)
0254           ietamax_ = hid.ieta();
0255         if (hid.ieta() < ietamin_)
0256           ietamin_ = hid.ieta();
0257       } else {
0258         geomtowers_[hid.ieta()]++;
0259       }
0260     }
0261   }
0262 
0263   for (auto const& gt : geomtowers_) {
0264     int it = gt.first;
0265     int vi = it - 1;
0266 
0267     if (it < 0)
0268       vi = nEtaTow_ + it;
0269 
0270     if (vi < nEtaTow_)
0271       vngeom_[vi] = gt.second;
0272   }
0273 }
0274 
0275 void HiPuRhoProducer::calculatePedestal(std::vector<fastjet::PseudoJet> const& coll) {
0276   LogDebug("PileUpSubtractor") << "The subtractor calculating pedestals...\n";
0277 
0278   std::map<int, double> emean2;
0279   std::map<int, int> ntowers;
0280 
0281   int ietaold = -10000;
0282   // Initial values for emean_, emean2, esigma_, ntowers
0283 
0284   for (int vi = 0; vi < nEtaTow_; ++vi) {
0285     vntow_[vi] = initialValue_;
0286     vmean1_[vi] = initialValue_;
0287     vrms1_[vi] = initialValue_;
0288     vrho1_[vi] = initialValue_;
0289 
0290     if (setInitialValue_) {
0291       vmean0_[vi] = initialValue_;
0292       vrms0_[vi] = initialValue_;
0293       vrho0_[vi] = initialValue_;
0294     }
0295   }
0296 
0297   for (int i = ietamin_; i < ietamax_ + 1; i++) {
0298     emean_[i] = 0.;
0299     emean2[i] = 0.;
0300     esigma_[i] = 0.;
0301     ntowers[i] = 0;
0302 
0303     eTop4_[i] = {{0., 0., 0., 0.}};
0304   }
0305 
0306   for (auto const& input_object : coll) {
0307     const CaloTower* originalTower = inputs_[input_object.user_index()];
0308     double original_Et = originalTower->et();
0309     int ieta0 = originalTower->ieta();
0310 
0311     if (original_Et > eTop4_[ieta0][0]) {
0312       eTop4_[ieta0][3] = eTop4_[ieta0][2];
0313       eTop4_[ieta0][2] = eTop4_[ieta0][1];
0314       eTop4_[ieta0][1] = eTop4_[ieta0][0];
0315       eTop4_[ieta0][0] = original_Et;
0316     } else if (original_Et > eTop4_[ieta0][1]) {
0317       eTop4_[ieta0][3] = eTop4_[ieta0][2];
0318       eTop4_[ieta0][2] = eTop4_[ieta0][1];
0319       eTop4_[ieta0][1] = original_Et;
0320     } else if (original_Et > eTop4_[ieta0][2]) {
0321       eTop4_[ieta0][3] = eTop4_[ieta0][2];
0322       eTop4_[ieta0][2] = original_Et;
0323     } else if (original_Et > eTop4_[ieta0][3]) {
0324       eTop4_[ieta0][3] = original_Et;
0325     }
0326 
0327     emean_[ieta0] = emean_[ieta0] + original_Et;
0328     emean2[ieta0] = emean2[ieta0] + original_Et * original_Et;
0329     if (ieta0 - ietaold != 0) {
0330       ntowers[ieta0] = 1;
0331       ietaold = ieta0;
0332     } else {
0333       ntowers[ieta0]++;
0334     }
0335   }
0336 
0337   for (auto const& gt : geomtowers_) {
0338     int it = gt.first;
0339 
0340     int vi = it - 1;
0341 
0342     if (it < 0)
0343       vi = nEtaTow_ + it;
0344 
0345     double e1 = emean_[it];
0346     double e2 = emean2[it];
0347     int nt = gt.second - ntowersWithJets_[it];
0348 
0349     if (vi < nEtaTow_) {
0350       vntow_[vi] = nt;
0351     }
0352 
0353     LogDebug("PileUpSubtractor") << " ieta: " << it << " number of towers: " << nt << " e1: " << e1 << " e2: " << e2
0354                                  << "\n";
0355 
0356     if (nt > 0) {
0357       if (postOrphan_) {
0358         if (nt > (int)minimumTowersFraction_ * (gt.second)) {
0359           emean_[it] = e1 / (double)nt;
0360           double eee = e2 / (double)nt - e1 * e1 / (double)(nt * nt);
0361           if (eee < 0.)
0362             eee = 0.;
0363           esigma_[it] = nSigmaPU_ * sqrt(eee);
0364 
0365           uint32_t numToCheck = std::min(int(eTop4_[it].size()), nt - (int)minimumTowersFraction_ * (gt.second));
0366 
0367           for (unsigned int lI = 0; lI < numToCheck; ++lI) {
0368             if (eTop4_[it][lI] >= emean_[it] + towSigmaCut_ * esigma_[it] && towSigmaCut_ > 0) {
0369               e1 -= eTop4_[it][lI];
0370               nt -= 1;
0371             } else
0372               break;
0373           }
0374 
0375           if (e1 < .000000001)
0376             e1 = 0;
0377         }
0378       }
0379 
0380       emean_[it] = e1 / (double)nt;
0381       double eee = e2 / nt - e1 * e1 / (nt * nt);
0382       if (eee < 0.)
0383         eee = 0.;
0384       esigma_[it] = nSigmaPU_ * sqrt(eee);
0385 
0386       double etaWidth = std::abs(hi::etaedge[abs(it)] - hi::etaedge[abs(it) - 1]);
0387 
0388       int sign = (it < 0) ? -1 : 1;
0389       auto absIt = std::abs(it);
0390       if (it < 0) {
0391         etaEdgeLow_.push_back(sign * hi::etaedge[absIt]);
0392         etaEdgeHi_.push_back(sign * hi::etaedge[absIt - 1]);
0393       } else {
0394         etaEdgeHi_.push_back(sign * hi::etaedge[absIt]);
0395         etaEdgeLow_.push_back(sign * hi::etaedge[absIt - 1]);
0396       }
0397 
0398       if (vi < nEtaTow_) {
0399         vmean1_[vi] = emean_[it];
0400         vrho1_[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom_[vi]));
0401         rho_.push_back(vrho1_[vi]);
0402         rhoM_.push_back(0);
0403         vrms1_[vi] = esigma_[it];
0404         if (vngeom_[vi] == vntow_[vi]) {
0405           vmean0_[vi] = emean_[it];
0406           vrho0_[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom_[vi]));
0407           vrms0_[vi] = esigma_[it];
0408         }
0409         rhoExtra_.push_back(vrho0_[vi]);
0410         nTow_.push_back(vntow_[vi]);
0411       }
0412     } else {
0413       emean_[it] = 0.;
0414       esigma_[it] = 0.;
0415     }
0416     LogDebug("PileUpSubtractor") << " ieta: " << it << " Pedestals: " << emean_[it] << "  " << esigma_[it] << "\n";
0417   }
0418 
0419   postOrphan_ = false;
0420 }
0421 
0422 void HiPuRhoProducer::subtractPedestal(std::vector<fastjet::PseudoJet>& coll) {
0423   LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
0424 
0425   std::vector<fastjet::PseudoJet> newcoll;
0426   for (auto& input_object : coll) {
0427     int index = input_object.user_index();
0428     CaloTower const* itow = inputs_[index];
0429     int it = itow->ieta();
0430 
0431     double original_Et = itow->et();
0432     double etnew = original_Et - emean_.at(it) - esigma_.at(it);
0433     float mScale = etnew / input_object.Et();
0434     if (etnew < 0.)
0435       mScale = 0.;
0436 
0437     math::XYZTLorentzVectorD towP4(
0438         input_object.px() * mScale, input_object.py() * mScale, input_object.pz() * mScale, input_object.e() * mScale);
0439 
0440     input_object.reset_momentum(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
0441     input_object.set_user_index(index);
0442 
0443     if (etnew > 0. && dropZeroTowers_)
0444       newcoll.push_back(input_object);
0445   }
0446 
0447   if (dropZeroTowers_)
0448     coll = newcoll;
0449 }
0450 
0451 void HiPuRhoProducer::calculateOrphanInput(std::vector<fastjet::PseudoJet>& orphanInput) {
0452   LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
0453 
0454   fjInputs_ = fjOriginalInputs_;
0455 
0456   std::vector<int> jettowers;                         // vector of towers indexed by "user_index"
0457   std::map<std::pair<int, int>, int> excludedTowers;  // map of excluded ieta, iphi values
0458   for (auto const& pseudojetTMP : fjJets_) {
0459     EtaPhi jet_etaphi(pseudojetTMP.eta(), pseudojetTMP.phi());
0460 
0461     if (pseudojetTMP.perp() < puPtMin_)
0462       continue;
0463 
0464     for (auto const& im : towermap_) {
0465       double dr2 = reco::deltaR2(im.eta, im.phi, jet_etaphi.first, jet_etaphi.second);
0466       if (dr2 < radiusPU_ * radiusPU_ && !excludedTowers[std::pair(im.ieta, im.iphi)] &&
0467           (geomtowers_[im.ieta] - ntowersWithJets_[im.ieta]) > minimumTowersFraction_ * geomtowers_[im.ieta]) {
0468         ntowersWithJets_[im.ieta]++;
0469         excludedTowers[std::pair(im.ieta, im.iphi)] = 1;
0470       }
0471     }
0472 
0473     for (auto const& input : fjInputs_) {
0474       int index = input.user_index();
0475       const CaloTower* originalTower = inputs_[index];
0476       int ie = originalTower->ieta();
0477       int ip = originalTower->iphi();
0478 
0479       if (excludedTowers[std::pair<int, int>(ie, ip)]) {
0480         jettowers.push_back(index);
0481       }
0482     }
0483   }
0484   // Create a new collections from the towers not included in jets
0485   for (auto const& input : fjInputs_) {
0486     int index = input.user_index();
0487     const CaloTower* originalTower = inputs_[index];
0488     auto itjet = std::find(jettowers.begin(), jettowers.end(), index);
0489     if (itjet == jettowers.end()) {
0490       orphanInput.emplace_back(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy());
0491       orphanInput.back().set_user_index(index);
0492     } else {
0493       towExcludePt_.push_back(originalTower->pt());
0494       towExcludePhi_.push_back(originalTower->phi());
0495       towExcludeEta_.push_back(originalTower->eta());
0496     }
0497   }
0498 
0499   postOrphan_ = true;
0500 }
0501 
0502 void HiPuRhoProducer::putRho(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0503   std::size_t size = etaEdgeLow_.size();
0504 
0505   std::vector<std::pair<std::size_t, double>> order;
0506   for (std::size_t i = 0; i < size; ++i) {
0507     order.emplace_back(i, etaEdgeLow_[i]);
0508   }
0509 
0510   std::sort(
0511       order.begin(), order.end(), [](auto const& pair0, auto const& pair1) { return pair0.second < pair1.second; });
0512 
0513   std::vector<double> sortedEtaEdgeLow(size);
0514   std::vector<double> sortedEtaEdgeHigh(size);
0515 
0516   auto mapToRhoOut = std::make_unique<std::vector<double>>(size);
0517   auto mapToRhoExtraOut = std::make_unique<std::vector<double>>(size);
0518   auto mapToRhoMOut = std::make_unique<std::vector<double>>(size);
0519   auto mapToNTowOut = std::make_unique<std::vector<int>>(size);
0520 
0521   for (std::size_t i = 0; i < size; ++i) {
0522     auto const& index = order[i].first;
0523 
0524     sortedEtaEdgeLow[i] = etaEdgeLow_[index];
0525     sortedEtaEdgeHigh[i] = etaEdgeHi_[index];
0526 
0527     (*mapToRhoOut)[i] = rho_[index];
0528     (*mapToRhoExtraOut)[i] = rhoExtra_[index];
0529     (*mapToRhoMOut)[i] = rhoM_[index];
0530     (*mapToNTowOut)[i] = nTow_[index];
0531   }
0532 
0533   auto mapToRhoMedianOut = std::make_unique<std::vector<double>>(size);
0534 
0535   for (uint32_t index = medianWindowWidth_; index < size - medianWindowWidth_; ++index) {
0536     auto centre = mapToRhoOut->begin() + index;
0537     std::vector<float> rhoRange(centre - medianWindowWidth_, centre + medianWindowWidth_);
0538     std::nth_element(rhoRange.begin(), rhoRange.begin() + medianWindowWidth_, rhoRange.end());
0539     (*mapToRhoMedianOut)[index] = rhoRange[medianWindowWidth_];
0540   }
0541 
0542   auto mapEtaRangesOut = std::make_unique<std::vector<double>>();
0543 
0544   mapEtaRangesOut->assign(sortedEtaEdgeLow.begin(), sortedEtaEdgeLow.end());
0545   mapEtaRangesOut->push_back(sortedEtaEdgeHigh.back());
0546 
0547   auto mapToTowExcludePtOut = std::make_unique<std::vector<double>>(std::move(towExcludePt_));
0548   auto mapToTowExcludePhiOut = std::make_unique<std::vector<double>>(std::move(towExcludePhi_));
0549   auto mapToTowExcludeEtaOut = std::make_unique<std::vector<double>>(std::move(towExcludeEta_));
0550 
0551   iEvent.put(std::move(mapEtaRangesOut), "mapEtaEdges");
0552   iEvent.put(std::move(mapToRhoOut), "mapToRho");
0553   iEvent.put(std::move(mapToRhoMedianOut), "mapToRhoMedian");
0554   iEvent.put(std::move(mapToRhoExtraOut), "mapToRhoExtra");
0555   iEvent.put(std::move(mapToRhoMOut), "mapToRhoM");
0556   iEvent.put(std::move(mapToNTowOut), "mapToNTow");
0557   iEvent.put(std::move(mapToTowExcludePtOut), "mapToTowExcludePt");
0558   iEvent.put(std::move(mapToTowExcludePhiOut), "mapToTowExcludePhi");
0559   iEvent.put(std::move(mapToTowExcludeEtaOut), "mapToTowExcludeEta");
0560 }
0561 
0562 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0563 void HiPuRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0564   edm::ParameterSetDescription desc;
0565   desc.add<edm::InputTag>("src", edm::InputTag("PFTowers"));
0566   desc.add<int>("medianWindowWidth", 2);
0567   desc.add<double>("towSigmaCut", 5.0);
0568   desc.add<double>("puPtMin", 15.0);
0569   desc.add<double>("rParam", 0.3);
0570   desc.add<double>("nSigmaPU", 1.0);
0571   desc.add<double>("radiusPU", 0.5);
0572   desc.add<double>("minimumTowersFraction", 0.7);
0573   desc.add<bool>("dropZeroTowers", true);
0574   descriptions.add("hiPuRhoProducer", desc);
0575 }
0576 
0577 //define this as a plug-in
0578 DEFINE_FWK_MODULE(HiPuRhoProducer);