#include <RecoJets/FFTJetProducers/plugins/FFTJetPileupProcessor.cc>

Inheritance diagram for FFTJetPileupProcessor:

Public Member Functions
	FFTJetPileupProcessor (const edm::ParameterSet &)

	~FFTJetPileupProcessor () override

Public Member Functions inherited from fftjetcms::FFTJetInterface
	~FFTJetInterface () override

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

SerialTaskQueue *	globalLuminosityBlocksQueue ()

SerialTaskQueue *	globalRunsQueue ()

ModuleDescription const &	moduleDescription () const

	~EDProducer () override

Public Member Functions inherited from edm::ProducerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

std::vector< edm::ProductResolverIndex > const &	indiciesForPutProducts (BranchType iBranchType) const

	ProducerBase ()

std::vector< edm::ProductResolverIndex > const &	putTokenIndexToProductResolverIndex () const

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

std::function< void(BranchDescription const &)>	registrationCallback () const
	used by the fwk to register list of products More...

void	resolvePutIndicies (BranchType iBranchType, ModuleToResolverIndicies const &iIndicies, std::string const &moduleLabel)

	~ProducerBase () noexcept(false) override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Protected Member Functions
void	beginJob () override

void	endJob () override

void	produce (edm::Event &, const edm::EventSetup &) override

Protected Member Functions inherited from fftjetcms::FFTJetInterface
template<class Ptr >
void	checkConfig (const Ptr &ptr, const char *message)

void	discretizeEnergyFlow ()

	FFTJetInterface (const edm::ParameterSet &)

double	getEventScale () const

void	loadInputCollection (const edm::Event &)

bool	storeInSinglePrecision () const

const reco::Particle::Point &	vertexUsed () const

Protected Member Functions inherited from edm::ProducerBase
ProducesCollector	producesCollector ()

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Private Member Functions
void	buildKernelConvolver (const edm::ParameterSet &)

	FFTJetPileupProcessor ()=delete

	FFTJetPileupProcessor (const FFTJetPileupProcessor &)=delete

void	loadFlatteningFactors (const edm::EventSetup &iSetup)

void	mixExtraGrid ()

FFTJetPileupProcessor &	operator= (const FFTJetPileupProcessor &)=delete

Private Attributes
std::unique_ptr< fftjet::Grid2d< fftjetcms::Real > >	convolvedFlow

std::unique_ptr< fftjet::AbsConvolverBase< Real > >	convolver

unsigned	convolverMaxBin

unsigned	convolverMinBin

unsigned	currentFileNum

std::unique_ptr< MyFFTEngine >	engine

std::vector< double >	etaFlatteningFactors

std::vector< std::string >	externalGridFiles

double	externalGridMaxEnergy

std::unique_ptr< fftjet::EquidistantInLogSpace >	filterScales

std::string	flatteningTableCategory

std::string	flatteningTableName

std::string	flatteningTableRecord

std::ifstream	gridStream

std::unique_ptr< fftjet::AbsFrequencyKernel >	kernel2d

bool	loadFlatteningFactorsFromDB

unsigned	nPercentiles

std::vector< double >	percentileData

double	pileupEtaPhiArea

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

Public Types inherited from edm::ProducerBase
using	ModuleToResolverIndicies = std::unordered_multimap< std::string, std::tuple< edm::TypeID const , const char , edm::ProductResolverIndex >>

typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

static bool	wantsGlobalLuminosityBlocks ()

static bool	wantsGlobalRuns ()

static bool	wantsStreamLuminosityBlocks ()

static bool	wantsStreamRuns ()

Protected Attributes inherited from fftjetcms::FFTJetInterface
const AnomalousTower	anomalous

std::vector< unsigned >	candidateIndex

const bool	doPVCorrection

std::unique_ptr< fftjet::Grid2d< fftjetcms::Real > >	energyFlow

const std::vector< double >	etaDependentMagnutideFactors

std::vector< fftjetcms::VectorLike >	eventData

edm::Handle< reco::CandidateView >	inputCollection

const edm::InputTag	inputLabel

const JetType	jetType

const std::string	outputLabel

const edm::InputTag	srcPVs

Detailed Description

Description: Runs FFTJet multiscale pileup filtering code and saves the results

Implementation: [Notes on implementation]

Definition at line 52 of file FFTJetPileupProcessor.cc.

Constructor & Destructor Documentation

FFTJetPileupProcessor::FFTJetPileupProcessor ( const edm::ParameterSet & ps )

explicit

Definition at line 123 of file FFTJetPileupProcessor.cc.

References buildKernelConvolver(), fftjetcms::FFTJetInterface::checkConfig(), convolvedFlow, fftjetcms::FFTJetInterface::energyFlow, etaFlatteningFactors, Exception, fftjetcms::fftjet_Grid2d_parser(), filterScales, edm::ParameterSet::getParameter(), fftjetcommon_cfi::maxScale, fftjetcommon_cfi::minScale, nPercentiles, fftjetcommon_cfi::nScales, fftjetcms::FFTJetInterface::outputLabel, percentileData, and mitigatedMETSequence_cff::U.

     : FFTJetInterface(ps),
       etaFlatteningFactors(ps.getParameter<std::vector<double>>("etaFlatteningFactors")),
       nPercentiles(ps.getParameter<unsigned>("nPercentiles")),
       convolverMinBin(ps.getParameter<unsigned>("convolverMinBin")),
       convolverMaxBin(ps.getParameter<unsigned>("convolverMaxBin")),
       pileupEtaPhiArea(ps.getParameter<double>("pileupEtaPhiArea")),
       externalGridFiles(ps.getParameter<std::vector<std::string>>("externalGridFiles")),
       externalGridMaxEnergy(ps.getParameter<double>("externalGridMaxEnergy")),
       currentFileNum(externalGridFiles.size() + 1U),
       flatteningTableRecord(ps.getParameter<std::string>("flatteningTableRecord")),
       flatteningTableName(ps.getParameter<std::string>("flatteningTableName")),
       flatteningTableCategory(ps.getParameter<std::string>("flatteningTableCategory")),
       loadFlatteningFactorsFromDB(ps.getParameter<bool>("loadFlatteningFactorsFromDB")) {
   // Build the discretization grid
   energyFlow = fftjet_Grid2d_parser(ps.getParameter<edm::ParameterSet>("GridConfiguration"));
   checkConfig(energyFlow, "invalid discretization grid");
 
   // Copy of the grid which will be used for convolutions
   convolvedFlow = std::unique_ptr<fftjet::Grid2d<fftjetcms::Real>>(new fftjet::Grid2d<fftjetcms::Real>(*energyFlow));
 
   // Make sure the size of flattening factors is appropriate
   if (!etaFlatteningFactors.empty()) {
     if (etaFlatteningFactors.size() != convolvedFlow->nEta())
       throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor constructor:"
                                                  " number of elements in the \"etaFlatteningFactors\""
                                                  " vector is inconsistent with the discretization grid binning"
                                               << std::endl;
   }
 
   // Build the FFT engine(s), pattern recognition kernel(s),
   // and the kernel convolver
   buildKernelConvolver(ps);
 
   // Build the set of pattern recognition scales
   const unsigned nScales = ps.getParameter<unsigned>("nScales");
   const double minScale = ps.getParameter<double>("minScale");
   const double maxScale = ps.getParameter<double>("maxScale");
   if (minScale <= 0.0 || maxScale < minScale || nScales == 0U)
     throw cms::Exception("FFTJetBadConfig") << "invalid filter scales" << std::endl;
 
   filterScales =
       std::unique_ptr<fftjet::EquidistantInLogSpace>(new fftjet::EquidistantInLogSpace(minScale, maxScale, nScales));
 
   percentileData.resize(nScales * nPercentiles);
 
   produces<reco::DiscretizedEnergyFlow>(outputLabel);
   produces<std::pair<double, double>>(outputLabel);
 }

FFTJetPileupProcessor::~FFTJetPileupProcessor ( )

override

Definition at line 200 of file FFTJetPileupProcessor.cc.

200 {}

FFTJetPileupProcessor::FFTJetPileupProcessor ( )

privatedelete

FFTJetPileupProcessor::FFTJetPileupProcessor ( const FFTJetPileupProcessor & )

privatedelete

Member Function Documentation

void FFTJetPileupProcessor::beginJob ( void )

overrideprotectedvirtual

Reimplemented from edm::EDProducer.

Definition at line 322 of file FFTJetPileupProcessor.cc.

322 {}

void FFTJetPileupProcessor::buildKernelConvolver ( const edm::ParameterSet & ps )

private

Definition at line 173 of file FFTJetPileupProcessor.cc.

References convolver, convolverMaxBin, convolverMinBin, fftjetcms::FFTJetInterface::energyFlow, engine, Exception, edm::ParameterSet::getParameter(), kernel2d, fftjeteflowsmoother_cfi::kernelEtaScale, fftjeteflowsmoother_cfi::kernelPhiScale, and M_PI.

Referenced by FFTJetPileupProcessor().

                                                                           {
   // Get the eta and phi scales for the kernel(s)
   double kernelEtaScale = ps.getParameter<double>("kernelEtaScale");
   const double kernelPhiScale = ps.getParameter<double>("kernelPhiScale");
   if (kernelEtaScale <= 0.0 || kernelPhiScale <= 0.0)
     throw cms::Exception("FFTJetBadConfig") << "invalid kernel scales" << std::endl;
 
   if (convolverMinBin || convolverMaxBin)
     if (convolverMinBin >= convolverMaxBin || convolverMaxBin > energyFlow->nEta())
       throw cms::Exception("FFTJetBadConfig") << "invalid convolver bin range" << std::endl;
 
   // FFT assumes that the grid extent in eta is 2*Pi. Adjust the
   // kernel scale in eta to compensate.
   kernelEtaScale *= (2.0 * M_PI / (energyFlow->etaMax() - energyFlow->etaMin()));
 
   // Build the FFT engine
   engine = std::unique_ptr<MyFFTEngine>(new MyFFTEngine(energyFlow->nEta(), energyFlow->nPhi()));
 
   // 2d kernel
   kernel2d = std::unique_ptr<fftjet::AbsFrequencyKernel>(
       new fftjet::DiscreteGauss2d(kernelEtaScale, kernelPhiScale, energyFlow->nEta(), energyFlow->nPhi()));
 
   // 2d convolver
   convolver = std::unique_ptr<fftjet::AbsConvolverBase<Real>>(new fftjet::FrequencyKernelConvolver<Real, Complex>(
       engine.get(), kernel2d.get(), convolverMinBin, convolverMaxBin));
 }

void FFTJetPileupProcessor::endJob ( void )

overrideprotectedvirtual

Reimplemented from edm::EDProducer.

Definition at line 325 of file FFTJetPileupProcessor.cc.

325 {}

void FFTJetPileupProcessor::loadFlatteningFactors ( const edm::EventSetup & iSetup )

private

Definition at line 327 of file FFTJetPileupProcessor.cc.

References DEFINE_FWK_MODULE, fftjetcms::FFTJetInterface::energyFlow, PVValHelper::eta, etaFlatteningFactors, f, flatteningTableCategory, flatteningTableName, flatteningTableRecord, h, LEDCalibrationChannels::ieta, StaticFFTJetRcdMapper< Mapper >::instance(), HLT_2018_cff::nEta, and mitigatedMETSequence_cff::U.

Referenced by produce().

                                                                              {
   edm::ESHandle<FFTJetLookupTableSequence> h;
   StaticFFTJetLookupTableSequenceLoader::instance().load(iSetup, flatteningTableRecord, h);
   std::shared_ptr<npstat::StorableMultivariateFunctor> f = (*h)[flatteningTableCategory][flatteningTableName];
 
   // Fill out the table of flattening factors as a function of eta
   const unsigned nEta = energyFlow->nEta();
   etaFlatteningFactors.clear();
   etaFlatteningFactors.reserve(nEta);
   for (unsigned ieta = 0; ieta < nEta; ++ieta) {
     const double eta = energyFlow->etaBinCenter(ieta);
     etaFlatteningFactors.push_back((*f)(&eta, 1U));
   }
 }

void FFTJetPileupProcessor::mixExtraGrid ( )

private

Definition at line 268 of file FFTJetPileupProcessor.cc.

References fftjetcms::add_Grid2d_data(), currentFileNum, fftjetcms::FFTJetInterface::energyFlow, Exception, externalGridFiles, externalGridMaxEnergy, g, gridStream, recoMuon::in, cmsBatch::nFiles, mergeAndRegister::ntries, readEcalDQMStatus::read, and mitigatedMETSequence_cff::U.

Referenced by produce().

                                          {
   const unsigned nFiles = externalGridFiles.size();
   if (currentFileNum > nFiles) {
     // This is the first time this function is called
     currentFileNum = 0;
     gridStream.open(externalGridFiles[currentFileNum].c_str(), std::ios_base::in | std::ios_base::binary);
     if (!gridStream.is_open())
       throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
                                                  " failed to open external grid file "
                                               << externalGridFiles[currentFileNum] << std::endl;
   }
 
   const fftjet::Grid2d<float>* g = nullptr;
   const unsigned maxFail = 100U;
   unsigned nEnergyRejected = 0;
 
   while (!g) {
     g = fftjet::Grid2d<float>::read(gridStream);
 
     // If we can't read the grid, we need to switch to another file
     for (unsigned ntries = 0; ntries < nFiles && g == nullptr; ++ntries) {
       gridStream.close();
       currentFileNum = (currentFileNum + 1U) % nFiles;
       gridStream.open(externalGridFiles[currentFileNum].c_str(), std::ios_base::in | std::ios_base::binary);
       if (!gridStream.is_open())
         throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
                                                    " failed to open external grid file "
                                                 << externalGridFiles[currentFileNum] << std::endl;
       g = fftjet::Grid2d<float>::read(gridStream);
     }
 
     if (g)
       if (g->sum() > externalGridMaxEnergy) {
         delete g;
         g = nullptr;
         if (++nEnergyRejected >= maxFail)
           throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
                                                      " too many grids in a row ("
                                                   << nEnergyRejected << ") failed the maximum energy cut" << std::endl;
       }
   }
 
   if (g) {
     add_Grid2d_data(energyFlow.get(), *g);
     delete g;
   } else {
     // Too bad, no useful file found
     throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
                                                " no valid grid records found"
                                             << std::endl;
   }
 }

FFTJetPileupProcessor& FFTJetPileupProcessor::operator= ( const FFTJetPileupProcessor & )

privatedelete

void FFTJetPileupProcessor::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

overrideprotected

Definition at line 203 of file FFTJetPileupProcessor.cc.

References convolvedFlow, convolver, convolverMaxBin, convolverMinBin, fftjetcms::FFTJetInterface::discretizeEnergyFlow(), fftjetcms::FFTJetInterface::energyFlow, etaFlatteningFactors, externalGridFiles, filterScales, g, loadFlatteningFactors(), loadFlatteningFactorsFromDB, fftjetcms::FFTJetInterface::loadInputCollection(), mixExtraGrid(), nPercentiles, fftjetcommon_cfi::nScales, fftjetcms::FFTJetInterface::outputLabel, percentileData, pileupEtaPhiArea, edm::Event::put(), data-class-funcs::q, and mitigatedMETSequence_cff::U.

                                                                                  {
   loadInputCollection(iEvent);
   discretizeEnergyFlow();
 
   // Determine the average Et density for this event.
   // Needs to be done here, before mixing in another grid.
   const fftjet::Grid2d<Real>& g(*energyFlow);
   const double densityBeforeMixing = g.sum() / pileupEtaPhiArea;
 
   // Mix an extra grid (if requested)
   double densityAfterMixing = -1.0;
   if (!externalGridFiles.empty()) {
     mixExtraGrid();
     densityAfterMixing = g.sum() / pileupEtaPhiArea;
   }
 
   // Various useful variables
   const unsigned nScales = filterScales->size();
   const double* scales = &(*filterScales)[0];
   Real* convData = const_cast<Real*>(convolvedFlow->data());
   Real* sortData = convData + convolverMinBin * convolvedFlow->nPhi();
   const unsigned dataLen = convolverMaxBin ? (convolverMaxBin - convolverMinBin) * convolvedFlow->nPhi()
                                            : convolvedFlow->nEta() * convolvedFlow->nPhi();
 
   // Load flattenning factors from DB (if requested)
   if (loadFlatteningFactorsFromDB)
     loadFlatteningFactors(iSetup);
 
   // Go over all scales and perform the convolutions
   convolver->setEventData(g.data(), g.nEta(), g.nPhi());
   for (unsigned iscale = 0; iscale < nScales; ++iscale) {
     // Perform the convolution
     convolver->convolveWithKernel(scales[iscale], convData, convolvedFlow->nEta(), convolvedFlow->nPhi());
 
     // Apply the flattening factors
     if (!etaFlatteningFactors.empty())
       convolvedFlow->scaleData(&etaFlatteningFactors[0], etaFlatteningFactors.size());
 
     // Sort the convolved data
     std::sort(sortData, sortData + dataLen);
 
     // Determine the percentile points
     for (unsigned iper = 0; iper < nPercentiles; ++iper) {
       // Map percentile 0 into point number 0,
       // 1 into point number dataLen-1
       const double q = (iper + 0.5) / nPercentiles;
       const double dindex = q * (dataLen - 1U);
       const unsigned ilow = static_cast<unsigned>(std::floor(dindex));
       const double percentile =
           fftjet::lin_interpolate_1d(ilow, ilow + 1U, sortData[ilow], sortData[ilow + 1U], dindex);
 
       // Store the calculated percentile
       percentileData[iscale * nPercentiles + iper] = percentile;
     }
   }
 
   // Convert percentile data into a more convenient storable object
   // and put it into the event record
   iEvent.put(std::make_unique<reco::DiscretizedEnergyFlow>(
                  &percentileData[0], "FFTJetPileupProcessor", -0.5, nScales - 0.5, 0.0, nScales, nPercentiles),
              outputLabel);
 
   iEvent.put(std::make_unique<std::pair<double, double>>(densityBeforeMixing, densityAfterMixing), outputLabel);
 }