JetProducers/plugins/VirtualJetProducer.cc

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // VirtualJetProducer
0004 // ------------------
0005 //
0006 //            04/21/2009 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
0007 ////////////////////////////////////////////////////////////////////////////////
0008
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0011 #include "RecoJets/JetProducers/plugins/VirtualJetProducer.h"
0012 #include "RecoJets/JetProducers/interface/BackgroundEstimator.h"
0013 #include "RecoJets/JetProducers/interface/VirtualJetProducerHelper.h"
0014
0015 #include "DataFormats/Common/interface/RefProd.h"
0016 #include "DataFormats/Common/interface/Ref.h"
0017 #include "DataFormats/Common/interface/RefVector.h"
0018
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/ESHandle.h"
0022 #include "FWCore/Utilities/interface/Exception.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/Utilities/interface/isFinite.h"
0026 #include "FWCore/Framework/interface/ConsumesCollector.h"
0027
0028 #include "DataFormats/Common/interface/View.h"
0029 #include "DataFormats/Common/interface/Handle.h"
0030 #include "DataFormats/VertexReco/interface/Vertex.h"
0031 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0032 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0033 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0034 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0035 #include "DataFormats/JetReco/interface/BasicJetCollection.h"
0036 #include "DataFormats/JetReco/interface/TrackJetCollection.h"
0037 #include "DataFormats/JetReco/interface/PFClusterJetCollection.h"
0038 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0039 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0040 #include "DataFormats/Math/interface/deltaR.h"
0041 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0042
0043 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0044 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0045
0046 #include "fastjet/SISConePlugin.hh"
0047 #include "fastjet/CMSIterativeConePlugin.hh"
0048 #include "fastjet/ATLASConePlugin.hh"
0049 #include "fastjet/CDFMidPointPlugin.hh"
0050
0051 #include <iostream>
0052 #include <memory>
0053 #include <algorithm>
0054 #include <limits>
0055 #include <cmath>
0056 #include <vdt/vdtMath.h>
0057
0058 using namespace std;
0059 using namespace edm;
0060
0061 namespace reco {
0062   namespace helper {
0063     struct GreaterByPtPseudoJet {
0064       bool operator()(const fastjet::PseudoJet& t1, const fastjet::PseudoJet& t2) const {
0065         return t1.perp2() > t2.perp2();
0066       }
0067     };
0068
0069   }  // namespace helper
0070 }  // namespace reco
0071
0072 //______________________________________________________________________________
0073 const char* const VirtualJetProducer::JetType::names[] = {
0074     "BasicJet", "GenJet", "CaloJet", "PFJet", "TrackJet", "PFClusterJet"};
0075
0076 //______________________________________________________________________________
0077 VirtualJetProducer::JetType::Type VirtualJetProducer::JetType::byName(const string& name) {
0078   const char* const* pos = std::find(names, names + LastJetType, name);
0079   if (pos == names + LastJetType) {
0080     std::string errorMessage = "Requested jetType not supported: " + name + "\n";
0081     throw cms::Exception("Configuration", errorMessage);
0082   }
0083   return (Type)(pos - names);
0084 }
0085
0086 void VirtualJetProducer::makeProduces(std::string alias, std::string tag) {
0087   if (writeCompound_) {
0088     produces<reco::BasicJetCollection>();
0089   }
0090
0091   if (writeJetsWithConst_) {
0092     produces<reco::PFCandidateCollection>(tag).setBranchAlias(alias);
0093     produces<reco::PFJetCollection>();
0094   } else {
0095     if (makeCaloJet(jetTypeE)) {
0096       produces<reco::CaloJetCollection>(tag).setBranchAlias(alias);
0097     } else if (makePFJet(jetTypeE)) {
0098       produces<reco::PFJetCollection>(tag).setBranchAlias(alias);
0099     } else if (makeGenJet(jetTypeE)) {
0100       produces<reco::GenJetCollection>(tag).setBranchAlias(alias);
0101     } else if (makeTrackJet(jetTypeE)) {
0102       produces<reco::TrackJetCollection>(tag).setBranchAlias(alias);
0103     } else if (makePFClusterJet(jetTypeE)) {
0104       produces<reco::PFClusterJetCollection>(tag).setBranchAlias(alias);
0105     } else if (makeBasicJet(jetTypeE)) {
0106       produces<reco::BasicJetCollection>(tag).setBranchAlias(alias);
0107     }
0108   }
0109 }
0110
0111 ////////////////////////////////////////////////////////////////////////////////
0112 // construction / destruction
0113 ////////////////////////////////////////////////////////////////////////////////
0114
0115 //______________________________________________________________________________
0116 VirtualJetProducer::VirtualJetProducer(const edm::ParameterSet& iConfig) {
0117   moduleLabel_ = iConfig.getParameter<string>("@module_label");
0118   src_ = iConfig.getParameter<edm::InputTag>("src");
0119   srcPVs_ = iConfig.getParameter<edm::InputTag>("srcPVs");
0120   jetType_ = iConfig.getParameter<string>("jetType");
0121   jetAlgorithm_ = iConfig.getParameter<string>("jetAlgorithm");
0122   rParam_ = iConfig.getParameter<double>("rParam");
0123   inputEtMin_ = iConfig.getParameter<double>("inputEtMin");
0124   inputEMin_ = iConfig.getParameter<double>("inputEMin");
0125   jetPtMin_ = iConfig.getParameter<double>("jetPtMin");
0126   doPVCorrection_ = iConfig.getParameter<bool>("doPVCorrection");
0127   doAreaFastjet_ = iConfig.getParameter<bool>("doAreaFastjet");
0128   doRhoFastjet_ = iConfig.getParameter<bool>("doRhoFastjet");
0129   jetCollInstanceName_ = iConfig.getParameter<string>("jetCollInstanceName");
0130   doPUOffsetCorr_ = iConfig.getParameter<bool>("doPUOffsetCorr");
0131   puSubtractorName_ = iConfig.getParameter<string>("subtractorName");
0132   useExplicitGhosts_ =
0133       iConfig.getParameter<bool>("useExplicitGhosts");  // use explicit ghosts in the fastjet clustering sequence?
0134   doAreaDiskApprox_ = iConfig.getParameter<bool>("doAreaDiskApprox");
0135   voronoiRfact_ = iConfig.getParameter<double>(
0136       "voronoiRfact");  // Voronoi-based area calculation allows for an empirical scale factor
0137   rhoEtaMax_ = iConfig.getParameter<double>(
0138       "Rho_EtaMax");  // do fasjet area / rho calcluation? => accept corresponding parameters
0139   ghostEtaMax_ = iConfig.getParameter<double>("Ghost_EtaMax");
0140   activeAreaRepeats_ = iConfig.getParameter<int>("Active_Area_Repeats");
0141   ghostArea_ = iConfig.getParameter<double>("GhostArea");
0142   restrictInputs_ = iConfig.getParameter<bool>("restrictInputs");  // restrict inputs to first "maxInputs" towers?
0143   maxInputs_ = iConfig.getParameter<unsigned int>("maxInputs");
0144   writeCompound_ = iConfig.getParameter<bool>(
0145       "writeCompound");  // Check to see if we are writing compound jets for substructure and jet grooming
0146   writeJetsWithConst_ = iConfig.getParameter<bool>("writeJetsWithConst");  //write subtracted jet constituents
0147   doFastJetNonUniform_ = iConfig.getParameter<bool>("doFastJetNonUniform");
0148   puCenters_ = iConfig.getParameter<vector<double>>("puCenters");
0149   puWidth_ = iConfig.getParameter<double>("puWidth");
0150   nExclude_ = iConfig.getParameter<unsigned int>("nExclude");
0151   useDeterministicSeed_ = iConfig.getParameter<bool>("useDeterministicSeed");
0152   minSeed_ = iConfig.getParameter<unsigned int>("minSeed");
0153   verbosity_ = iConfig.getParameter<int>("verbosity");
0154   applyWeight_ = iConfig.getParameter<bool>("applyWeight");
0155   if (applyWeight_) {
0156     edm::InputTag srcWeights = iConfig.getParameter<edm::InputTag>("srcWeights");
0157     if (srcWeights.label() == src_.label())
0158       LogWarning("VirtualJetProducer")
0159           << "Particle and weights collection have the same label. You may be applying the same weights twice.\n";
0160     if (!srcWeights.label().empty())
0161       input_weights_token_ = consumes<edm::ValueMap<float>>(srcWeights);
0162   }
0163
0164   anomalousTowerDef_ = std::make_unique<AnomalousTower>(iConfig);
0165
0166   input_vertex_token_ = consumes<reco::VertexCollection>(srcPVs_);
0167   input_candidateview_token_ = consumes<reco::CandidateView>(src_);
0168   input_candidatefwdptr_token_ =
0169       consumes<vector<edm::FwdPtr<reco::PFCandidate>>>(iConfig.getParameter<edm::InputTag>("src"));
0170   input_packedcandidatefwdptr_token_ =
0171       consumes<vector<edm::FwdPtr<pat::PackedCandidate>>>(iConfig.getParameter<edm::InputTag>("src"));
0172   input_gencandidatefwdptr_token_ =
0173       consumes<vector<edm::FwdPtr<reco::GenParticle>>>(iConfig.getParameter<edm::InputTag>("src"));
0174   input_packedgencandidatefwdptr_token_ =
0175       consumes<vector<edm::FwdPtr<pat::PackedGenParticle>>>(iConfig.getParameter<edm::InputTag>("src"));
0176
0177   //
0178   // additional parameters to think about:
0179   // - overlap threshold (set to 0.75 for the time being)
0180   // - p parameter for generalized kT (set to -2 for the time being)
0181   // - fastjet PU subtraction parameters (not yet considered)
0182   //
0183   if (jetAlgorithm_ == "Kt")
0184     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::kt_algorithm, rParam_);
0185
0186   else if (jetAlgorithm_ == "CambridgeAachen")
0187     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::cambridge_algorithm, rParam_);
0188
0189   else if (jetAlgorithm_ == "AntiKt")
0190     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, rParam_);
0191
0192   else if (jetAlgorithm_ == "GeneralizedKt")
0193     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(fastjet::genkt_algorithm, rParam_, -2);
0194
0195   else if (jetAlgorithm_ == "SISCone") {
0196     fjPlugin_ = PluginPtr(new fastjet::SISConePlugin(rParam_, 0.75, 0, 0.0, false, fastjet::SISConePlugin::SM_pttilde));
0197     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(&*fjPlugin_);
0198
0199   } else if (jetAlgorithm_ == "IterativeCone") {
0200     fjPlugin_ = PluginPtr(new fastjet::CMSIterativeConePlugin(rParam_, 1.0));
0201     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(&*fjPlugin_);
0202
0203   } else if (jetAlgorithm_ == "CDFMidPoint") {
0204     fjPlugin_ = PluginPtr(new fastjet::CDFMidPointPlugin(rParam_, 0.75));
0205     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(&*fjPlugin_);
0206
0207   } else if (jetAlgorithm_ == "ATLASCone") {
0208     fjPlugin_ = PluginPtr(new fastjet::ATLASConePlugin(rParam_));
0209     fjJetDefinition_ = std::make_shared<fastjet::JetDefinition>(&*fjPlugin_);
0210
0211   } else {
0212     throw cms::Exception("Invalid jetAlgorithm") << "Jet algorithm for VirtualJetProducer is invalid, Abort!\n";
0213   }
0214
0215   jetTypeE = JetType::byName(jetType_);
0216
0217   if (JetType::CaloJet == jetTypeE) {
0218     geometry_token_ = esConsumes();
0219     topology_token_ = esConsumes();
0220   }
0221
0222   if (doPUOffsetCorr_) {
0223     if (puSubtractorName_.empty()) {
0224       LogWarning("VirtualJetProducer")
0225           << "Pile Up correction on; however, pile up type is not specified. Using default... \n";
0226       subtractor_ = std::make_shared<PileUpSubtractor>(iConfig, consumesCollector());
0227     } else
0228       subtractor_ = std::shared_ptr<PileUpSubtractor>(
0229           PileUpSubtractorFactory::get()->create(puSubtractorName_, iConfig, consumesCollector()));
0230   }
0231
0232   // do approximate disk-based area calculation => warn if conflicting request
0233   if (doAreaDiskApprox_ && doAreaFastjet_)
0234     throw cms::Exception("Conflicting area calculations")
0235         << "Both the calculation of jet area via fastjet and via an analytical disk approximation have been requested. "
0236            "Please decide on one.\n";
0237
0238   if (doAreaFastjet_ || doRhoFastjet_) {
0239     if (voronoiRfact_ <= 0) {
0240       fjActiveArea_ = std::make_shared<fastjet::GhostedAreaSpec>(ghostEtaMax_, activeAreaRepeats_, ghostArea_);
0241
0242       if (!useExplicitGhosts_) {
0243         fjAreaDefinition_ = std::make_shared<fastjet::AreaDefinition>(fastjet::active_area, *fjActiveArea_);
0244       } else {
0245         fjAreaDefinition_ =
0246             std::make_shared<fastjet::AreaDefinition>(fastjet::active_area_explicit_ghosts, *fjActiveArea_);
0247       }
0248     }
0249     fjSelector_ = std::make_shared<fastjet::Selector>(fastjet::SelectorAbsRapMax(rhoEtaMax_));
0250   }
0251
0252   if ((doFastJetNonUniform_) && (puCenters_.empty()))
0253     throw cms::Exception("doFastJetNonUniform")
0254         << "Parameter puCenters for doFastJetNonUniform is not defined." << std::endl;
0255
0256   // make the "produces" statements
0257   makeProduces(moduleLabel_, jetCollInstanceName_);
0258   produces<vector<double>>("rhos");
0259   produces<vector<double>>("sigmas");
0260   produces<double>("rho");
0261   produces<double>("sigma");
0262 }
0263
0264 //______________________________________________________________________________
0265 VirtualJetProducer::~VirtualJetProducer() {}
0266
0267 ////////////////////////////////////////////////////////////////////////////////
0268 // implementation of member functions
0269 ////////////////////////////////////////////////////////////////////////////////
0270
0271 //______________________________________________________________________________
0272 void VirtualJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0273   // If requested, set the fastjet random seed to a deterministic function
0274   // of the run/lumi/event.
0275   // NOTE!!! The fastjet random number sequence is a global singleton.
0276   // Thus, we have to create an object and get access to the global singleton
0277   // in order to change it.
0278   if (useDeterministicSeed_) {
0279     fastjet::GhostedAreaSpec gas;
0280     std::vector<int> seeds(2);
0281     unsigned int runNum_uint = static_cast<unsigned int>(iEvent.id().run());
0282     unsigned int evNum_uint = static_cast<unsigned int>(iEvent.id().event());
0283     seeds[0] = std::max(runNum_uint, minSeed_ + 3) + 3 * evNum_uint;
0284     seeds[1] = std::max(runNum_uint, minSeed_ + 5) + 5 * evNum_uint;
0285     gas.set_random_status(seeds);
0286   }
0287
0288   LogDebug("VirtualJetProducer") << "Entered produce\n";
0289   //determine signal vertex2
0290   vertex_ = reco::Jet::Point(0, 0, 0);
0291   if ((makeCaloJet(jetTypeE) || makePFJet(jetTypeE)) && doPVCorrection_) {
0292     LogDebug("VirtualJetProducer") << "Adding PV info\n";
0293     edm::Handle<reco::VertexCollection> pvCollection;
0294     iEvent.getByToken(input_vertex_token_, pvCollection);
0295     if (!pvCollection->empty())
0296       vertex_ = pvCollection->begin()->position();
0297   }
0298
0299   // Get Weights Collection
0300   if ((applyWeight_) && (!input_weights_token_.isUninitialized()))
0301     weights_ = iEvent.get(input_weights_token_);
0302
0303   // For Pileup subtraction using offset correction:
0304   // set up geometry map
0305   if (doPUOffsetCorr_) {
0306     subtractor_->setupGeometryMap(iEvent, iSetup);
0307   }
0308
0309   // clear data
0310   LogDebug("VirtualJetProducer") << "Clear data\n";
0311   fjInputs_.clear();
0312   fjJets_.clear();
0313   inputs_.clear();
0314
0315   // get inputs and convert them to the fastjet format (fastjet::PeudoJet)
0316   edm::Handle<reco::CandidateView> inputsHandle;
0317
0318   edm::Handle<std::vector<edm::FwdPtr<reco::PFCandidate>>> pfinputsHandleAsFwdPtr;
0319   edm::Handle<std::vector<edm::FwdPtr<pat::PackedCandidate>>> packedinputsHandleAsFwdPtr;
0320   edm::Handle<std::vector<edm::FwdPtr<reco::GenParticle>>> geninputsHandleAsFwdPtr;
0321   edm::Handle<std::vector<edm::FwdPtr<pat::PackedGenParticle>>> packedgeninputsHandleAsFwdPtr;
0322
0323   bool isView = iEvent.getByToken(input_candidateview_token_, inputsHandle);
0324   if (isView) {
0325     if (inputsHandle->empty()) {
0326       output(iEvent, iSetup);
0327       return;
0328     }
0329     for (size_t i = 0; i < inputsHandle->size(); ++i) {
0330       inputs_.push_back(inputsHandle->ptrAt(i));
0331     }
0332   } else {
0333     bool isPF = iEvent.getByToken(input_candidatefwdptr_token_, pfinputsHandleAsFwdPtr);
0334     bool isPFFwdPtr = iEvent.getByToken(input_packedcandidatefwdptr_token_, packedinputsHandleAsFwdPtr);
0335     bool isGen = iEvent.getByToken(input_gencandidatefwdptr_token_, geninputsHandleAsFwdPtr);
0336     bool isGenFwdPtr = iEvent.getByToken(input_packedgencandidatefwdptr_token_, packedgeninputsHandleAsFwdPtr);
0337
0338     if (isPF) {
0339       if (pfinputsHandleAsFwdPtr->empty()) {
0340         output(iEvent, iSetup);
0341         return;
0342       }
0343       for (size_t i = 0; i < pfinputsHandleAsFwdPtr->size(); ++i) {
0344         if ((*pfinputsHandleAsFwdPtr)[i].ptr().isAvailable()) {
0345           inputs_.push_back((*pfinputsHandleAsFwdPtr)[i].ptr());
0346         } else if ((*pfinputsHandleAsFwdPtr)[i].backPtr().isAvailable()) {
0347           inputs_.push_back((*pfinputsHandleAsFwdPtr)[i].backPtr());
0348         }
0349       }
0350     } else if (isPFFwdPtr) {
0351       if (packedinputsHandleAsFwdPtr->empty()) {
0352         output(iEvent, iSetup);
0353         return;
0354       }
0355       for (size_t i = 0; i < packedinputsHandleAsFwdPtr->size(); ++i) {
0356         if ((*packedinputsHandleAsFwdPtr)[i].ptr().isAvailable()) {
0357           inputs_.push_back((*packedinputsHandleAsFwdPtr)[i].ptr());
0358         } else if ((*packedinputsHandleAsFwdPtr)[i].backPtr().isAvailable()) {
0359           inputs_.push_back((*packedinputsHandleAsFwdPtr)[i].backPtr());
0360         }
0361       }
0362     } else if (isGen) {
0363       if (geninputsHandleAsFwdPtr->empty()) {
0364         output(iEvent, iSetup);
0365         return;
0366       }
0367       for (size_t i = 0; i < geninputsHandleAsFwdPtr->size(); ++i) {
0368         if ((*geninputsHandleAsFwdPtr)[i].ptr().isAvailable()) {
0369           inputs_.push_back((*geninputsHandleAsFwdPtr)[i].ptr());
0370         } else if ((*geninputsHandleAsFwdPtr)[i].backPtr().isAvailable()) {
0371           inputs_.push_back((*geninputsHandleAsFwdPtr)[i].backPtr());
0372         }
0373       }
0374     } else if (isGenFwdPtr) {
0375       if (geninputsHandleAsFwdPtr->empty()) {
0376         output(iEvent, iSetup);
0377         return;
0378       }
0379       for (size_t i = 0; i < packedgeninputsHandleAsFwdPtr->size(); ++i) {
0380         if ((*packedgeninputsHandleAsFwdPtr)[i].ptr().isAvailable()) {
0381           inputs_.push_back((*packedgeninputsHandleAsFwdPtr)[i].ptr());
0382         } else if ((*packedgeninputsHandleAsFwdPtr)[i].backPtr().isAvailable()) {
0383           inputs_.push_back((*packedgeninputsHandleAsFwdPtr)[i].backPtr());
0384         }
0385       }
0386     } else {
0387       throw cms::Exception("Invalid Jet Inputs")
0388           << "Did not specify appropriate inputs for VirtualJetProducer, Abort!\n";
0389     }
0390   }
0391   LogDebug("VirtualJetProducer") << "Got inputs\n";
0392
0393   // Convert candidates to fastjet::PseudoJets.
0394   // Also correct to Primary Vertex. Will modify fjInputs_
0395   // and use inputs_
0396   fjInputs_.reserve(inputs_.size());
0397   inputTowers();
0398   LogDebug("VirtualJetProducer") << "Inputted towers\n";
0399
0400   // For Pileup subtraction using offset correction:
0401   // Subtract pedestal.
0402   if (doPUOffsetCorr_) {
0403     subtractor_->setDefinition(fjJetDefinition_);
0404     subtractor_->reset(inputs_, fjInputs_, fjJets_);
0405     subtractor_->calculatePedestal(fjInputs_);
0406     subtractor_->subtractPedestal(fjInputs_);
0407     LogDebug("VirtualJetProducer") << "Subtracted pedestal\n";
0408   }
0409   // Run algorithm. Will modify fjJets_ and allocate fjClusterSeq_.
0410   // This will use fjInputs_
0411   runAlgorithm(iEvent, iSetup);
0412
0413   // if ( doPUOffsetCorr_ ) {
0414   //    subtractor_->setAlgorithm(fjClusterSeq_);
0415   // }
0416
0417   LogDebug("VirtualJetProducer") << "Ran algorithm\n";
0418   // For Pileup subtraction using offset correction:
0419   // Now we find jets and need to recalculate their energy,
0420   // mark towers participated in jet,
0421   // remove occupied towers from the list and recalculate mean and sigma
0422   // put the initial towers collection to the jet,
0423   // and subtract from initial towers in jet recalculated mean and sigma of towers
0424   if (doPUOffsetCorr_) {
0425     LogDebug("VirtualJetProducer") << "Do PUOffsetCorr\n";
0426     vector<fastjet::PseudoJet> orphanInput;
0427     subtractor_->calculateOrphanInput(orphanInput);
0428     subtractor_->calculatePedestal(orphanInput);
0429     subtractor_->offsetCorrectJets();
0430   }
0431   // Write the output jets.
0432   // This will (by default) call the member function template
0433   // "writeJets", but can be overridden.
0434   // this will use inputs_
0435   output(iEvent, iSetup);
0436   LogDebug("VirtualJetProducer") << "Wrote jets\n";
0437
0438   // Clear the work vectors so that memory is free for other modules.
0439   // Use the trick of swapping with an empty vector so that the memory
0440   // is actually given back rather than silently kept.
0441   decltype(fjInputs_)().swap(fjInputs_);
0442   decltype(fjJets_)().swap(fjJets_);
0443   decltype(inputs_)().swap(inputs_);
0444
0445   return;
0446 }
0447
0448 //______________________________________________________________________________
0449
0450 void VirtualJetProducer::inputTowers() {
0451   auto inBegin = inputs_.begin(), inEnd = inputs_.end(), i = inBegin;
0452   for (; i != inEnd; ++i) {
0453     auto const& input = **i;
0454     // std::cout << "CaloTowerVI jets " << input->pt() << " " << input->et() << ' '<< input->energy() << ' ' << (isAnomalousTower(input) ? " bad" : " ok") << std::endl;
0455     if (edm::isNotFinite(input.pt()))
0456       continue;
0457     if (input.et() < inputEtMin_)
0458       continue;
0459     if (input.energy() < inputEMin_)
0460       continue;
0461     if (isAnomalousTower(*i))
0462       continue;
0463     // Change by SRR : this is no longer an error nor warning, this can happen with PU mitigation algos.
0464     // Also switch to something more numerically safe. (VI: 10^-42GeV????)
0465     if (input.pt() < 100 * std::numeric_limits<double>::epsilon()) {
0466       continue;
0467     }
0468     if (makeCaloJet(jetTypeE) && doPVCorrection_) {
0469       const CaloTower& tower = dynamic_cast<const CaloTower&>(input);
0470       auto const& ct = tower.p4(vertex_);  // very expensive as computed in eta/phi
0471       fjInputs_.emplace_back(ct.px(), ct.py(), ct.pz(), ct.energy());
0472       fjInputs_.back().set_user_index(i - inBegin);
0473       //std::cout << "tower:" << *tower << '\n';
0474     } else {
0475       /*
0476       if(makePFJet(jetTypeE)) {
0477     reco::PFCandidate& pfc = (reco::PFCandidate&)input;
0478     std::cout << "PF cand:" << pfc << '\n';
0479       }
0480       */
0481       if (!applyWeight_) {
0482         fjInputs_.emplace_back(input.px(), input.py(), input.pz(), input.energy());
0483         fjInputs_.back().set_user_index(i - inBegin);
0484       } else {
0485         if (input_weights_token_.isUninitialized())
0486           throw cms::Exception("InvalidInput")
0487               << "applyWeight set to True, but no weights given in VirtualJetProducer\n";
0488         float w = weights_[*i];
0489         if (w > 0) {
0490           fjInputs_.emplace_back(input.px() * w, input.py() * w, input.pz() * w, input.energy() * w);
0491           fjInputs_.back().set_user_index(i - inBegin);
0492         }
0493       }
0494     }
0495   }
0496
0497   if (restrictInputs_ && fjInputs_.size() > maxInputs_) {
0498     reco::helper::GreaterByPtPseudoJet pTComparator;
0499     std::sort(fjInputs_.begin(), fjInputs_.end(), pTComparator);
0500     fjInputs_.resize(maxInputs_);
0501     edm::LogWarning("JetRecoTooManyEntries")
0502         << "Too many inputs in the event, limiting to first " << maxInputs_ << ". Output is suspect.";
0503   }
0504 }
0505
0506 //______________________________________________________________________________
0507 bool VirtualJetProducer::isAnomalousTower(reco::CandidatePtr input) {
0508   if (!makeCaloJet(jetTypeE))
0509     return false;
0510   else
0511     return (*anomalousTowerDef_)(*input);
0512 }
0513
0514 //------------------------------------------------------------------------------
0515 // This is pure virtual.
0516 //______________________________________________________________________________
0517 // void VirtualJetProducer::runAlgorithm( edm::Event & iEvent, edm::EventSetup const& iSetup,
0518 //                                        std::vector<edm::Ptr<reco::Candidate> > const & inputs_);
0519
0520 //______________________________________________________________________________
0521 void VirtualJetProducer::copyConstituents(const vector<fastjet::PseudoJet>& fjConstituents, reco::Jet* jet) {
0522   for (unsigned int i = 0; i < fjConstituents.size(); ++i) {
0523     int index = fjConstituents[i].user_index();
0524     if (index >= 0 && static_cast<unsigned int>(index) < inputs_.size())
0525       jet->addDaughter(inputs_[index]);
0526   }
0527 }
0528
0529 //______________________________________________________________________________
0530 vector<reco::CandidatePtr> VirtualJetProducer::getConstituents(const vector<fastjet::PseudoJet>& fjConstituents) {
0531   vector<reco::CandidatePtr> result;
0532   result.reserve(fjConstituents.size() / 2);
0533   for (unsigned int i = 0; i < fjConstituents.size(); i++) {
0534     auto index = fjConstituents[i].user_index();
0535     if (index >= 0 && static_cast<unsigned int>(index) < inputs_.size()) {
0536       result.emplace_back(inputs_[index]);
0537     }
0538   }
0539   return result;
0540 }
0541
0542 //_____________________________________________________________________________
0543
0544 void VirtualJetProducer::output(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0545   // Write jets and constitutents. Will use fjJets_, inputs_
0546   // and fjClusterSeq_
0547
0548   if (writeCompound_) {
0549     // Write jets and subjets
0550     switch (jetTypeE) {
0551       case JetType::CaloJet:
0552         writeCompoundJets<reco::CaloJet>(iEvent, iSetup);
0553         break;
0554       case JetType::PFJet:
0555         writeCompoundJets<reco::PFJet>(iEvent, iSetup);
0556         break;
0557       case JetType::GenJet:
0558         writeCompoundJets<reco::GenJet>(iEvent, iSetup);
0559         break;
0560       case JetType::BasicJet:
0561         writeCompoundJets<reco::BasicJet>(iEvent, iSetup);
0562         break;
0563       default:
0564         throw cms::Exception("InvalidInput") << "invalid jet type in CompoundJetProducer\n";
0565         break;
0566     };
0567   } else if (writeJetsWithConst_) {
0568     // Write jets and new constituents.
0569     writeJetsWithConstituents<reco::PFJet>(iEvent, iSetup);
0570   } else {
0571     switch (jetTypeE) {
0572       case JetType::CaloJet:
0573         writeJets<reco::CaloJet>(iEvent, iSetup);
0574         break;
0575       case JetType::PFJet:
0576         writeJets<reco::PFJet>(iEvent, iSetup);
0577         break;
0578       case JetType::GenJet:
0579         writeJets<reco::GenJet>(iEvent, iSetup);
0580         break;
0581       case JetType::TrackJet:
0582         writeJets<reco::TrackJet>(iEvent, iSetup);
0583         break;
0584       case JetType::PFClusterJet:
0585         writeJets<reco::PFClusterJet>(iEvent, iSetup);
0586         break;
0587       case JetType::BasicJet:
0588         writeJets<reco::BasicJet>(iEvent, iSetup);
0589         break;
0590       default:
0591         throw cms::Exception("InvalidInput") << "invalid jet type in VirtualJetProducer\n";
0592         break;
0593     };
0594   }
0595 }
0596
0597 namespace {
0598   template <typename T>
0599   struct Area {
0600     static float get(T const&) { return 0; }
0601   };
0602
0603   template <>
0604   struct Area<reco::CaloJet> {
0605     static float get(reco::CaloJet const& jet) { return jet.getSpecific().mTowersArea; }
0606   };
0607 }  // namespace
0608
0609 template <typename T>
0610 void VirtualJetProducer::writeJets(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0611   // std::cout << "writeJets " << typeid(T).name()
0612   //          << (doRhoFastjet_ ? " doRhoFastjet " : "")
0613   //          << (doAreaFastjet_ ? " doAreaFastjet " : "")
0614   //          << (doAreaDiskApprox_ ? " doAreaDiskApprox " : "")
0615   //          << std::endl;
0616
0617   if (doRhoFastjet_) {
0618     // declare jet collection without the two jets,
0619     // for unbiased background estimation.
0620     std::vector<fastjet::PseudoJet> fjexcluded_jets;
0621     fjexcluded_jets = fjJets_;
0622
0623     if (fjexcluded_jets.size() > 2)
0624       fjexcluded_jets.resize(nExclude_);
0625
0626     if (doFastJetNonUniform_) {
0627       auto rhos = std::make_unique<std::vector<double>>();
0628       auto sigmas = std::make_unique<std::vector<double>>();
0629       int nEta = puCenters_.size();
0630       rhos->reserve(nEta);
0631       sigmas->reserve(nEta);
0632       fastjet::ClusterSequenceAreaBase const* clusterSequenceWithArea =
0633           dynamic_cast<fastjet::ClusterSequenceAreaBase const*>(fjClusterSeq_.get());
0634
0635       if (clusterSequenceWithArea == nullptr) {
0636         if (!fjJets_.empty()) {
0637           throw cms::Exception("LogicError") << "fjClusterSeq is not initialized while inputs are present\n ";
0638         }
0639       } else {
0640         for (int ie = 0; ie < nEta; ++ie) {
0641           double eta = puCenters_[ie];
0642           double etamin = eta - puWidth_;
0643           double etamax = eta + puWidth_;
0644           fastjet::RangeDefinition range_rho(etamin, etamax);
0645           fastjet::BackgroundEstimator bkgestim(*clusterSequenceWithArea, range_rho);
0646           bkgestim.set_excluded_jets(fjexcluded_jets);
0647           rhos->push_back(bkgestim.rho());
0648           sigmas->push_back(bkgestim.sigma());
0649         }
0650       }
0651       iEvent.put(std::move(rhos), "rhos");
0652       iEvent.put(std::move(sigmas), "sigmas");
0653     } else {
0654       auto rho = std::make_unique<double>(0.0);
0655       auto sigma = std::make_unique<double>(0.0);
0656       double mean_area = 0;
0657
0658       fastjet::ClusterSequenceAreaBase const* clusterSequenceWithArea =
0659           dynamic_cast<fastjet::ClusterSequenceAreaBase const*>(fjClusterSeq_.get());
0660       if (clusterSequenceWithArea == nullptr) {
0661         if (!fjJets_.empty()) {
0662           throw cms::Exception("LogicError") << "fjClusterSeq is not initialized while inputs are present\n ";
0663         }
0664       } else {
0665         clusterSequenceWithArea->get_median_rho_and_sigma(*fjSelector_, false, *rho, *sigma, mean_area);
0666         if ((*rho < 0) || (edm::isNotFinite(*rho))) {
0667           edm::LogError("BadRho") << "rho value is " << *rho << " area:" << mean_area
0668                                   << " and n_empty_jets: " << clusterSequenceWithArea->n_empty_jets(*fjSelector_)
0669                                   << " with range " << fjSelector_->description() << ". Setting rho to rezo.";
0670           *rho = 0;
0671         }
0672       }
0673       iEvent.put(std::move(rho), "rho");
0674       iEvent.put(std::move(sigma), "sigma");
0675     }
0676   }  // doRhoFastjet_
0677
0678   // produce output jet collection
0679
0680   using namespace reco;
0681
0682   // allocate fjJets_.size() Ts in vector
0683   auto jets = std::make_unique<std::vector<T>>(fjJets_.size());
0684
0685   if (!fjJets_.empty()) {
0686     // Distance between jet centers and overlap area -- for disk-based area calculation
0687     using RIJ = std::pair<double, double>;
0688     std::vector<RIJ> rijStorage(fjJets_.size() == 1 ? 1 : fjJets_.size() * (fjJets_.size() >> 1));
0689     std::vector<RIJ*> rij(fjJets_.size());
0690     unsigned int k = 0;
0691     for (unsigned int ijet = 0; ijet < fjJets_.size(); ++ijet) {
0692       rij[ijet] = &rijStorage[k];
0693       k += ijet;
0694     }
0695
0696     float etaJ[fjJets_.size()], phiJ[fjJets_.size()];
0697
0698     auto orParam_ = 1. / rParam_;
0699     // fill jets
0700     [[maybe_unused]] const CaloGeometry* pGeometry = nullptr;
0701     [[maybe_unused]] const HcalTopology* pTopology = nullptr;
0702     if constexpr (std::is_same_v<T, reco::CaloJet>) {
0703       pGeometry = &getGeometry(iSetup);
0704       pTopology = &getTopology(iSetup);
0705     }
0706     for (unsigned int ijet = 0; ijet < fjJets_.size(); ++ijet) {
0707       auto& jet = (*jets)[ijet];
0708
0709       // get the fastjet jet
0710       const fastjet::PseudoJet& fjJet = fjJets_[ijet];
0711       // get the constituents from fastjet
0712       std::vector<fastjet::PseudoJet> const& fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
0713       // convert them to CandidatePtr vector
0714       std::vector<CandidatePtr> const& constituents = getConstituents(fjConstituents);
0715
0716       // write the specifics to the jet (simultaneously sets 4-vector, vertex).
0717       // These are overridden functions that will call the appropriate
0718       // specific allocator.
0719       if ((applyWeight_) && (makePFJet(jetTypeE)))
0720         writeSpecific(dynamic_cast<reco::PFJet&>(jet),
0721                       Particle::LorentzVector(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E()),
0722                       vertex_,
0723                       constituents,
0724                       &weights_);
0725       else {
0726         if constexpr (std::is_same_v<T, reco::CaloJet>) {
0727           writeSpecific(jet,
0728                         Particle::LorentzVector(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E()),
0729                         vertex_,
0730                         constituents,
0731                         *pGeometry,
0732                         *pTopology);
0733         } else {
0734           writeSpecific(
0735               jet, Particle::LorentzVector(fjJet.px(), fjJet.py(), fjJet.pz(), fjJet.E()), vertex_, constituents);
0736         }
0737       }
0738       phiJ[ijet] = jet.phi();
0739       etaJ[ijet] = jet.eta();
0740       jet.setIsWeighted(applyWeight_);
0741     }
0742
0743     // calcuate the jet area
0744     for (unsigned int ijet = 0; ijet < fjJets_.size(); ++ijet) {
0745       // calcuate the jet area
0746       double jetArea = 0.0;
0747       // get the fastjet jet
0748       const auto& fjJet = fjJets_[ijet];
0749       if (doAreaFastjet_ && fjJet.has_area()) {
0750         jetArea = fjJet.area();
0751       } else if (doAreaDiskApprox_) {
0752         // Here it is assumed that fjJets_ is in decreasing order of pT,
0753         // which should happen in FastjetJetProducer::runAlgorithm()
0754         jetArea = M_PI;
0755         RIJ* distance = rij[ijet];
0756         for (unsigned jJet = 0; jJet < ijet; ++jJet) {
0757           [[clang::suppress]] distance[jJet].first =
0758               std::sqrt(reco::deltaR2(etaJ[ijet], phiJ[ijet], etaJ[jJet], phiJ[jJet])) * orParam_;
0759           distance[jJet].second = reco::helper::VirtualJetProducerHelper::intersection(distance[jJet].first);
0760           jetArea -= distance[jJet].second;
0761           for (unsigned kJet = 0; kJet < jJet; ++kJet) {
0762             jetArea += reco::helper::VirtualJetProducerHelper::intersection(distance[jJet].first,
0763                                                                             distance[kJet].first,
0764                                                                             rij[jJet][kJet].first,
0765                                                                             distance[jJet].second,
0766                                                                             distance[kJet].second,
0767                                                                             rij[jJet][kJet].second);
0768           }  // end loop over harder jets
0769         }  // end loop over harder jets
0770         jetArea *= (rParam_ * rParam_);
0771       }
0772       auto& jet = (*jets)[ijet];
0773       jet.setJetArea(jetArea);
0774
0775       if (doPUOffsetCorr_) {
0776         jet.setPileup(subtractor_->getPileUpEnergy(ijet));
0777       } else {
0778         jet.setPileup(0.0);
0779       }
0780
0781       // std::cout << "area " << ijet << " " << jetArea << " " << Area<T>::get(jet) << std::endl;
0782       // std::cout << "JetVI " << ijet << ' '<< jet.pt() << " " << jet.et() << ' '<< jet.energy() << ' '<< jet.mass() << std::endl;
0783     }
0784   }
0785   // put the jets in the collection
0786   iEvent.put(std::move(jets), jetCollInstanceName_);
0787 }
0788
0789 /// function template to write out the outputs
0790 template <class T>
0791 void VirtualJetProducer::writeCompoundJets(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0792   if (verbosity_ >= 1) {
0793     std::cout << "<VirtualJetProducer::writeCompoundJets (moduleLabel = " << moduleLabel_ << ")>:" << std::endl;
0794   }
0795
0796   // get a list of output jets
0797   auto jetCollection = std::make_unique<reco::BasicJetCollection>();
0798   // get a list of output subjets
0799   auto subjetCollection = std::make_unique<std::vector<T>>();
0800
0801   // This will store the handle for the subjets after we write them
0802   edm::OrphanHandle<std::vector<T>> subjetHandleAfterPut;
0803   // this is the mapping of subjet to hard jet
0804   std::vector<std::vector<int>> indices;
0805   // this is the list of hardjet 4-momenta
0806   std::vector<math::XYZTLorentzVector> p4_hardJets;
0807   // this is the hardjet areas
0808   std::vector<double> area_hardJets;
0809
0810   // Loop over the hard jets
0811   std::vector<fastjet::PseudoJet>::const_iterator it = fjJets_.begin(), iEnd = fjJets_.end(), iBegin = fjJets_.begin();
0812   indices.resize(fjJets_.size());
0813   for (; it != iEnd; ++it) {
0814     fastjet::PseudoJet const& localJet = *it;
0815     unsigned int jetIndex = it - iBegin;
0816     // Get the 4-vector for the hard jet
0817     p4_hardJets.push_back(math::XYZTLorentzVector(localJet.px(), localJet.py(), localJet.pz(), localJet.e()));
0818     double localJetArea = 0.0;
0819     if (doAreaFastjet_ && localJet.has_area()) {
0820       localJetArea = localJet.area();
0821     }
0822     area_hardJets.push_back(localJetArea);
0823
0824     // create the subjet list
0825     std::vector<fastjet::PseudoJet> constituents;
0826     if (it->has_pieces()) {
0827       constituents = it->pieces();
0828     } else if (it->has_constituents()) {
0829       constituents = it->constituents();
0830     }
0831
0832     std::vector<fastjet::PseudoJet>::const_iterator itSubJetBegin = constituents.begin(), itSubJet = itSubJetBegin,
0833                                                     itSubJetEnd = constituents.end();
0834     for (; itSubJet != itSubJetEnd; ++itSubJet) {
0835       fastjet::PseudoJet const& subjet = *itSubJet;
0836       if (verbosity_ >= 1) {
0837         std::cout << "subjet #" << (itSubJet - itSubJetBegin) << ": Pt = " << subjet.pt() << ", eta = " << subjet.eta()
0838                   << ", phi = " << subjet.phi() << ", mass = " << subjet.m()
0839                   << " (#constituents = " << subjet.constituents().size() << ")" << std::endl;
0840         std::vector<fastjet::PseudoJet> subjet_constituents = subjet.constituents();
0841         int idx_constituent = 0;
0842         for (std::vector<fastjet::PseudoJet>::const_iterator constituent = subjet_constituents.begin();
0843              constituent != subjet_constituents.end();
0844              ++constituent) {
0845           if (constituent->pt() < 1.e-3)
0846             continue;  // CV: skip ghosts
0847           std::cout << "  constituent #" << idx_constituent << ": Pt = " << constituent->pt()
0848                     << ", eta = " << constituent->eta() << ", phi = " << constituent->phi() << ","
0849                     << " mass = " << constituent->m() << std::endl;
0850           ++idx_constituent;
0851         }
0852       }
0853
0854       if (verbosity_ >= 1) {
0855         std::cout << "subjet #" << (itSubJet - itSubJetBegin) << ": Pt = " << subjet.pt() << ", eta = " << subjet.eta()
0856                   << ", phi = " << subjet.phi() << ", mass = " << subjet.m()
0857                   << " (#constituents = " << subjet.constituents().size() << ")" << std::endl;
0858         std::vector<fastjet::PseudoJet> subjet_constituents = subjet.constituents();
0859         int idx_constituent = 0;
0860         for (std::vector<fastjet::PseudoJet>::const_iterator constituent = subjet_constituents.begin();
0861              constituent != subjet_constituents.end();
0862              ++constituent) {
0863           if (constituent->pt() < 1.e-3)
0864             continue;  // CV: skip ghosts
0865           std::cout << " constituent #" << idx_constituent << ": Pt = " << constituent->pt()
0866                     << ", eta = " << constituent->eta() << ", phi = " << constituent->phi() << ","
0867                     << " mass = " << constituent->m() << std::endl;
0868           ++idx_constituent;
0869         }
0870       }
0871
0872       math::XYZTLorentzVector p4Subjet(subjet.px(), subjet.py(), subjet.pz(), subjet.e());
0873       reco::Particle::Point point(0, 0, 0);
0874
0875       // This will hold ptr's to the subjets
0876       std::vector<reco::CandidatePtr> subjetConstituents;
0877
0878       // Get the transient subjet constituents from fastjet
0879       std::vector<fastjet::PseudoJet> subjetFastjetConstituents = subjet.constituents();
0880       std::vector<reco::CandidatePtr> constituents = getConstituents(subjetFastjetConstituents);
0881
0882       indices[jetIndex].push_back(subjetCollection->size());
0883
0884       // Add the concrete subjet type to the subjet list to write to event record
0885       subjetCollection->push_back(*std::make_unique<T>());
0886       auto& jet = subjetCollection->back();
0887       if ((applyWeight_) && (makePFJet(jetTypeE)))
0888         reco::writeSpecific(dynamic_cast<reco::PFJet&>(jet), p4Subjet, point, constituents, &weights_);
0889       else if constexpr (std::is_same_v<T, reco::CaloJet>) {
0890         reco::writeSpecific(
0891             jet, p4Subjet, point, constituents, iSetup.getData(geometry_token_), iSetup.getData(topology_token_));
0892       } else {
0893         reco::writeSpecific(jet, p4Subjet, point, constituents);
0894       }
0895       jet.setIsWeighted(applyWeight_);
0896       double subjetArea = 0.0;
0897       if (doAreaFastjet_ && itSubJet->has_area()) {
0898         subjetArea = itSubJet->area();
0899       }
0900       jet.setJetArea(subjetArea);
0901     }
0902   }
0903
0904   // put subjets into event record
0905   subjetHandleAfterPut = iEvent.put(std::move(subjetCollection), jetCollInstanceName_);
0906
0907   // Now create the hard jets with ptr's to the subjets as constituents
0908   std::vector<math::XYZTLorentzVector>::const_iterator ip4 = p4_hardJets.begin(), ip4Begin = p4_hardJets.begin(),
0909                                                        ip4End = p4_hardJets.end();
0910
0911   for (; ip4 != ip4End; ++ip4) {
0912     int p4_index = ip4 - ip4Begin;
0913     std::vector<int>& ind = indices[p4_index];
0914     std::vector<reco::CandidatePtr> i_hardJetConstituents;
0915     // Add the subjets to the hard jet
0916     for (std::vector<int>::const_iterator isub = ind.begin(); isub != ind.end(); ++isub) {
0917       reco::CandidatePtr candPtr(subjetHandleAfterPut, *isub, false);
0918       i_hardJetConstituents.push_back(candPtr);
0919     }
0920     reco::Particle::Point point(0, 0, 0);
0921     reco::BasicJet toput(*ip4, point, i_hardJetConstituents);
0922     toput.setJetArea(area_hardJets[ip4 - ip4Begin]);
0923     jetCollection->push_back(toput);
0924   }
0925
0926   // put hard jets into event record
0927   // Store the Orphan handle for adding HTT information
0928   edm::OrphanHandle<reco::BasicJetCollection> oh = iEvent.put(std::move(jetCollection));
0929
0930   if (fromHTTTopJetProducer_) {
0931     addHTTTopJetTagInfoCollection(iEvent, iSetup, oh);
0932   }
0933 }
0934
0935 /// function template to write out the outputs
0936 template <class T>
0937 void VirtualJetProducer::writeJetsWithConstituents(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0938   if (verbosity_ >= 1) {
0939     std::cout << "<VirtualJetProducer::writeJetsWithConstituents (moduleLabel = " << moduleLabel_ << ")>:" << std::endl;
0940   }
0941
0942   // get a list of output jets  MV: make this compatible with template
0943   auto jetCollection = std::make_unique<reco::PFJetCollection>();
0944
0945   // this is the mapping of jet to constituents
0946   std::vector<std::vector<int>> indices;
0947   // this is the list of jet 4-momenta
0948   std::vector<math::XYZTLorentzVector> p4_Jets;
0949   // this is the jet areas
0950   std::vector<double> area_Jets;
0951
0952   // get a list of output constituents
0953   auto constituentCollection = std::make_unique<reco::PFCandidateCollection>();
0954
0955   // This will store the handle for the constituents after we write them
0956   edm::OrphanHandle<reco::PFCandidateCollection> constituentHandleAfterPut;
0957
0958   // Loop over the jets and extract constituents
0959   std::vector<fastjet::PseudoJet> constituentsSub;
0960   std::vector<fastjet::PseudoJet>::const_iterator it = fjJets_.begin(), iEnd = fjJets_.end(), iBegin = fjJets_.begin();
0961   indices.resize(fjJets_.size());
0962
0963   for (; it != iEnd; ++it) {
0964     fastjet::PseudoJet const& localJet = *it;
0965     unsigned int jetIndex = it - iBegin;
0966     // Get the 4-vector for the hard jet
0967     p4_Jets.push_back(math::XYZTLorentzVector(localJet.px(), localJet.py(), localJet.pz(), localJet.e()));
0968     double localJetArea = 0.0;
0969     if (doAreaFastjet_ && localJet.has_area()) {
0970       localJetArea = localJet.area();
0971     }
0972     area_Jets.push_back(localJetArea);
0973
0974     // create the constituent list
0975     std::vector<fastjet::PseudoJet> constituents, ghosts;
0976     if (it->has_pieces())
0977       constituents = it->pieces();
0978     else if (it->has_constituents())
0979       fastjet::SelectorIsPureGhost().sift(it->constituents(), ghosts, constituents);  //filter out ghosts
0980
0981     //loop over constituents of jet (can be subjets or normal constituents)
0982     indices[jetIndex].reserve(constituents.size());
0983     constituentsSub.reserve(constituentsSub.size() + constituents.size());
0984     for (fastjet::PseudoJet const& constit : constituents) {
0985       indices[jetIndex].push_back(constituentsSub.size());
0986       constituentsSub.push_back(constit);
0987     }
0988   }
0989
0990   //Loop over constituents and store in the event
0991   static const reco::PFCandidate dummySinceTranslateIsNotStatic;
0992   for (fastjet::PseudoJet const& constit : constituentsSub) {
0993     auto orig = inputs_[constit.user_index()];
0994     auto id = dummySinceTranslateIsNotStatic.translatePdgIdToType(orig->pdgId());
0995     reco::PFCandidate pCand(reco::PFCandidate(orig->charge(), orig->p4(), id));
0996     math::XYZTLorentzVector pVec;
0997     pVec.SetPxPyPzE(constit.px(), constit.py(), constit.pz(), constit.e());
0998     pCand.setP4(pVec);
0999     pCand.setSourceCandidatePtr(orig->sourceCandidatePtr(0));
1000     constituentCollection->push_back(pCand);
1001   }
1002   // put constituents into event record
1003   constituentHandleAfterPut = iEvent.put(std::move(constituentCollection), jetCollInstanceName_);
1004
1005   // Now create the jets with ptr's to the constituents
1006   std::vector<math::XYZTLorentzVector>::const_iterator ip4 = p4_Jets.begin(), ip4Begin = p4_Jets.begin(),
1007                                                        ip4End = p4_Jets.end();
1008
1009   for (; ip4 != ip4End; ++ip4) {
1010     int p4_index = ip4 - ip4Begin;
1011     std::vector<int>& ind = indices[p4_index];
1012     std::vector<reco::CandidatePtr> i_jetConstituents;
1013     // Add the constituents to the jet
1014     for (std::vector<int>::const_iterator iconst = ind.begin(); iconst != ind.end(); ++iconst) {
1015       reco::CandidatePtr candPtr(constituentHandleAfterPut, *iconst, false);
1016       i_jetConstituents.push_back(candPtr);
1017     }
1018     if (!i_jetConstituents.empty()) {  //only keep jets which have constituents after subtraction
1019       reco::Particle::Point point(0, 0, 0);
1020       reco::PFJet jet;
1021       reco::writeSpecific(jet, *ip4, point, i_jetConstituents);
1022       jet.setJetArea(area_Jets[ip4 - ip4Begin]);
1023       jetCollection->emplace_back(jet);
1024     }
1025   }
1026
1027   // put jets into event record
1028   iEvent.put(std::move(jetCollection));
1029 }
1030
1031 CaloGeometry const& VirtualJetProducer::getGeometry(edm::EventSetup const& iSetup) const {
1032   return iSetup.getData(geometry_token_);
1033 }
1034
1035 HcalTopology const& VirtualJetProducer::getTopology(edm::EventSetup const& iSetup) const {
1036   return iSetup.getData(topology_token_);
1037 }
1038
1039 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
1040 void VirtualJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1041   edm::ParameterSetDescription desc;
1042   fillDescriptionsFromVirtualJetProducer(desc);
1043   desc.add<string>("jetCollInstanceName", "");
1044
1045   // addDefault must be used here instead of add unless
1046   // all the classes that inherit from this class redefine
1047   // the fillDescriptions function. Otherwise, the autogenerated
1048   // cfi filenames are the same and conflict.
1049   descriptions.addDefault(desc);
1050 }
1051
1052 void VirtualJetProducer::fillDescriptionsFromVirtualJetProducer(edm::ParameterSetDescription& desc) {
1053   desc.add<edm::InputTag>("src", edm::InputTag("particleFlow"));
1054   desc.add<edm::InputTag>("srcPVs", edm::InputTag(""));
1055   desc.add<string>("jetType", "PFJet");
1056   desc.add<string>("jetAlgorithm", "AntiKt");
1057   desc.add<double>("rParam", 0.4);
1058   desc.add<double>("inputEtMin", 0.0);
1059   desc.add<double>("inputEMin", 0.0);
1060   desc.add<double>("jetPtMin", 5.);
1061   desc.add<bool>("doPVCorrection", false);
1062   desc.add<bool>("doAreaFastjet", false);
1063   desc.add<bool>("doRhoFastjet", false);
1064   desc.add<bool>("doPUOffsetCorr", false);
1065   desc.add<double>("puPtMin", 10.);
1066   desc.add<double>("nSigmaPU", 1.0);
1067   desc.add<double>("radiusPU", 0.5);
1068   desc.add<string>("subtractorName", "");
1069   desc.add<bool>("useExplicitGhosts", false);
1070   desc.add<bool>("doAreaDiskApprox", false);
1071   desc.add<double>("voronoiRfact", -0.9);
1072   desc.add<double>("Rho_EtaMax", 4.4);
1073   desc.add<double>("Ghost_EtaMax", 5.);
1074   desc.add<int>("Active_Area_Repeats", 1);
1075   desc.add<double>("GhostArea", 0.01);
1076   desc.add<bool>("restrictInputs", false);
1077   desc.add<unsigned int>("maxInputs", 1);
1078   desc.add<bool>("writeCompound", false);
1079   desc.add<bool>("writeJetsWithConst", false);
1080   desc.add<bool>("doFastJetNonUniform", false);
1081   desc.add<bool>("useDeterministicSeed", false);
1082   desc.add<unsigned int>("minSeed", 14327);
1083   desc.add<int>("verbosity", 0);
1084   desc.add<double>("puWidth", 0.);
1085   desc.add<unsigned int>("nExclude", 0);
1086   desc.add<unsigned int>("maxBadEcalCells", 9999999);
1087   desc.add<unsigned int>("maxBadHcalCells", 9999999);
1088   desc.add<unsigned int>("maxProblematicEcalCells", 9999999);
1089   desc.add<unsigned int>("maxProblematicHcalCells", 9999999);
1090   desc.add<unsigned int>("maxRecoveredEcalCells", 9999999);
1091   desc.add<unsigned int>("maxRecoveredHcalCells", 9999999);
1092   vector<double> puCentersDefault;
1093   desc.add<vector<double>>("puCenters", puCentersDefault);
1094   desc.add<bool>("applyWeight", false);
1095   desc.add<edm::InputTag>("srcWeights", edm::InputTag(""));
1096   desc.add<double>("minimumTowersFraction", 0.);
1097 }