#include <RecoJets/FFTJetProducer/plugins/FFTJetPatRecoProducer.cc>

Inheritance diagram for FFTJetPatRecoProducer:

Public Member Functions
	FFTJetPatRecoProducer (const edm::ParameterSet &)
	~FFTJetPatRecoProducer ()
Protected Types
typedef fftjet::ProximityClusteringTree < fftjet::Peak, long >	ClusteringTree
typedef fftjet::ClusteringSequencer < Real >	Sequencer
typedef fftjet::SparseClusteringTree < fftjet::Peak, long >	SparseTree
typedef fftjet::ClusteringTreeSparsifier < fftjet::Peak, long >	Sparsifier
Protected Member Functions
void	beginJob ()
template<class Real >
void	buildDenseProduct (edm::Event &) const
void	buildKernelConvolver (const edm::ParameterSet &)
fftjet::PeakFinder	buildPeakFinder (const edm::ParameterSet &)
template<class Real >
void	buildSparseProduct (edm::Event &) const
void	endJob ()
void	produce (edm::Event &, const edm::EventSetup &)
Protected Attributes
std::auto_ptr< MyFFTEngine >	anotherEngine
ClusteringTree *	clusteringTree
const double	completeEventDataCutoff
std::auto_ptr < fftjet::AbsConvolverBase < Real > >	convolver
std::auto_ptr < fftjet::AbsDistanceCalculator < fftjet::Peak > >	distanceCalc
std::auto_ptr< MyFFTEngine >	engine
std::auto_ptr < fftjet::AbsFrequencyKernel1d >	etaKernel
std::auto_ptr < fftjet::AbsFrequencyKernel >	kernel2d
const bool	makeClusteringTree
std::auto_ptr < fftjet::Functor1< bool, fftjet::Peak > >	peakSelector
std::auto_ptr < fftjet::AbsFrequencyKernel1d >	phiKernel
std::auto_ptr< Sequencer >	sequencer
SparseTree	sparseTree
std::auto_ptr< Sparsifier >	sparsifier
const bool	sparsify
const bool	storeDiscretizationGrid
const bool	verifyDataConversion
Private Member Functions
	FFTJetPatRecoProducer ()
	FFTJetPatRecoProducer (const FFTJetPatRecoProducer &)
FFTJetPatRecoProducer &	operator= (const FFTJetPatRecoProducer &)
Private Attributes
std::ofstream	externalGridStream
fftjet::Grid2d< float > *	extGrid
bool	storeGridsExternally

Detailed Description

Description: Runs FFTJet pattern recognition stage and saves the results

Implementation: [Notes on implementation]

Definition at line 67 of file FFTJetPatRecoProducer.cc.

Member Typedef Documentation

typedef fftjet::ProximityClusteringTree<fftjet::Peak,long> FFTJetPatRecoProducer::ClusteringTree [protected]

Definition at line 75 of file FFTJetPatRecoProducer.cc.

typedef fftjet::ClusteringSequencer<Real> FFTJetPatRecoProducer::Sequencer [protected]

Definition at line 77 of file FFTJetPatRecoProducer.cc.

typedef fftjet::SparseClusteringTree<fftjet::Peak,long> FFTJetPatRecoProducer::SparseTree [protected]

Definition at line 76 of file FFTJetPatRecoProducer.cc.

typedef fftjet::ClusteringTreeSparsifier<fftjet::Peak,long> FFTJetPatRecoProducer::Sparsifier [protected]

Definition at line 78 of file FFTJetPatRecoProducer.cc.

Constructor & Destructor Documentation

FFTJetPatRecoProducer::FFTJetPatRecoProducer ( const edm::ParameterSet & ps ) [explicit]

Definition at line 156 of file FFTJetPatRecoProducer.cc.

References buildKernelConvolver(), buildPeakFinder(), fftjetcms::FFTJetInterface::checkConfig(), clusteringTree, convolver, fwrapper::cs, distanceCalc, fftjetcms::FFTJetInterface::energyFlow, Exception, externalGridStream, fftjetcms::fftjet_ClusteringTreeSparsifier_parser(), fftjetcms::fftjet_DistanceCalculator_parser(), fftjetcms::fftjet_Grid2d_parser(), fftjetcms::fftjet_PeakSelector_parser(), fftjetcms::fftjet_ScaleSet_parser(), fftjetcms::FFTJetInterface::getEventScale(), edm::ParameterSet::getParameter(), i, makeClusteringTree, dbtoconf::out, fftjetcms::FFTJetInterface::outputLabel, peakSelector, sequencer, sparsifier, storeDiscretizationGrid, storeGridsExternally, and fftjetcms::FFTJetInterface::storeInSinglePrecision().

    : FFTJetInterface(ps),
      clusteringTree(0),
      completeEventDataCutoff(ps.getParameter<double>("completeEventDataCutoff")),
      makeClusteringTree(ps.getParameter<bool>("makeClusteringTree")),
      verifyDataConversion(ps.getUntrackedParameter<bool>("verifyDataConversion",false)),
      storeDiscretizationGrid(ps.getParameter<bool>("storeDiscretizationGrid")),
      sparsify(ps.getParameter<bool>("sparsify")),
      extGrid(0)
{
    // register your products
    if (makeClusteringTree)
    {
        if (storeInSinglePrecision())
            produces<reco::PattRecoTree<float,reco::PattRecoPeak<float> > >(outputLabel);
        else
            produces<reco::PattRecoTree<double,reco::PattRecoPeak<double> > >(outputLabel);
    }
    if (storeDiscretizationGrid)
        produces<reco::DiscretizedEnergyFlow>(outputLabel);

    // Check if we want to write the grids into an external file
    const std::string externalGridFile(ps.getParameter<std::string>("externalGridFile"));
    storeGridsExternally = externalGridFile.size() > 0;
    if (storeGridsExternally)
    {
        externalGridStream.open(externalGridFile.c_str(), std::ios_base::out | 
                                                          std::ios_base::binary);
        if (!externalGridStream.is_open())
            throw cms::Exception("FFTJetBadConfig")
                << "FFTJetPatRecoProducer failed to open file "
                << externalGridFile << std::endl;
    }

    if (!makeClusteringTree && !storeDiscretizationGrid && !storeGridsExternally)
    {
        throw cms::Exception("FFTJetBadConfig")
            << "FFTJetPatRecoProducer is not configured to produce anything"
            << std::endl;
    }

    // now do whatever other initialization is needed

    // Build the discretization grid
    energyFlow = fftjet_Grid2d_parser(
        ps.getParameter<edm::ParameterSet>("GridConfiguration"));
    checkConfig(energyFlow, "invalid discretization grid");

    // Build the FFT engine(s), pattern recognition kernel(s),
    // and the kernel convolver
    buildKernelConvolver(ps);

    // Build the peak selector
    peakSelector = fftjet_PeakSelector_parser(
        ps.getParameter<edm::ParameterSet>("PeakSelectorConfiguration"));
    checkConfig(peakSelector, "invalid peak selector");

    // Build the initial set of pattern recognition scales
    std::auto_ptr<std::vector<double> > iniScales = fftjet_ScaleSet_parser(
        ps.getParameter<edm::ParameterSet>("InitialScales"));
    checkConfig(iniScales, "invalid set of scales");

    // Do we want to use the adaptive clustering tree algorithm?
    const unsigned maxAdaptiveScales = 
        ps.getParameter<unsigned>("maxAdaptiveScales");
    const double minAdaptiveRatioLog = 
        ps.getParameter<double>("minAdaptiveRatioLog");
    if (minAdaptiveRatioLog <= 0.0)
        throw cms::Exception("FFTJetBadConfig")
            << "bad adaptive ratio logarithm limit" << std::endl;

    // Make sure that all standard scales are larger than the
    // complete event scale
    if (getEventScale() > 0.0)
    {
      const double cs = getEventScale();
      const unsigned nscales = iniScales->size();
      for (unsigned i=0; i<nscales; ++i)
        if (cs >= (*iniScales)[i])
          throw cms::Exception("FFTJetBadConfig")
            << "incompatible scale for complete event" << std::endl;
    }

    // At this point we are ready to build the clustering sequencer
    sequencer = std::auto_ptr<Sequencer>(new Sequencer(
        convolver.get(), peakSelector.get(), buildPeakFinder(ps),
        *iniScales, maxAdaptiveScales, minAdaptiveRatioLog));

    // Build the clustering tree sparsifier
    const edm::ParameterSet& SparsifierConfiguration(
        ps.getParameter<edm::ParameterSet>("SparsifierConfiguration"));
    sparsifier = fftjet_ClusteringTreeSparsifier_parser(
        SparsifierConfiguration);
    checkConfig(sparsifier, "invalid sparsifier parameters");

    // Build the distance calculator for the clustering tree
    const edm::ParameterSet& TreeDistanceCalculator(
        ps.getParameter<edm::ParameterSet>("TreeDistanceCalculator"));
    distanceCalc = fftjet_DistanceCalculator_parser(TreeDistanceCalculator);
    checkConfig(distanceCalc, "invalid tree distance calculator");

    // Build the clustering tree itself
    clusteringTree = new ClusteringTree(distanceCalc.get());
}

FFTJetPatRecoProducer::~FFTJetPatRecoProducer ( )

Definition at line 360 of file FFTJetPatRecoProducer.cc.

References clusteringTree, and extGrid.

{
    // do anything here that needs to be done at desctruction time
    // (e.g. close files, deallocate resources etc.)
    delete clusteringTree;
    delete extGrid;
}

FFTJetPatRecoProducer::FFTJetPatRecoProducer ( ) [private]

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

Member Function Documentation

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

Reimplemented from edm::EDProducer.

Definition at line 501 of file FFTJetPatRecoProducer.cc.

{
}

template<class Real >

void FFTJetPatRecoProducer::buildDenseProduct ( edm::Event & ev ) const [protected]

Definition at line 400 of file FFTJetPatRecoProducer.cc.

References CastorDataFrameFilter_impl::check(), clusteringTree, fftjetcms::densePeakTreeFromStorable(), fftjetcms::densePeakTreeToStorable(), distanceCalc, Exception, fftjetcms::FFTJetInterface::getEventScale(), fftjetcms::FFTJetInterface::outputLabel, edm::Event::put(), sequencer, fftjetcms::FFTJetInterface::storeInSinglePrecision(), diffTreeTool::tree, and verifyDataConversion.

{
    typedef reco::PattRecoTree<Real,reco::PattRecoPeak<Real> > StoredTree;

    std::auto_ptr<StoredTree> tree(new StoredTree());

    densePeakTreeToStorable(*clusteringTree,
                            sequencer->maxAdaptiveScales(),
                            tree.get());

    // Check that we can restore the tree
    if (verifyDataConversion && !storeInSinglePrecision())
    {
        ClusteringTree check(distanceCalc.get());
        const std::vector<double>& scalesUsed(sequencer->getInitialScales());
        densePeakTreeFromStorable(*tree, &scalesUsed, getEventScale(), &check);
        if (*clusteringTree != check)
            throw cms::Exception("FFTJetInterface")
                << "Data conversion failed for dense clustering tree"
                << std::endl;
    }

    ev.put(tree, outputLabel);
}

void FFTJetPatRecoProducer::buildKernelConvolver ( const edm::ParameterSet & ps ) [protected]

Definition at line 262 of file FFTJetPatRecoProducer.cc.

References anotherEngine, convolver, fftjetcms::FFTJetInterface::energyFlow, engine, etaKernel, Exception, edm::ParameterSet::getParameter(), kernel2d, M_PI, and phiKernel.

Referenced by FFTJetPatRecoProducer().

{
    // Check the parameter named "etaDependentScaleFactors". If the vector
    // of scales is empty we will use 2d kernel, otherwise use 1d kernels
    const std::vector<double> etaDependentScaleFactors(
        ps.getParameter<std::vector<double> >("etaDependentScaleFactors"));

    // Make sure that the number of scale factors provided is correct
    const bool use2dKernel = etaDependentScaleFactors.empty();
    if (!use2dKernel)
        if (etaDependentScaleFactors.size() != energyFlow->nEta())
            throw cms::Exception("FFTJetBadConfig")
                << "invalid number of eta-dependent scale factors"
                << std::endl;

    // 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 scale" << 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()));

    // Are we going to try to fix the efficiency near detector edges?
    const bool fixEfficiency = ps.getParameter<bool>("fixEfficiency");

    // Minimum and maximum eta bin for the convolver
    unsigned convolverMinBin = 0, convolverMaxBin = 0;
    if (fixEfficiency || !use2dKernel)
    {
        convolverMinBin = ps.getParameter<unsigned>("convolverMinBin");
        convolverMaxBin = ps.getParameter<unsigned>("convolverMaxBin");
    }

    if (use2dKernel)
    {
        // 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));
    }
    else
    {
        // Need separate FFT engines for eta and phi
        engine = std::auto_ptr<MyFFTEngine>(
            new MyFFTEngine(1, energyFlow->nEta()));
        anotherEngine = std::auto_ptr<MyFFTEngine>(
            new MyFFTEngine(1, energyFlow->nPhi()));

        // 1d kernels
        etaKernel = std::auto_ptr<fftjet::AbsFrequencyKernel1d>(
            new fftjet::DiscreteGauss1d(kernelEtaScale, energyFlow->nEta()));

        phiKernel = std::auto_ptr<fftjet::AbsFrequencyKernel1d>(
            new fftjet::DiscreteGauss1d(kernelPhiScale, energyFlow->nPhi()));

        // Sequential convolver
        convolver = std::auto_ptr<fftjet::AbsConvolverBase<Real> >(
            new fftjet::FrequencySequentialConvolver<Real,Complex>(
                engine.get(), anotherEngine.get(),
                etaKernel.get(), phiKernel.get(),
                etaDependentScaleFactors, convolverMinBin,
                convolverMaxBin, fixEfficiency));
    }
}

fftjet::PeakFinder FFTJetPatRecoProducer::buildPeakFinder ( const edm::ParameterSet & ps ) [protected]

Definition at line 343 of file FFTJetPatRecoProducer.cc.

References fftjetcms::FFTJetInterface::energyFlow, Exception, and edm::ParameterSet::getParameter().

Referenced by FFTJetPatRecoProducer().

{
    const double peakFinderMaxEta = ps.getParameter<double>("peakFinderMaxEta");
    if (peakFinderMaxEta <= 0.0)
        throw cms::Exception("FFTJetBadConfig")
            << "invalid peak finder eta cut" << std::endl;
    const double maxMagnitude = ps.getParameter<double>("peakFinderMaxMagnitude");
    int minBin = energyFlow->getEtaBin(-peakFinderMaxEta);
    if (minBin < 0)
        minBin = 0;
    int maxBin = energyFlow->getEtaBin(peakFinderMaxEta) + 1;
    if (maxBin < 0)
        maxBin = 0;
    return fftjet::PeakFinder(maxMagnitude, true, minBin, maxBin);
}

template<class Real >

void FFTJetPatRecoProducer::buildSparseProduct ( edm::Event & ev ) const [protected]

Definition at line 373 of file FFTJetPatRecoProducer.cc.

References CastorDataFrameFilter_impl::check(), Exception, fftjetcms::FFTJetInterface::getEventScale(), fftjetcms::FFTJetInterface::outputLabel, edm::Event::put(), sequencer, fftjetcms::sparsePeakTreeFromStorable(), fftjetcms::sparsePeakTreeToStorable(), sparseTree, fftjetcms::FFTJetInterface::storeInSinglePrecision(), diffTreeTool::tree, and verifyDataConversion.

{
    typedef reco::PattRecoTree<Real,reco::PattRecoPeak<Real> > StoredTree;

    std::auto_ptr<StoredTree> tree(new StoredTree());

    sparsePeakTreeToStorable(sparseTree,
                             sequencer->maxAdaptiveScales(),
                             tree.get());

    // Check that we can restore the tree
    if (verifyDataConversion && !storeInSinglePrecision())
    {
        SparseTree check;
        const std::vector<double>& scalesUsed(sequencer->getInitialScales());
        sparsePeakTreeFromStorable(*tree, &scalesUsed, getEventScale(), &check);
        if (sparseTree != check)
            throw cms::Exception("FFTJetInterface")
                << "Data conversion failed for sparse clustering tree"
                << std::endl;
    }

    ev.put(tree, outputLabel);
}

void FFTJetPatRecoProducer::endJob ( void ) [protected, virtual]

Reimplemented from edm::EDProducer.

Definition at line 507 of file FFTJetPatRecoProducer.cc.

References externalGridStream, and storeGridsExternally.

{
    if (storeGridsExternally)
        externalGridStream.close();
}

FFTJetPatRecoProducer& FFTJetPatRecoProducer::operator= ( const FFTJetPatRecoProducer & ) [private]

void FFTJetPatRecoProducer::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)		`[protected, virtual]`

Implements edm::EDProducer.

Definition at line 427 of file FFTJetPatRecoProducer.cc.

References CastorDataFrameFilter_impl::check(), clusteringTree, completeEventDataCutoff, fftjetcms::convert_Grid2d_to_float(), fftjetcms::copy_Grid2d_data(), fftjetcms::FFTJetInterface::discretizeEnergyFlow(), fftjetcms::FFTJetInterface::energyFlow, Exception, externalGridStream, extGrid, g, fftjetcms::FFTJetInterface::getEventScale(), iEvent, fftjetcms::FFTJetInterface::loadInputCollection(), makeClusteringTree, fftjetcms::FFTJetInterface::outputLabel, edm::Event::put(), sequencer, sparseTree, sparsifier, sparsify, storeDiscretizationGrid, storeGridsExternally, fftjetcms::FFTJetInterface::storeInSinglePrecision(), and verifyDataConversion.

{
    loadInputCollection(iEvent);
    discretizeEnergyFlow();

    if (makeClusteringTree)
    {
        sequencer->run(*energyFlow, clusteringTree);
        if (getEventScale() > 0.0)
            sequencer->insertCompleteEvent(getEventScale(), *energyFlow,
                                           clusteringTree, completeEventDataCutoff);

        if (sparsify)
        {
            sparsifier->sparsify(*clusteringTree, &sparseTree);

            // Do not call the "sortNodes" method of the sparse tree here.
            // Currently, the nodes are sorted by daughter number.
            // This is the way we want it in storage because the stored
            // tree does not include daughter ordering info explicitly.

            if (storeInSinglePrecision())
                buildSparseProduct<float>(iEvent);
            else
                buildSparseProduct<double>(iEvent);
        }
        else
        {
            if (storeInSinglePrecision())
                buildDenseProduct<float>(iEvent);
            else
                buildDenseProduct<double>(iEvent);
        }
    }

    if (storeDiscretizationGrid)
    {
        const fftjet::Grid2d<Real>& g(*energyFlow);

        std::auto_ptr<reco::DiscretizedEnergyFlow> flow(
            new reco::DiscretizedEnergyFlow(
                g.data(), g.title(), g.etaMin(), g.etaMax(),
                g.phiBin0Edge(), g.nEta(), g.nPhi()));

        if (verifyDataConversion)
        {
            fftjet::Grid2d<Real> check(
                flow->nEtaBins(), flow->etaMin(), flow->etaMax(),
                flow->nPhiBins(), flow->phiBin0Edge(), flow->title());
            check.blockSet(flow->data(), flow->nEtaBins(), flow->nPhiBins());
            assert(g == check);
        }

        iEvent.put(flow, outputLabel);
    }

    if (storeGridsExternally)
    {
        if (extGrid)
            copy_Grid2d_data(extGrid, *energyFlow);
        else
            extGrid = convert_Grid2d_to_float(*energyFlow);
        if (!extGrid->write(externalGridStream))
        {
            throw cms::Exception("FFTJetPatRecoProducer::produce")
                << "Failed to write grid data into an external file"
                << std::endl;
        }
    }
}