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

Inheritance diagram for FFTJetPileupProcessor:

Public Member Functions
	FFTJetPileupProcessor (const edm::ParameterSet &)
	~FFTJetPileupProcessor ()
Protected Member Functions
void	beginJob ()
void	endJob ()
void	produce (edm::Event &, const edm::EventSetup &)
Private Member Functions
void	buildKernelConvolver (const edm::ParameterSet &)
	FFTJetPileupProcessor ()
	FFTJetPileupProcessor (const FFTJetPileupProcessor &)
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
std::auto_ptr < fftjet::EquidistantInLogSpace >	filterScales
std::ifstream	gridStream
std::auto_ptr < fftjet::AbsFrequencyKernel >	kernel2d
unsigned	nPercentiles
double	pileupEtaPhiArea

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 110 of file FFTJetPileupProcessor.cc.

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

    : 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")),
      currentFileNum(externalGridFiles.size() + 1U)
{
    // 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 iscretization 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));

    produces<TH2D>(outputLabel);
    produces<std::pair<double,double> >(outputLabel);
}

FFTJetPileupProcessor::~FFTJetPileupProcessor ( )

Definition at line 198 of file FFTJetPileupProcessor.cc.

{
}

FFTJetPileupProcessor::FFTJetPileupProcessor ( ) [private]

FFTJetPileupProcessor::FFTJetPileupProcessor ( const FFTJetPileupProcessor & ) [private]

Member Function Documentation

void FFTJetPileupProcessor::beginJob ( void ) [protected, virtual]

Reimplemented from edm::EDProducer.

Definition at line 332 of file FFTJetPileupProcessor.cc.

{
}

void FFTJetPileupProcessor::buildKernelConvolver ( const edm::ParameterSet & ps ) [private]

Definition at line 161 of file FFTJetPileupProcessor.cc.

References convolver, convolverMaxBin, convolverMinBin, fftjetcms::FFTJetInterface::energyFlow, engine, 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 ) [protected, virtual]

Reimplemented from edm::EDProducer.

Definition at line 338 of file FFTJetPileupProcessor.cc.

{
}

void FFTJetPileupProcessor::mixExtraGrid ( ) [private]

Definition at line 283 of file FFTJetPileupProcessor.cc.

References fftjetcms::add_Grid2d_data(), currentFileNum, fftjetcms::FFTJetInterface::energyFlow, Exception, externalGridFiles, 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 = 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)
    {
        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
	)		`[protected, virtual]`

Implements edm::EDProducer.

Definition at line 204 of file FFTJetPileupProcessor.cc.

References convolvedFlow, convolver, convolverMaxBin, convolverMinBin, fftjetcms::FFTJetInterface::discretizeEnergyFlow(), fftjetcms::FFTJetInterface::energyFlow, etaFlatteningFactors, externalGridFiles, filterScales, g, fftjetcms::FFTJetInterface::loadInputCollection(), mixExtraGrid(), nPercentiles, fftjetcms::FFTJetInterface::outputLabel, 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();

    // We will fill the following histo
    std::auto_ptr<TH2D> pTable(new TH2D("FFTJetPileupProcessor",
                                        "FFTJetPileupProcessor",
                                        nScales, -0.5, nScales-0.5,
                                        nPercentiles, 0.0, 1.0));
    pTable->SetDirectory(0);
    pTable->GetXaxis()->SetTitle("Filter Number");
    pTable->GetYaxis()->SetTitle("Et CDF");
    pTable->GetZaxis()->SetTitle("Et Density");

    // 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
            pTable->SetBinContent(iscale+1U, iper+1U, percentile);
        }
    }

    iEvent.put(pTable, outputLabel);

    std::auto_ptr<std::pair<double,double> > etSum(
        new std::pair<double,double>(densityBeforeMixing, densityAfterMixing));
    iEvent.put(etSum, outputLabel);
}