Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:19 

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 
0136   int ietamax_;                                 // maximum eta in geometry
0137   int ietamin_;                                 // minimum eta in geometry
0138   std::map<int, int> ntowersWithJets_;          // number of towers with jets
0139   std::map<int, int> geomtowers_;               // map of geometry towers to det id
0140   std::map<int, double> esigma_;                // energy sigma
0141   std::map<int, double> emean_;                 // energy mean
0142   std::map<int, std::array<double, 4>> eTop4_;  // energy mean
0143 
0144   typedef std::pair<double, double> EtaPhi;
0145   std::vector<EtaPhiTower> towermap_;
0146 };
0147 
0148 HiPuRhoProducer::HiPuRhoProducer(const edm::ParameterSet& iConfig)
0149     : dropZeroTowers_(iConfig.getParameter<bool>("dropZeroTowers")),
0150       medianWindowWidth_(iConfig.getParameter<int>("medianWindowWidth")),
0151       minimumTowersFraction_(iConfig.getParameter<double>("minimumTowersFraction")),
0152       nSigmaPU_(iConfig.getParameter<double>("nSigmaPU")),
0153       puPtMin_(iConfig.getParameter<double>("puPtMin")),
0154       radiusPU_(iConfig.getParameter<double>("radiusPU")),
0155       rParam_(iConfig.getParameter<double>("rParam")),
0156       towSigmaCut_(iConfig.getParameter<double>("towSigmaCut")),
0157       caloTowerToken_(consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("src"))),
0158       caloGeometryToken_(esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})) {
0159   //register your products
0160   produces<std::vector<double>>("mapEtaEdges");
0161   produces<std::vector<double>>("mapToRho");
0162   produces<std::vector<double>>("mapToRhoMedian");
0163   produces<std::vector<double>>("mapToRhoExtra");
0164   produces<std::vector<double>>("mapToRhoM");
0165   produces<std::vector<int>>("mapToNTow");
0166   produces<std::vector<double>>("mapToTowExcludePt");
0167   produces<std::vector<double>>("mapToTowExcludePhi");
0168   produces<std::vector<double>>("mapToTowExcludeEta");
0169 }
0170 
0171 // ------------ method called to produce the data  ------------
0172 void HiPuRhoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0173   setupGeometryMap(iEvent, iSetup);
0174 
0175   for (int i = ietamin_; i < ietamax_ + 1; i++) {
0176     ntowersWithJets_[i] = 0;
0177   }
0178 
0179   auto const& inputView = iEvent.get(caloTowerToken_);
0180   inputs_.reserve(inputView.size());
0181   for (auto const& input : inputView)
0182     inputs_.push_back(&input);
0183 
0184   fjInputs_.reserve(inputs_.size());
0185   inputTowers();
0186   fjOriginalInputs_ = fjInputs_;
0187   setInitialValue_ = true;
0188   calculatePedestal(fjInputs_);
0189   subtractPedestal(fjInputs_);
0190 
0191   fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, rParam_);
0192   fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, *fjJetDefinition_);
0193   fjJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(puPtMin_));
0194 
0195   etaEdgeLow_.clear();
0196   etaEdgeHi_.clear();
0197   etaEdges_.clear();
0198 
0199   rho_.clear();
0200   rhoExtra_.clear();
0201   rhoM_.clear();
0202   nTow_.clear();
0203 
0204   towExcludePt_.clear();
0205   towExcludePhi_.clear();
0206   towExcludeEta_.clear();
0207 
0208   setInitialValue_ = false;
0209   std::vector<fastjet::PseudoJet> orphanInput;
0210   calculateOrphanInput(orphanInput);
0211   calculatePedestal(orphanInput);
0212   putRho(iEvent, iSetup);
0213 
0214   inputs_.clear();
0215   fjInputs_.clear();
0216   fjJets_.clear();
0217 }
0218 
0219 void HiPuRhoProducer::inputTowers() {
0220   int index = -1;
0221   for (auto const& input : inputs_) {
0222     index++;
0223 
0224     if (edm::isNotFinite(input->pt()))
0225       continue;
0226     if (input->pt() < 100 * std::numeric_limits<double>::epsilon())
0227       continue;
0228 
0229     fjInputs_.push_back(fastjet::PseudoJet(input->px(), input->py(), input->pz(), input->energy()));
0230     fjInputs_.back().set_user_index(index);
0231   }
0232 }
0233 
0234 void HiPuRhoProducer::setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0235   LogDebug("PileUpSubtractor") << "The subtractor setting up geometry...\n";
0236   const auto& pG = iSetup.getData(caloGeometryToken_);
0237   geo_ = &pG;
0238   std::vector<DetId> alldid = geo_->getValidDetIds();
0239   int ietaold = -10000;
0240   ietamax_ = -10000;
0241   ietamin_ = 10000;
0242   towermap_.clear();
0243 
0244   for (auto const& did : alldid) {
0245     if (did.det() == DetId::Hcal) {
0246       HcalDetId hid = HcalDetId(did);
0247       towermap_.push_back({hid.ieta(), hid.iphi(), geo_->getPosition(did).eta(), geo_->getPosition(did).phi()});
0248       if (hid.ieta() != ietaold) {
0249         ietaold = hid.ieta();
0250         geomtowers_[hid.ieta()] = 1;
0251         if (hid.ieta() > ietamax_)
0252           ietamax_ = hid.ieta();
0253         if (hid.ieta() < ietamin_)
0254           ietamin_ = hid.ieta();
0255       } else {
0256         geomtowers_[hid.ieta()]++;
0257       }
0258     }
0259   }
0260 
0261   for (auto const& gt : geomtowers_) {
0262     int it = gt.first;
0263     int vi = it - 1;
0264 
0265     if (it < 0)
0266       vi = nEtaTow_ + it;
0267 
0268     if (vi < nEtaTow_)
0269       vngeom_[vi] = gt.second;
0270   }
0271 }
0272 
0273 void HiPuRhoProducer::calculatePedestal(std::vector<fastjet::PseudoJet> const& coll) {
0274   LogDebug("PileUpSubtractor") << "The subtractor calculating pedestals...\n";
0275 
0276   std::map<int, double> emean2;
0277   std::map<int, int> ntowers;
0278 
0279   int ietaold = -10000;
0280   // Initial values for emean_, emean2, esigma_, ntowers
0281 
0282   for (int vi = 0; vi < nEtaTow_; ++vi) {
0283     vntow_[vi] = initialValue_;
0284     vmean1_[vi] = initialValue_;
0285     vrms1_[vi] = initialValue_;
0286     vrho1_[vi] = initialValue_;
0287 
0288     if (setInitialValue_) {
0289       vmean0_[vi] = initialValue_;
0290       vrms0_[vi] = initialValue_;
0291       vrho0_[vi] = initialValue_;
0292     }
0293   }
0294 
0295   for (int i = ietamin_; i < ietamax_ + 1; i++) {
0296     emean_[i] = 0.;
0297     emean2[i] = 0.;
0298     esigma_[i] = 0.;
0299     ntowers[i] = 0;
0300 
0301     eTop4_[i] = {{0., 0., 0., 0.}};
0302   }
0303 
0304   for (auto const& input_object : coll) {
0305     const CaloTower* originalTower = inputs_[input_object.user_index()];
0306     double original_Et = originalTower->et();
0307     int ieta0 = originalTower->ieta();
0308 
0309     if (original_Et > eTop4_[ieta0][0]) {
0310       eTop4_[ieta0][3] = eTop4_[ieta0][2];
0311       eTop4_[ieta0][2] = eTop4_[ieta0][1];
0312       eTop4_[ieta0][1] = eTop4_[ieta0][0];
0313       eTop4_[ieta0][0] = original_Et;
0314     } else if (original_Et > eTop4_[ieta0][1]) {
0315       eTop4_[ieta0][3] = eTop4_[ieta0][2];
0316       eTop4_[ieta0][2] = eTop4_[ieta0][1];
0317       eTop4_[ieta0][1] = original_Et;
0318     } else if (original_Et > eTop4_[ieta0][2]) {
0319       eTop4_[ieta0][3] = eTop4_[ieta0][2];
0320       eTop4_[ieta0][2] = original_Et;
0321     } else if (original_Et > eTop4_[ieta0][3]) {
0322       eTop4_[ieta0][3] = original_Et;
0323     }
0324 
0325     emean_[ieta0] = emean_[ieta0] + original_Et;
0326     emean2[ieta0] = emean2[ieta0] + original_Et * original_Et;
0327     if (ieta0 - ietaold != 0) {
0328       ntowers[ieta0] = 1;
0329       ietaold = ieta0;
0330     } else {
0331       ntowers[ieta0]++;
0332     }
0333   }
0334 
0335   for (auto const& gt : geomtowers_) {
0336     int it = gt.first;
0337 
0338     int vi = it - 1;
0339 
0340     if (it < 0)
0341       vi = nEtaTow_ + it;
0342 
0343     double e1 = emean_[it];
0344     double e2 = emean2[it];
0345     int nt = gt.second - ntowersWithJets_[it];
0346 
0347     if (vi < nEtaTow_) {
0348       vntow_[vi] = nt;
0349     }
0350 
0351     LogDebug("PileUpSubtractor") << " ieta: " << it << " number of towers: " << nt << " e1: " << e1 << " e2: " << e2
0352                                  << "\n";
0353 
0354     if (nt > 0) {
0355       if (postOrphan_) {
0356         if (nt > (int)minimumTowersFraction_ * (gt.second)) {
0357           emean_[it] = e1 / (double)nt;
0358           double eee = e2 / (double)nt - e1 * e1 / (double)(nt * nt);
0359           if (eee < 0.)
0360             eee = 0.;
0361           esigma_[it] = nSigmaPU_ * sqrt(eee);
0362 
0363           uint32_t numToCheck = std::min(int(eTop4_[it].size()), nt - (int)minimumTowersFraction_ * (gt.second));
0364 
0365           for (unsigned int lI = 0; lI < numToCheck; ++lI) {
0366             if (eTop4_[it][lI] >= emean_[it] + towSigmaCut_ * esigma_[it] && towSigmaCut_ > 0) {
0367               e1 -= eTop4_[it][lI];
0368               nt -= 1;
0369             } else
0370               break;
0371           }
0372 
0373           if (e1 < .000000001)
0374             e1 = 0;
0375         }
0376       }
0377 
0378       emean_[it] = e1 / (double)nt;
0379       double eee = e2 / nt - e1 * e1 / (nt * nt);
0380       if (eee < 0.)
0381         eee = 0.;
0382       esigma_[it] = nSigmaPU_ * sqrt(eee);
0383 
0384       double etaWidth = std::abs(hi::etaedge[abs(it)] - hi::etaedge[abs(it) - 1]);
0385 
0386       int sign = (it < 0) ? -1 : 1;
0387       auto absIt = std::abs(it);
0388       if (it < 0) {
0389         etaEdgeLow_.push_back(sign * hi::etaedge[absIt]);
0390         etaEdgeHi_.push_back(sign * hi::etaedge[absIt - 1]);
0391       } else {
0392         etaEdgeHi_.push_back(sign * hi::etaedge[absIt]);
0393         etaEdgeLow_.push_back(sign * hi::etaedge[absIt - 1]);
0394       }
0395 
0396       if (vi < nEtaTow_) {
0397         vmean1_[vi] = emean_[it];
0398         vrho1_[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom_[vi]));
0399         rho_.push_back(vrho1_[vi]);
0400         rhoM_.push_back(0);
0401         vrms1_[vi] = esigma_[it];
0402         if (vngeom_[vi] == vntow_[vi]) {
0403           vmean0_[vi] = emean_[it];
0404           vrho0_[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom_[vi]));
0405           vrms0_[vi] = esigma_[it];
0406         }
0407         rhoExtra_.push_back(vrho0_[vi]);
0408         nTow_.push_back(vntow_[vi]);
0409       }
0410     } else {
0411       emean_[it] = 0.;
0412       esigma_[it] = 0.;
0413     }
0414     LogDebug("PileUpSubtractor") << " ieta: " << it << " Pedestals: " << emean_[it] << "  " << esigma_[it] << "\n";
0415   }
0416 
0417   postOrphan_ = false;
0418 }
0419 
0420 void HiPuRhoProducer::subtractPedestal(std::vector<fastjet::PseudoJet>& coll) {
0421   LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
0422 
0423   std::vector<fastjet::PseudoJet> newcoll;
0424   for (auto& input_object : coll) {
0425     int index = input_object.user_index();
0426     CaloTower const* itow = inputs_[index];
0427     int it = itow->ieta();
0428 
0429     double original_Et = itow->et();
0430     double etnew = original_Et - emean_.at(it) - esigma_.at(it);
0431     float mScale = etnew / input_object.Et();
0432     if (etnew < 0.)
0433       mScale = 0.;
0434 
0435     math::XYZTLorentzVectorD towP4(
0436         input_object.px() * mScale, input_object.py() * mScale, input_object.pz() * mScale, input_object.e() * mScale);
0437 
0438     input_object.reset_momentum(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
0439     input_object.set_user_index(index);
0440 
0441     if (etnew > 0. && dropZeroTowers_)
0442       newcoll.push_back(input_object);
0443   }
0444 
0445   if (dropZeroTowers_)
0446     coll = newcoll;
0447 }
0448 
0449 void HiPuRhoProducer::calculateOrphanInput(std::vector<fastjet::PseudoJet>& orphanInput) {
0450   LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
0451 
0452   fjInputs_ = fjOriginalInputs_;
0453 
0454   std::vector<int> jettowers;                         // vector of towers indexed by "user_index"
0455   std::map<std::pair<int, int>, int> excludedTowers;  // map of excluded ieta, iphi values
0456   for (auto const& pseudojetTMP : fjJets_) {
0457     EtaPhi jet_etaphi(pseudojetTMP.eta(), pseudojetTMP.phi());
0458 
0459     if (pseudojetTMP.perp() < puPtMin_)
0460       continue;
0461 
0462     for (auto const& im : towermap_) {
0463       double dr2 = reco::deltaR2(im.eta, im.phi, jet_etaphi.first, jet_etaphi.second);
0464       if (dr2 < radiusPU_ * radiusPU_ && !excludedTowers[std::pair(im.ieta, im.iphi)] &&
0465           (geomtowers_[im.ieta] - ntowersWithJets_[im.ieta]) > minimumTowersFraction_ * geomtowers_[im.ieta]) {
0466         ntowersWithJets_[im.ieta]++;
0467         excludedTowers[std::pair(im.ieta, im.iphi)] = 1;
0468       }
0469     }
0470 
0471     for (auto const& input : fjInputs_) {
0472       int index = input.user_index();
0473       const CaloTower* originalTower = inputs_[index];
0474       int ie = originalTower->ieta();
0475       int ip = originalTower->iphi();
0476 
0477       if (excludedTowers[std::pair<int, int>(ie, ip)]) {
0478         jettowers.push_back(index);
0479       }
0480     }
0481   }
0482   // Create a new collections from the towers not included in jets
0483   for (auto const& input : fjInputs_) {
0484     int index = input.user_index();
0485     const CaloTower* originalTower = inputs_[index];
0486     auto itjet = std::find(jettowers.begin(), jettowers.end(), index);
0487     if (itjet == jettowers.end()) {
0488       orphanInput.emplace_back(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy());
0489       orphanInput.back().set_user_index(index);
0490     } else {
0491       towExcludePt_.push_back(originalTower->pt());
0492       towExcludePhi_.push_back(originalTower->phi());
0493       towExcludeEta_.push_back(originalTower->eta());
0494     }
0495   }
0496 
0497   postOrphan_ = true;
0498 }
0499 
0500 void HiPuRhoProducer::putRho(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0501   std::size_t size = etaEdgeLow_.size();
0502 
0503   std::vector<std::pair<std::size_t, double>> order;
0504   for (std::size_t i = 0; i < size; ++i) {
0505     order.emplace_back(i, etaEdgeLow_[i]);
0506   }
0507 
0508   std::sort(
0509       order.begin(), order.end(), [](auto const& pair0, auto const& pair1) { return pair0.second < pair1.second; });
0510 
0511   std::vector<double> sortedEtaEdgeLow(size);
0512   std::vector<double> sortedEtaEdgeHigh(size);
0513 
0514   auto mapToRhoOut = std::make_unique<std::vector<double>>(size);
0515   auto mapToRhoExtraOut = std::make_unique<std::vector<double>>(size);
0516   auto mapToRhoMOut = std::make_unique<std::vector<double>>(size);
0517   auto mapToNTowOut = std::make_unique<std::vector<int>>(size);
0518 
0519   for (std::size_t i = 0; i < size; ++i) {
0520     auto const& index = order[i].first;
0521 
0522     sortedEtaEdgeLow[i] = etaEdgeLow_[index];
0523     sortedEtaEdgeHigh[i] = etaEdgeHi_[index];
0524 
0525     (*mapToRhoOut)[i] = rho_[index];
0526     (*mapToRhoExtraOut)[i] = rhoExtra_[index];
0527     (*mapToRhoMOut)[i] = rhoM_[index];
0528     (*mapToNTowOut)[i] = nTow_[index];
0529   }
0530 
0531   auto mapToRhoMedianOut = std::make_unique<std::vector<double>>(size);
0532 
0533   for (uint32_t index = medianWindowWidth_; index < size - medianWindowWidth_; ++index) {
0534     auto centre = mapToRhoOut->begin() + index;
0535     std::vector<float> rhoRange(centre - medianWindowWidth_, centre + medianWindowWidth_);
0536     std::nth_element(rhoRange.begin(), rhoRange.begin() + medianWindowWidth_, rhoRange.end());
0537     (*mapToRhoMedianOut)[index] = rhoRange[medianWindowWidth_];
0538   }
0539 
0540   auto mapEtaRangesOut = std::make_unique<std::vector<double>>();
0541 
0542   mapEtaRangesOut->assign(sortedEtaEdgeLow.begin(), sortedEtaEdgeLow.end());
0543   mapEtaRangesOut->push_back(sortedEtaEdgeHigh.back());
0544 
0545   auto mapToTowExcludePtOut = std::make_unique<std::vector<double>>(std::move(towExcludePt_));
0546   auto mapToTowExcludePhiOut = std::make_unique<std::vector<double>>(std::move(towExcludePhi_));
0547   auto mapToTowExcludeEtaOut = std::make_unique<std::vector<double>>(std::move(towExcludeEta_));
0548 
0549   iEvent.put(std::move(mapEtaRangesOut), "mapEtaEdges");
0550   iEvent.put(std::move(mapToRhoOut), "mapToRho");
0551   iEvent.put(std::move(mapToRhoMedianOut), "mapToRhoMedian");
0552   iEvent.put(std::move(mapToRhoExtraOut), "mapToRhoExtra");
0553   iEvent.put(std::move(mapToRhoMOut), "mapToRhoM");
0554   iEvent.put(std::move(mapToNTowOut), "mapToNTow");
0555   iEvent.put(std::move(mapToTowExcludePtOut), "mapToTowExcludePt");
0556   iEvent.put(std::move(mapToTowExcludePhiOut), "mapToTowExcludePhi");
0557   iEvent.put(std::move(mapToTowExcludeEtaOut), "mapToTowExcludeEta");
0558 }
0559 
0560 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0561 void HiPuRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0562   edm::ParameterSetDescription desc;
0563   desc.add<edm::InputTag>("src", edm::InputTag("PFTowers"));
0564   desc.add<int>("medianWindowWidth", 2);
0565   desc.add<double>("towSigmaCut", 5.0);
0566   desc.add<double>("puPtMin", 15.0);
0567   desc.add<double>("rParam", 0.3);
0568   desc.add<double>("nSigmaPU", 1.0);
0569   desc.add<double>("radiusPU", 0.5);
0570   desc.add<double>("minimumTowersFraction", 0.7);
0571   desc.add<bool>("dropZeroTowers", true);
0572   descriptions.add("hiPuRhoProducer", desc);
0573 }
0574 
0575 //define this as a plug-in
0576 DEFINE_FWK_MODULE(HiPuRhoProducer);
[ 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!