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

Inheritance diagram for FFTJetPileupProcessor:

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

	~FFTJetPileupProcessor ()

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

ModuleDescription const &	moduleDescription () const

virtual	~EDProducer ()

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

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

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

virtual	~ProducerBase ()

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

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

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

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

std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

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

Protected Member Functions
void	beginJob () override

void	endJob () override

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

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

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 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)

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

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

	FFTJetPileupProcessor ()

	FFTJetPileupProcessor (const FFTJetPileupProcessor &)

void	loadFlatteningFactors (const edm::EventSetup &iSetup)

void	mixExtraGrid ()

FFTJetPileupProcessor &	operator= (const FFTJetPileupProcessor &)

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

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

unsigned	convolverMaxBin

unsigned	convolverMinBin

unsigned	currentFileNum

std::auto_ptr< MyFFTEngine >	engine

std::vector< double >	etaFlatteningFactors

std::vector< std::string >	externalGridFiles

double	externalGridMaxEnergy

std::auto_ptr < fftjet::EquidistantInLogSpace >	filterScales

std::string	flatteningTableCategory

std::string	flatteningTableName

std::string	flatteningTableRecord

std::ifstream	gridStream

std::auto_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
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Attributes inherited from fftjetcms::FFTJetInterface
const AnomalousTower	anomalous

std::vector< unsigned >	candidateIndex

const bool	doPVCorrection

std::auto_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 56 of file FFTJetPileupProcessor.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 128 of file FFTJetPileupProcessor.cc.

References buildKernelConvolver(), fftjetcms::FFTJetInterface::checkConfig(), convolvedFlow, fftjetcms::FFTJetInterface::energyFlow, etaFlatteningFactors, edm::hlt::Exception, fftjetcms::fftjet_Grid2d_parser(), filterScales, edm::ParameterSet::getParameter(), nPercentiles, fftjetcms::FFTJetInterface::outputLabel, and percentileData.

     : 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::auto_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::auto_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 ( )

Definition at line 223 of file FFTJetPileupProcessor.cc.

224 {

225 }

FFTJetPileupProcessor::FFTJetPileupProcessor ( )

private

FFTJetPileupProcessor::FFTJetPileupProcessor ( const FFTJetPileupProcessor & )

private

Member Function Documentation

void FFTJetPileupProcessor::beginJob ( void )

overrideprotectedvirtual

Reimplemented from edm::EDProducer.

Definition at line 376 of file FFTJetPileupProcessor.cc.

377 {

378 }

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

private

Definition at line 186 of file FFTJetPileupProcessor.cc.

References convolver, convolverMaxBin, convolverMinBin, fftjetcms::FFTJetInterface::energyFlow, engine, edm::hlt::Exception, edm::ParameterSet::getParameter(), kernel2d, 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::auto_ptr<MyFFTEngine>(
         new MyFFTEngine(energyFlow->nEta(), energyFlow->nPhi()));
 
     // 2d kernel
     kernel2d = std::auto_ptr<fftjet::AbsFrequencyKernel>(
         new fftjet::DiscreteGauss2d(
             kernelEtaScale, kernelPhiScale,
             energyFlow->nEta(), energyFlow->nPhi()));
 
     // 2d convolver
     convolver = std::auto_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 382 of file FFTJetPileupProcessor.cc.

383 {

384 }

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

private

Definition at line 387 of file FFTJetPileupProcessor.cc.

References fftjetcms::FFTJetInterface::energyFlow, eta(), etaFlatteningFactors, f, flatteningTableCategory, flatteningTableName, flatteningTableRecord, h, and StaticFFTJetRcdMapper< Mapper >::instance().

Referenced by produce().

 {
     edm::ESHandle<FFTJetLookupTableSequence> h;
     StaticFFTJetLookupTableSequenceLoader::instance().load(
         iSetup, flatteningTableRecord, h);
     boost::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 308 of file FFTJetPileupProcessor.cc.

References fftjetcms::add_Grid2d_data(), currentFileNum, fftjetcms::FFTJetInterface::energyFlow, edm::hlt::Exception, externalGridFiles, externalGridMaxEnergy, g, gridStream, recoMuon::in, and SiPixelLorentzAngle_cfi::read.

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 = 0;
     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 == 0; ++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 = 0;
                 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 & )

private

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

overrideprotectedvirtual

Implements edm::EDProducer.

Definition at line 229 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, fftjetcms::FFTJetInterface::outputLabel, percentileData, pileupEtaPhiArea, edm::Event::put(), lumiQueryAPI::q, and python.multivaluedict::sort().

 {
     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
     std::auto_ptr<reco::DiscretizedEnergyFlow> pTable(
         new reco::DiscretizedEnergyFlow(
             &percentileData[0], "FFTJetPileupProcessor",
             -0.5, nScales-0.5, 0.0, nScales, nPercentiles));
     iEvent.put(pTable, outputLabel);
 
     std::auto_ptr<std::pair<double,double> > etSum(
         new std::pair<double,double>(densityBeforeMixing, densityAfterMixing));
     iEvent.put(etSum, outputLabel);
 }