#include <RecoJets/FFTJetProducers/plugins/FFTJetProducer.h>

Inheritance diagram for FFTJetProducer:

Public Types
typedef fftjet::RecombinedJet < fftjetcms::VectorLike >	RecoFFTJet
enum	Resolution { FIXED = 0, MAXIMALLY_STABLE, GLOBALLY_ADAPTIVE, LOCALLY_ADAPTIVE }
typedef fftjet::SparseClusteringTree < fftjet::Peak, long >	SparseTree
Public Member Functions
	FFTJetProducer (const edm::ParameterSet &)
virtual	~FFTJetProducer ()
Static Public Member Functions
static Resolution	parse_resolution (const std::string &name)
Protected Member Functions
virtual void	assignMembershipFunctions (std::vector< fftjet::Peak > *preclusters)
virtual void	beginJob ()
virtual void	endJob ()
virtual std::auto_ptr < fftjetcms::AbsBgFunctor >	parse_bgMembershipFunction (const edm::ParameterSet &)
virtual std::auto_ptr < fftjet::Functor2< double, RecoFFTJet, RecoFFTJet > >	parse_jetDistanceCalc (const edm::ParameterSet &)
virtual std::auto_ptr < fftjet::ScaleSpaceKernel >	parse_jetMembershipFunction (const edm::ParameterSet &)
virtual std::auto_ptr < fftjet::Functor1< double, RecoFFTJet > >	parse_memberFactorCalcJet (const edm::ParameterSet &)
virtual std::auto_ptr < fftjet::Functor1< double, fftjet::Peak > >	parse_memberFactorCalcPeak (const edm::ParameterSet &)
virtual std::auto_ptr < fftjet::Functor1< bool, fftjet::Peak > >	parse_peakSelector (const edm::ParameterSet &)
virtual std::auto_ptr < fftjet::Functor1< double, RecoFFTJet > >	parse_recoScaleCalcJet (const edm::ParameterSet &)
virtual std::auto_ptr < fftjet::Functor1< double, fftjet::Peak > >	parse_recoScaleCalcPeak (const edm::ParameterSet &)
virtual std::auto_ptr < fftjet::Functor1< double, RecoFFTJet > >	parse_recoScaleRatioCalcJet (const edm::ParameterSet &)
virtual std::auto_ptr < fftjet::Functor1< double, fftjet::Peak > >	parse_recoScaleRatioCalcPeak (const edm::ParameterSet &)
virtual void	produce (edm::Event &, const edm::EventSetup &)
virtual void	selectPreclusters (const SparseTree &tree, const fftjet::Functor1< bool, fftjet::Peak > &peakSelector, std::vector< fftjet::Peak > *preclusters)
void	selectTreeNodes (const SparseTree &tree, const fftjet::Functor1< bool, fftjet::Peak > &peakSelect, std::vector< SparseTree::NodeId > *nodes)
Private Types
typedef fftjet::AbsRecombinationAlg < fftjetcms::Real, fftjetcms::VectorLike, fftjetcms::BgData >	GridAlg
typedef fftjet::AbsVectorRecombinationAlg < fftjetcms::VectorLike, fftjetcms::BgData >	RecoAlg
Private Member Functions
void	buildGridAlg ()
bool	checkConvergence (const std::vector< RecoFFTJet > &previousIterResult, std::vector< RecoFFTJet > &thisIterResult)
void	determineGriddedConstituents ()
void	determineVectorConstituents ()
	FFTJetProducer ()
	FFTJetProducer (const FFTJetProducer &)
unsigned	iterateJetReconstruction ()
void	loadEnergyFlow (const edm::Event &iEvent)
template<class Real >
void	loadSparseTreeData (const edm::Event &)
template<typename Jet >
void	makeProduces (const std::string &alias, const std::string &tag)
FFTJetProducer &	operator= (const FFTJetProducer &)
void	prepareRecombinationScales ()
void	saveResults (edm::Event &iEvent, const edm::EventSetup &)
template<typename Jet >
void	writeJets (edm::Event &iEvent, const edm::EventSetup &)
Private Attributes
const bool	assignConstituents
std::auto_ptr < fftjetcms::AbsBgFunctor >	bgMembershipFunction
std::vector< std::vector < reco::CandidatePtr > >	constituents
const double	convergenceDistance
const double	fixedScale
std::auto_ptr< GridAlg >	gridAlg
const double	gridScanMaxEta
std::auto_ptr< std::vector < double > >	iniScales
const bool	isCrisp
unsigned	iterationsPerformed
std::vector< RecoFFTJet >	iterJets
std::vector< fftjet::Peak >	iterPreclusters
std::auto_ptr < fftjet::Functor2< double, RecoFFTJet, RecoFFTJet > >	jetDistanceCalc
std::auto_ptr < fftjet::ScaleSpaceKernel >	jetMembershipFunction
const unsigned	maxIterations
unsigned	maxLevel
const double	maxStableScale
std::auto_ptr < fftjet::Functor1< double, RecoFFTJet > >	memberFactorCalcJet
std::auto_ptr < fftjet::Functor1< double, fftjet::Peak > >	memberFactorCalcPeak
unsigned	minLevel
const double	minStableScale
const edm::ParameterSet	myConfiguration
const unsigned	nClustersRequested
const unsigned	nJetsRequiredToConverge
std::vector< SparseTree::NodeId >	nodes
const double	noiseLevel
std::vector< unsigned >	occupancy
std::auto_ptr < fftjet::Functor1< bool, fftjet::Peak > >	peakSelector
std::vector< fftjet::Peak >	preclusters
std::auto_ptr< RecoAlg >	recoAlg
std::vector< RecoFFTJet >	recoJets
const std::string	recombinationAlgorithm
const double	recombinationDataCutoff
std::auto_ptr < fftjet::Functor1< double, RecoFFTJet > >	recoScaleCalcJet
std::auto_ptr < fftjet::Functor1< double, fftjet::Peak > >	recoScaleCalcPeak
std::auto_ptr < fftjet::Functor1< double, RecoFFTJet > >	recoScaleRatioCalcJet
std::auto_ptr < fftjet::Functor1< double, fftjet::Peak > >	recoScaleRatioCalcPeak
Resolution	resolution
const bool	resumConstituents
const bool	reuseExistingGrid
SparseTree	sparseTree
const double	stabilityAlpha
std::vector< double >	thresholds
const edm::InputTag	treeLabel
fftjetcms::VectorLike	unclustered
const double	unlikelyBgWeight
double	unused
unsigned	usedLevel
const bool	useGriddedAlgorithm

Detailed Description

Description: makes jets using FFTJet clustering tree

Implementation: [Notes on implementation]

Description: makes jets using FFTJet clustering tree

Implementation: If you want to change the jet algorithm functionality (for example, by providing your own jet membership function), derive from this class and override the appropriate parser method (for example, parse_jetMembershipFunction). At the end of your own parser, don't forget to call the parser of the base class in order to get the default behavior when your special configuration is not provided (this is known as the "chain-of-responsibility" design pattern). If you also need to override "beginJob" and/or "produce" methods, the first thing to do in your method is to call the corresponding method of this base.

Definition at line 57 of file FFTJetProducer.h.

Member Typedef Documentation

typedef fftjet::AbsRecombinationAlg< fftjetcms::Real,fftjetcms::VectorLike,fftjetcms::BgData> FFTJetProducer::GridAlg [private]

Definition at line 150 of file FFTJetProducer.h.

typedef fftjet::AbsVectorRecombinationAlg< fftjetcms::VectorLike,fftjetcms::BgData> FFTJetProducer::RecoAlg [private]

Definition at line 148 of file FFTJetProducer.h.

typedef fftjet::RecombinedJet<fftjetcms::VectorLike> FFTJetProducer::RecoFFTJet

Definition at line 61 of file FFTJetProducer.h.

typedef fftjet::SparseClusteringTree<fftjet::Peak,long> FFTJetProducer::SparseTree

Definition at line 62 of file FFTJetProducer.h.

Member Enumeration Documentation

enum FFTJetProducer::Resolution

Enumerator:

FIXED
MAXIMALLY_STABLE
GLOBALLY_ADAPTIVE
LOCALLY_ADAPTIVE

Definition at line 64 of file FFTJetProducer.h.

    {
        FIXED = 0,
        MAXIMALLY_STABLE,
        GLOBALLY_ADAPTIVE,
        LOCALLY_ADAPTIVE
    };

Constructor & Destructor Documentation

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

Definition at line 123 of file FFTJetProducer.cc.

References assignConstituents, fftjetcms::FFTJetInterface::checkConfig(), Exception, fftjetcms::fftjet_ScaleSet_parser(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), iniScales, jet_type_switch, makeProduces(), fftjetcms::FFTJetInterface::outputLabel, resumConstituents, and python::multivaluedict::sort().

    : FFTJetInterface(ps),
      myConfiguration(ps),
      init_param(edm::InputTag, treeLabel),
      init_param(bool, useGriddedAlgorithm),
      init_param(bool, reuseExistingGrid),
      init_param(unsigned, maxIterations),
      init_param(unsigned, nJetsRequiredToConverge),
      init_param(double, convergenceDistance),
      init_param(bool, assignConstituents),
      init_param(bool, resumConstituents),
      init_param(double, fixedScale),
      init_param(double, minStableScale),
      init_param(double, maxStableScale),
      init_param(double, stabilityAlpha),
      init_param(double, noiseLevel),
      init_param(unsigned, nClustersRequested),
      init_param(double, gridScanMaxEta),
      init_param(std::string, recombinationAlgorithm),
      init_param(bool, isCrisp),
      init_param(double, unlikelyBgWeight),
      init_param(double, recombinationDataCutoff),
      resolution(parse_resolution(ps.getParameter<std::string>("resolution"))),

      minLevel(0),
      maxLevel(0),
      usedLevel(0),
      unused(0.0),
      iterationsPerformed(1U),
      constituents(200)
{
    // Check that the settings make sense
    if (resumConstituents && !assignConstituents)
        throw cms::Exception("FFTJetBadConfig") 
            << "Can't resum constituents if they are not assigned"
            << std::endl;

    produces<reco::FFTJetProducerSummary>(outputLabel);
    const std::string alias(ps.getUntrackedParameter<std::string>(
                                "alias", outputLabel));
    jet_type_switch(makeProduces, alias, outputLabel);

    // Build the set of pattern recognition scales.
    // This is needed in order to read the clustering tree
    // from the event record.
    iniScales = fftjet_ScaleSet_parser(
        ps.getParameter<edm::ParameterSet>("InitialScales"));
    checkConfig(iniScales, "invalid set of scales");
    std::sort(iniScales->begin(), iniScales->end(), std::greater<double>());

    // Most of the configuration has to be performed inside
    // the "beginJob" method. This is because chaining of the
    // parsers between this base class and the derived classes
    // can not work from the constructor of the base class.
}

FFTJetProducer::~FFTJetProducer ( ) [virtual]

Definition at line 180 of file FFTJetProducer.cc.

{
}

FFTJetProducer::FFTJetProducer ( ) [private]

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

Member Function Documentation

void FFTJetProducer::assignMembershipFunctions ( std::vector< fftjet::Peak > * preclusters ) [protected, virtual]

Definition at line 804 of file FFTJetProducer.cc.

Referenced by produce().

{
}

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

Reimplemented from edm::EDProducer.

Definition at line 810 of file FFTJetProducer.cc.

References assignConstituents, bgMembershipFunction, buildGridAlg(), fftjetcms::FFTJetInterface::checkConfig(), fftjetcms::FFTJetInterface::energyFlow, reco::tau::disc::Eta(), Exception, fftjetcms::fftjet_Grid2d_parser(), edm::ParameterSet::getParameter(), isCrisp, jetDistanceCalc, jetMembershipFunction, maxIterations, memberFactorCalcJet, memberFactorCalcPeak, myConfiguration, NULL, parse_bgMembershipFunction(), parse_jetDistanceCalc(), parse_jetMembershipFunction(), parse_memberFactorCalcJet(), parse_memberFactorCalcPeak(), parse_peakSelector(), parse_recoScaleCalcJet(), parse_recoScaleCalcPeak(), parse_recoScaleRatioCalcJet(), parse_recoScaleRatioCalcPeak(), peakSelector, colinearityKinematic::Phi, recoAlg, recombinationAlgorithm, recoScaleCalcJet, recoScaleCalcPeak, recoScaleRatioCalcJet, recoScaleRatioCalcPeak, reuseExistingGrid, unlikelyBgWeight, and useGriddedAlgorithm.

{
    const edm::ParameterSet& ps(myConfiguration);

    // Parse the peak selector definition
    peakSelector = parse_peakSelector(ps);
    checkConfig(peakSelector, "invalid peak selector");

    jetMembershipFunction = parse_jetMembershipFunction(ps);
    checkConfig(jetMembershipFunction, "invalid jet membership function");

    bgMembershipFunction = parse_bgMembershipFunction(ps);
    checkConfig(bgMembershipFunction, "invalid noise membership function");

    // Build the energy recombination algorithm
    if (!useGriddedAlgorithm)
    {
        fftjet::DefaultVectorRecombinationAlgFactory<
            VectorLike,BgData,VBuilder> factory;
        if (factory[recombinationAlgorithm] == NULL)
            throw cms::Exception("FFTJetBadConfig")
                << "Invalid vector recombination algorithm \""
                << recombinationAlgorithm << "\"" << std::endl;
        recoAlg = std::auto_ptr<RecoAlg>(
            factory[recombinationAlgorithm]->create(
                jetMembershipFunction.get(),
                &VectorLike::Et, &VectorLike::Eta, &VectorLike::Phi,
                bgMembershipFunction.get(),
                unlikelyBgWeight, isCrisp, false, assignConstituents));
    }
    else if (!reuseExistingGrid)
    {
        energyFlow = fftjet_Grid2d_parser(
            ps.getParameter<edm::ParameterSet>("GridConfiguration"));
        checkConfig(energyFlow, "invalid discretization grid");
        buildGridAlg();
    }

    // Parse the calculator of the recombination scale
    recoScaleCalcPeak = parse_recoScaleCalcPeak(ps);
    checkConfig(recoScaleCalcPeak, "invalid spec for the "
                "reconstruction scale calculator from peaks");

    // Parse the calculator of the recombination scale ratio
    recoScaleRatioCalcPeak = parse_recoScaleRatioCalcPeak(ps);
    checkConfig(recoScaleRatioCalcPeak, "invalid spec for the "
                "reconstruction scale ratio calculator from peaks");

    // Calculator for the membership function factor
    memberFactorCalcPeak = parse_memberFactorCalcPeak(ps);
    checkConfig(memberFactorCalcPeak, "invalid spec for the "
                "membership function factor calculator from peaks");

    if (maxIterations > 1)
    {
        // We are going to run iteratively. Make required objects.
        recoScaleCalcJet = parse_recoScaleCalcJet(ps);
        checkConfig(recoScaleCalcJet, "invalid spec for the "
                    "reconstruction scale calculator from jets");

        recoScaleRatioCalcJet = parse_recoScaleRatioCalcJet(ps);
        checkConfig(recoScaleRatioCalcJet, "invalid spec for the "
                    "reconstruction scale ratio calculator from jets");

        memberFactorCalcJet = parse_memberFactorCalcJet(ps);
        checkConfig(memberFactorCalcJet, "invalid spec for the "
                    "membership function factor calculator from jets");

        jetDistanceCalc = parse_jetDistanceCalc(ps);
        checkConfig(memberFactorCalcJet, "invalid spec for the "
                    "jet distance calculator");
    }
}

void FFTJetProducer::buildGridAlg ( ) [private]

Definition at line 361 of file FFTJetProducer.cc.

References assignConstituents, bgMembershipFunction, fftjetcms::FFTJetInterface::energyFlow, Exception, gridAlg, gridScanMaxEta, isCrisp, jetMembershipFunction, NULL, recombinationAlgorithm, recombinationDataCutoff, and unlikelyBgWeight.

Referenced by beginJob(), and loadEnergyFlow().

{
    int minBin = energyFlow->getEtaBin(-gridScanMaxEta);
    if (minBin < 0)
        minBin = 0;
    int maxBin = energyFlow->getEtaBin(gridScanMaxEta) + 1;
    if (maxBin < 0)
        maxBin = 0;
    
    fftjet::DefaultRecombinationAlgFactory<
        Real,VectorLike,BgData,VBuilder> factory;
    if (factory[recombinationAlgorithm] == NULL)
        throw cms::Exception("FFTJetBadConfig")
            << "Invalid grid recombination algorithm \""
            << recombinationAlgorithm << "\"" << std::endl;
    gridAlg = std::auto_ptr<GridAlg>(
        factory[recombinationAlgorithm]->create(
            jetMembershipFunction.get(),
            bgMembershipFunction.get(),
            unlikelyBgWeight, recombinationDataCutoff,
            isCrisp, false, assignConstituents, minBin, maxBin));
}

bool FFTJetProducer::checkConvergence	(	const std::vector< RecoFFTJet > &	previousIterResult,
		std::vector< RecoFFTJet > &	thisIterResult
	)		`[private]`

Definition at line 412 of file FFTJetProducer.cc.

References convergenceDistance, i, and jetDistanceCalc.

Referenced by iterateJetReconstruction().

{
    fftjet::Functor2<double,RecoFFTJet,RecoFFTJet>&
        distanceCalc(*jetDistanceCalc);

    const unsigned nJets = previous.size();
    const RecoFFTJet* prev = &previous[0];
    RecoFFTJet* next = &nextSet[0];

    // Calculate convergence distances for all jets
    bool converged = true;
    for (unsigned i=0; i<nJets; ++i)
    {
        const double d = distanceCalc(prev[i], next[i]);
        next[i].setConvergenceDistance(d);
        if (i < nJetsRequiredToConverge && d > convergenceDistance)
            converged = false;
    }

    return converged;
}

void FFTJetProducer::determineGriddedConstituents ( ) [private]

Definition at line 504 of file FFTJetProducer.cc.

References fftjetcms::FFTJetInterface::candidateIndex, constituents, fftjetcms::FFTJetInterface::energyFlow, fftjetcms::FFTJetInterface::eventData, g, gridAlg, i, fftjetcms::FFTJetInterface::inputCollection, and recoJets.

Referenced by produce().

{
    const unsigned nJets = recoJets.size();
    const unsigned* clusterMask = gridAlg->getClusterMask();
    const int nEta = gridAlg->getLastNEta();
    const int nPhi = gridAlg->getLastNPhi();
    const fftjet::Grid2d<Real>& g(*energyFlow);

    const unsigned nInputs = eventData.size();
    const VectorLike* inp = nInputs ? &eventData[0] : 0;
    const unsigned* candIdx = nInputs ? &candidateIndex[0] : 0;
    for (unsigned i=0; i<nInputs; ++i)
    {
        const VectorLike& item(inp[i]);
        const int iPhi = g.getPhiBin(item.Phi());
        const int iEta = g.getEtaBin(item.Eta());
        const unsigned mask = iEta >= 0 && iEta < nEta ?
            clusterMask[iEta*nPhi + iPhi] : 0;
        assert(mask <= nJets);
        constituents[mask].push_back(inputCollection->ptrAt(candIdx[i]));
    }
}

void FFTJetProducer::determineVectorConstituents ( ) [private]

Definition at line 528 of file FFTJetProducer.cc.

References fftjetcms::FFTJetInterface::candidateIndex, constituents, fftjetcms::FFTJetInterface::eventData, i, fftjetcms::FFTJetInterface::inputCollection, recoAlg, and recoJets.

Referenced by produce().

{
    const unsigned nJets = recoJets.size();
    const unsigned* clusterMask = recoAlg->getClusterMask();
    const unsigned maskLength = recoAlg->getLastNData();
    assert(maskLength == eventData.size());

    const unsigned* candIdx = maskLength ? &candidateIndex[0] : 0;
    for (unsigned i=0; i<maskLength; ++i)
    {
        // In FFTJet, the mask value of 0 corresponds to unclustered
        // energy. We will do the same here. Jet numbers are therefore
        // shifted by 1 wrt constituents vector, and constituents[1]
        // corresponds to jet number 0.
        const unsigned mask = clusterMask[i];
        assert(mask <= nJets);
        constituents[mask].push_back(inputCollection->ptrAt(candIdx[i]));
    }
}

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

Reimplemented from edm::EDProducer.

Definition at line 886 of file FFTJetProducer.cc.

{
}

unsigned FFTJetProducer::iterateJetReconstruction ( ) [private]

Definition at line 436 of file FFTJetProducer.cc.

References checkConvergence(), fftjetcms::FFTJetInterface::energyFlow, fftjetcms::FFTJetInterface::eventData, Exception, gridAlg, i, iterJets, iterPreclusters, metsig::jet, analyzePatCleaning_cfg::jets, maxIterations, memberFactorCalcJet, noiseLevel, L1TEmulatorMonitor_cff::p, preclusters, recoAlg, recoJets, recoScaleCalcJet, recoScaleRatioCalcJet, ntuplemaker::status, unclustered, unused, and useGriddedAlgorithm.

Referenced by produce().

{
    fftjet::Functor1<double,RecoFFTJet>& scaleCalc(*recoScaleCalcJet);
    fftjet::Functor1<double,RecoFFTJet>& ratioCalc(*recoScaleRatioCalcJet);
    fftjet::Functor1<double,RecoFFTJet>& factorCalc(*memberFactorCalcJet);

    const unsigned nJets = recoJets.size();

    unsigned iterNum = 1U;
    bool converged = false;
    for (; iterNum<maxIterations && !converged; ++iterNum)
    {
        // Recreate the vector of preclusters using the jets
        const RecoFFTJet* jets = &recoJets[0];
        iterPreclusters.clear();
        iterPreclusters.reserve(nJets);
        for (unsigned i=0; i<nJets; ++i)
        {
            const RecoFFTJet& jet(jets[i]);
            fftjet::Peak p(jet.precluster());
            p.setEtaPhi(jet.vec().Eta(), jet.vec().Phi());
            p.setRecoScale(scaleCalc(jet));
            p.setRecoScaleRatio(ratioCalc(jet));
            p.setMembershipFactor(factorCalc(jet));
            iterPreclusters.push_back(p);
        }

        // Run the algorithm
        int status = 0;
        if (useGriddedAlgorithm)
            status = gridAlg->run(iterPreclusters, *energyFlow,
                                  &noiseLevel, 1U, 1U,
                                  &iterJets, &unclustered, &unused);
        else
            status = recoAlg->run(iterPreclusters, eventData, &noiseLevel, 1U,
                                  &iterJets, &unclustered, &unused);
        if (status)
            throw cms::Exception("FFTJetInterface")
                << "FFTJet algorithm failed" << std::endl;
        assert(iterJets.size() == nJets);

        // Figure out if the iterations have converged
        converged = checkConvergence(recoJets, iterJets);

        // Prepare for the next cycle
        iterJets.swap(recoJets);
    }

    // Plug in the original precluster coordinates into the result
    assert(preclusters.size() == nJets);
    RecoFFTJet* jets = &recoJets[0];
    for (unsigned i=0; i<nJets; ++i)
    {
        const fftjet::Peak& oldp(preclusters[i]);
        jets[i].setPeakEtaPhi(oldp.eta(), oldp.phi());
    }

    // If we have converged on the last cycle, the result
    // would be indistinguishable from no convergence.
    // Because of this, raise the counter by one to indicate
    // the case when the convergence is not achieved.
    if (!converged)
        ++iterNum;

    return iterNum;
}

void FFTJetProducer::loadEnergyFlow ( const edm::Event & iEvent ) [private]

Definition at line 385 of file FFTJetProducer.cc.

References buildGridAlg(), fftjetcms::FFTJetInterface::energyFlow, edm::Event::getByLabel(), collect_tpl::input, NULL, and treeLabel.

Referenced by produce().

{
    edm::Handle<DiscretizedEnergyFlow> input;
    iEvent.getByLabel(treeLabel, input);

    // Make sure that the grid is compatible with the stored one
    bool rebuildGrid = energyFlow.get() == NULL;
    if (!rebuildGrid)
        rebuildGrid = 
            !(energyFlow->nEta() == input->nEtaBins() &&
              energyFlow->nPhi() == input->nPhiBins() &&
              energyFlow->etaMin() == input->etaMin() &&
              energyFlow->etaMax() == input->etaMax() &&
              energyFlow->phiBin0Edge() == input->phiBin0Edge());
    if (rebuildGrid)
    {
        // We should not get here very often...
        energyFlow = std::auto_ptr<fftjet::Grid2d<Real> >(
            new fftjet::Grid2d<Real>(
                input->nEtaBins(), input->etaMin(), input->etaMax(),
                input->nPhiBins(), input->phiBin0Edge(), input->title()));
        buildGridAlg();
    }
    energyFlow->blockSet(input->data(), input->nEtaBins(), input->nPhiBins());
}

template<class Real >

void FFTJetProducer::loadSparseTreeData ( const edm::Event & iEvent ) [private]

Definition at line 188 of file FFTJetProducer.cc.

References edm::Event::getByLabel(), fftjetcms::FFTJetInterface::getEventScale(), iniScales, collect_tpl::input, fftjetcms::sparsePeakTreeFromStorable(), sparseTree, and treeLabel.

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

    // Get the input
    edm::Handle<StoredTree> input;
    iEvent.getByLabel(treeLabel, input);

    if (!input->isSparse())
        throw cms::Exception("FFTJetBadConfig") 
            << "The stored clustering tree is not sparse" << std::endl;

    sparsePeakTreeFromStorable(*input, iniScales.get(), getEventScale(), &sparseTree);
    sparseTree.sortNodes();
}

template<typename T >

void FFTJetProducer::makeProduces	(	const std::string &	alias,
		const std::string &	tag
	)		`[private]`

Definition at line 114 of file FFTJetProducer.cc.

References GlobalPosition_Frontier_DevDB_cff::tag.

Referenced by FFTJetProducer().

{
    produces<std::vector<reco::FFTAnyJet<T> > >(tag).setBranchAlias(alias);
}

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

std::auto_ptr< AbsBgFunctor > FFTJetProducer::parse_bgMembershipFunction ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 737 of file FFTJetProducer.cc.

References fftjetcms::fftjet_BgFunctor_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_BgFunctor_parser(
        ps.getParameter<edm::ParameterSet>("bgMembershipFunction"));
}

std::auto_ptr< fftjet::Functor2< double, FFTJetProducer::RecoFFTJet, FFTJetProducer::RecoFFTJet > > FFTJetProducer::parse_jetDistanceCalc ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 797 of file FFTJetProducer.cc.

References fftjetcms::fftjet_JetDistance_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_JetDistance_parser(
        ps.getParameter<edm::ParameterSet>("jetDistanceCalc"));
}

std::auto_ptr< fftjet::ScaleSpaceKernel > FFTJetProducer::parse_jetMembershipFunction ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 728 of file FFTJetProducer.cc.

References fftjetcms::fftjet_MembershipFunction_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_MembershipFunction_parser(
        ps.getParameter<edm::ParameterSet>("jetMembershipFunction"));
}

std::auto_ptr< fftjet::Functor1< double, FFTJetProducer::RecoFFTJet > > FFTJetProducer::parse_memberFactorCalcJet ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 788 of file FFTJetProducer.cc.

References fftjetcms::fftjet_JetFunctor_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_JetFunctor_parser(
        ps.getParameter<edm::ParameterSet>("memberFactorCalcJet"));
}

std::auto_ptr< fftjet::Functor1< double, fftjet::Peak > > FFTJetProducer::parse_memberFactorCalcPeak ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 764 of file FFTJetProducer.cc.

References fftjetcms::fftjet_PeakFunctor_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_PeakFunctor_parser(
        ps.getParameter<edm::ParameterSet>("memberFactorCalcPeak"));
}

std::auto_ptr< fftjet::Functor1< bool, fftjet::Peak > > FFTJetProducer::parse_peakSelector ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 719 of file FFTJetProducer.cc.

References fftjetcms::fftjet_PeakSelector_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_PeakSelector_parser(
        ps.getParameter<edm::ParameterSet>("PeakSelectorConfiguration"));
}

std::auto_ptr< fftjet::Functor1< double, FFTJetProducer::RecoFFTJet > > FFTJetProducer::parse_recoScaleCalcJet ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 772 of file FFTJetProducer.cc.

References fftjetcms::fftjet_JetFunctor_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_JetFunctor_parser(
        ps.getParameter<edm::ParameterSet>("recoScaleCalcJet"));
}

std::auto_ptr< fftjet::Functor1< double, fftjet::Peak > > FFTJetProducer::parse_recoScaleCalcPeak ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 746 of file FFTJetProducer.cc.

References fftjetcms::fftjet_PeakFunctor_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_PeakFunctor_parser(
        ps.getParameter<edm::ParameterSet>("recoScaleCalcPeak"));
}

std::auto_ptr< fftjet::Functor1< double, FFTJetProducer::RecoFFTJet > > FFTJetProducer::parse_recoScaleRatioCalcJet ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 780 of file FFTJetProducer.cc.

References fftjetcms::fftjet_JetFunctor_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_JetFunctor_parser(
        ps.getParameter<edm::ParameterSet>("recoScaleRatioCalcJet"));
}

std::auto_ptr< fftjet::Functor1< double, fftjet::Peak > > FFTJetProducer::parse_recoScaleRatioCalcPeak ( const edm::ParameterSet & ps ) [protected, virtual]

Definition at line 755 of file FFTJetProducer.cc.

References fftjetcms::fftjet_PeakFunctor_parser(), and edm::ParameterSet::getParameter().

Referenced by beginJob().

{
    return fftjet_PeakFunctor_parser(
        ps.getParameter<edm::ParameterSet>("recoScaleRatioCalcPeak"));
}

FFTJetProducer::Resolution FFTJetProducer::parse_resolution ( const std::string & name ) [static]

Definition at line 95 of file FFTJetProducer.cc.

References Exception.

{
    if (!name.compare("fixed"))
        return FIXED;
    else if (!name.compare("maximallyStable"))
        return MAXIMALLY_STABLE;
    else if (!name.compare("globallyAdaptive"))
        return GLOBALLY_ADAPTIVE;
    else if (!name.compare("locallyAdaptive"))
        return LOCALLY_ADAPTIVE;
    else
        throw cms::Exception("FFTJetBadConfig")
            << "Invalid resolution specification \""
            << name << "\"" << std::endl;
}

void FFTJetProducer::prepareRecombinationScales ( ) [private]

Definition at line 338 of file FFTJetProducer.cc.

References i, memberFactorCalcPeak, preclusters, recoScaleCalcPeak, and recoScaleRatioCalcPeak.

Referenced by produce().

{
    const unsigned nClus = preclusters.size();
    if (nClus)
    {
        fftjet::Peak* clus = &preclusters[0];
        fftjet::Functor1<double,fftjet::Peak>& 
            scaleCalc(*recoScaleCalcPeak);
        fftjet::Functor1<double,fftjet::Peak>& 
            ratioCalc(*recoScaleRatioCalcPeak);
        fftjet::Functor1<double,fftjet::Peak>& 
            factorCalc(*memberFactorCalcPeak);

        for (unsigned i=0; i<nClus; ++i)
        {
            clus[i].setRecoScale(scaleCalc(clus[i]));
            clus[i].setRecoScaleRatio(ratioCalc(clus[i]));
            clus[i].setMembershipFactor(factorCalc(clus[i]));
        }
    }
}

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

Implements edm::EDProducer.

Definition at line 642 of file FFTJetProducer.cc.

References assignConstituents, assignMembershipFunctions(), constituents, determineGriddedConstituents(), determineVectorConstituents(), fftjetcms::FFTJetInterface::discretizeEnergyFlow(), fftjetcms::FFTJetInterface::energyFlow, fftjetcms::FFTJetInterface::eventData, Exception, gridAlg, i, iEvent, iterateJetReconstruction(), iterationsPerformed, loadEnergyFlow(), fftjetcms::FFTJetInterface::loadInputCollection(), maxIterations, noiseLevel, occupancy, peakSelector, preclusters, prepareRecombinationScales(), recoAlg, recoJets, reuseExistingGrid, saveResults(), selectPreclusters(), sparseTree, ntuplemaker::status, fftjetcms::FFTJetInterface::storeInSinglePrecision(), unclustered, unused, and useGriddedAlgorithm.

{
    // Load the clustering tree made by FFTJetPatRecoProducer
    if (storeInSinglePrecision())
        loadSparseTreeData<float>(iEvent);
    else
        loadSparseTreeData<double>(iEvent);

    // Do we need to load the candidate collection?
    if (assignConstituents || !(useGriddedAlgorithm && reuseExistingGrid))
        loadInputCollection(iEvent);

    // Do we need to have discretized energy flow?
    if (useGriddedAlgorithm)
    {
        if (reuseExistingGrid)
            loadEnergyFlow(iEvent);
        else
            discretizeEnergyFlow();
    }

    // Calculate cluster occupancy as a function of level number
    sparseTree.occupancyInScaleSpace(*peakSelector, &occupancy);

    // Select the preclusters using the requested resolution scheme
    preclusters.clear();
    selectPreclusters(sparseTree, *peakSelector, &preclusters);

    // Prepare to run the jet recombination procedure
    prepareRecombinationScales();

    // Assign membership functions to preclusters. If this function
    // is not overriden in a derived class, default algorithm membership
    // function will be used for every cluster.
    assignMembershipFunctions(&preclusters);

    // Run the recombination algorithm once
    int status = 0;
    if (useGriddedAlgorithm)
        status = gridAlg->run(preclusters, *energyFlow,
                              &noiseLevel, 1U, 1U,
                              &recoJets, &unclustered, &unused);
    else
        status = recoAlg->run(preclusters, eventData, &noiseLevel, 1U,
                              &recoJets, &unclustered, &unused);
    if (status)
        throw cms::Exception("FFTJetInterface")
            << "FFTJet algorithm failed (first iteration)" << std::endl;

    // If requested, iterate the jet recombination procedure
    if (maxIterations > 1U && !recoJets.empty())
        iterationsPerformed = iterateJetReconstruction();
    else
        iterationsPerformed = 1U;

    // Determine jet constituents. FFTJet returns a map
    // of constituents which is inverse to what we need here.
    const unsigned nJets = recoJets.size();
    if (constituents.size() <= nJets)
        constituents.resize(nJets + 1U);
    if (assignConstituents)
    {
        for (unsigned i=0; i<=nJets; ++i)
            constituents[i].clear();
        if (useGriddedAlgorithm)
            determineGriddedConstituents();
        else
            determineVectorConstituents();
    }

    // Write out the results
    saveResults(iEvent, iSetup);
}

void FFTJetProducer::saveResults	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)		`[private]`

Definition at line 212 of file SideBandSubtraction.cc.

References benchmark_cfg::cerr, i, combine::key, NULL, VarParsing::obj, convertSQLitetoXML_cfg::output, and stringify().

Referenced by produce().

{
  //saves the ephemeral stuff to a root file for future
  //use/modification (ie everything printed by printResults())

  TFile output(outname.c_str(),"UPDATE"); //open the output file,
                                          //create it if it doesn't
                                          //exist
  //Since keys are only available from files on disk, we need to write
  //out a new file.  If the file already existed, then we opened to
  //update, and are writing nothing new.  
  output.Write();
  TString dirname;
  TIter nextkey(output.GetListOfKeys());
  TKey *key;
  TDirectory* curDir=NULL;
  while((key=(TKey*)nextkey.Next()))
    {

      if(key==NULL)
        break;
      TObject *obj = key->ReadObj();
      if(obj->IsA()->InheritsFrom("TDirectory"))
          dirname=obj->GetName();

    }

  if(dirname=="")
    {
      //we didn't find any directories so, we'll make a new one
      curDir = output.mkdir("run0","Run 0");
      output.cd("run0");
    }
  else
    {
      //manipulate last dir string, make a new one, and get ready to fill
      dirname.Remove(0,3);
      Int_t run_num = dirname.Atoi();
      run_num++;
      dirname = "run" + stringify(run_num);
      curDir = output.mkdir(dirname.Data(),("Run "+stringify(run_num)).c_str());
      output.cd(dirname.Data());
    }
  if(curDir==NULL)
    curDir = output.GetDirectory("",kTRUE,"");

  //these should all be the same size, but to be pedantic we'll loop
  //over each one individually, also, we need to associate them with
  //the directory because by default they "float" in memory to avoid
  //conflicts with other root files the user has open. If they want to
  //write to those files, they need to close their file, pass the name
  //here, and then let us work.
  for(unsigned int i=0; i < SideBandHistos.size(); ++i)
    {
      SideBandHistos[i].SetDirectory(curDir);
      SideBandHistos[i].Write();
    }
  for(unsigned int i=0; i < RawHistos.size(); ++i)
    {
      RawHistos[i].SetDirectory(curDir);
      RawHistos[i].Write();
    }
  for(unsigned int i=0; i < SBSHistos.size(); ++i)
    {
      SBSHistos[i].SetDirectory(curDir);
      SBSHistos[i].Write();
    }
  if(Data!=NULL && ModelPDF!=NULL && BackgroundPDF!=NULL && SeparationVariable!=NULL)
    {
      RooPlot *sep_varFrame = SeparationVariable->frame();
      Data->plotOn(sep_varFrame);
      ModelPDF->plotOn(sep_varFrame);
      BackgroundPDF->plotOn(sep_varFrame);
      sep_varFrame->Write();
    }
  else
    cerr <<"ERROR: saveResults, did not save RooPlot of data and fit\n";
  output.Write();

}

void FFTJetProducer::selectPreclusters	(	const SparseTree &	tree,
		const fftjet::Functor1< bool, fftjet::Peak > &	peakSelector,
		std::vector< fftjet::Peak > *	preclusters
	)		`[protected, virtual]`

Definition at line 205 of file FFTJetProducer.cc.

References i, nodes, resolution, selectTreeNodes(), and sparseTree.

Referenced by produce().

{
    nodes.clear();
    selectTreeNodes(tree, peakSelector, &nodes);

    // Fill out the vector of preclusters using the tree node ids
    const unsigned nNodes = nodes.size();
    const SparseTree::NodeId* pnodes = nNodes ? &nodes[0] : 0;
    preclusters->reserve(nNodes);
    for (unsigned i=0; i<nNodes; ++i)
        preclusters->push_back(
            sparseTree.uncheckedNode(pnodes[i]).getCluster());

    // Set the status word to indicate the resolution scheme used
    fftjet::Peak* clusters = nNodes ? &(*preclusters)[0] : 0;
    for (unsigned i=0; i<nNodes; ++i)
        clusters[i].setStatus(resolution);
}

void FFTJetProducer::selectTreeNodes	(	const SparseTree &	tree,
		const fftjet::Functor1< bool, fftjet::Peak > &	peakSelect,
		std::vector< SparseTree::NodeId > *	nodes
	)		`[protected]`

Definition at line 228 of file FFTJetProducer.cc.

References delta, FIXED, fixedScale, GLOBALLY_ADAPTIVE, testEve_cfg::level, LOCALLY_ADAPTIVE, MAXIMALLY_STABLE, maxLevel, maxStableScale, minLevel, minStableScale, n, nClustersRequested, occupancy, resolution, stabilityAlpha, TopDecayID::stable, thresholds, and usedLevel.

Referenced by selectPreclusters().

{
    minLevel = maxLevel = usedLevel = 0;

    // Get the tree nodes which pass the cuts
    // (according to the selected resolution strategy)
    switch (resolution)
    {
    case FIXED:
    {
        usedLevel = tree.getLevel(fixedScale);
        tree.getPassingNodes(usedLevel, peakSelect, mynodes);
    }
    break;

    case MAXIMALLY_STABLE:
    {
        const unsigned minStartingLevel = maxStableScale > 0.0 ? 
            tree.getLevel(maxStableScale) : 0;
        const unsigned maxStartingLevel = minStableScale > 0.0 ?
            tree.getLevel(minStableScale) : UINT_MAX;

        if (tree.stableClusterCount(
                peakSelect, &minLevel, &maxLevel, stabilityAlpha,
                minStartingLevel, maxStartingLevel))
        {
            usedLevel = (minLevel + maxLevel)/2;
            tree.getPassingNodes(usedLevel, peakSelect, mynodes);
        }
    }
    break;

    case GLOBALLY_ADAPTIVE:
    {
        const bool stable = tree.clusterCountLevels(
            nClustersRequested, peakSelect, &minLevel, &maxLevel);
        if (minLevel || maxLevel)
        {
            usedLevel = (minLevel + maxLevel)/2;
            if (!stable)
            {
                const int maxlev = tree.maxStoredLevel();
                bool levelFound = false;
                for (int delta=0; delta<=maxlev && !levelFound; ++delta)
                    for (int ifac=1; ifac>-2 && !levelFound; ifac-=2)
                    {
                        const int level = usedLevel + ifac*delta;
                        if (level > 0 && level <= maxlev)
                            if (occupancy[level] == nClustersRequested)
                            {
                                usedLevel = level;
                                levelFound = true;
                            }
                    }
                assert(levelFound);
            }
        }
        else
        {
            // Can't find that exact number of preclusters.
            // Try to get the next best thing.
            usedLevel = 1;
            const unsigned occ1 = occupancy[1];
            if (nClustersRequested >= occ1)
            {
                const unsigned maxlev = tree.maxStoredLevel();
                if (nClustersRequested > occupancy[maxlev])
                    usedLevel = maxlev;
                else
                {
                    // It would be nice to use "lower_bound" here,
                    // but the occupancy is not necessarily monotonous.
                    unsigned bestDelta = nClustersRequested > occ1 ?
                        nClustersRequested - occ1 : occ1 - nClustersRequested;
                    for (unsigned level=2; level<=maxlev; ++level)
                    {
                        const unsigned n = occupancy[level];
                        const unsigned d = nClustersRequested > n ? 
                            nClustersRequested - n : n - nClustersRequested;
                        if (d < bestDelta)
                        {
                            bestDelta = d;
                            usedLevel = level;
                        }
                    }
                }
            }
        }
        tree.getPassingNodes(usedLevel, peakSelect, mynodes);
    }
    break;

    case LOCALLY_ADAPTIVE:
    {
        usedLevel = tree.getLevel(fixedScale);
        tree.getMagS2OptimalNodes(peakSelect, nClustersRequested,
                                  usedLevel, mynodes, &thresholds);
    }
    break;

    default:
        assert(!"ERROR in FFTJetProducer::selectTreeNodes : "
               "should never get here! This is a bug. Please report.");
    }
}

template<typename T >

void FFTJetProducer::writeJets	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)		`[private]`

Definition at line 552 of file FFTJetProducer.cc.

References constituents, fftjetcms::FFTJetInterface::energyFlow, i, metsig::jet, analyzePatCleaning_cfg::jets, fftjetcms::FFTJetInterface::outputLabel, p4, edm::Event::put(), dt_dqm_sourceclient_common_cff::reco, recoJets, recombinationDataCutoff, resumConstituents, useGriddedAlgorithm, fftjetcms::FFTJetInterface::vertexUsed(), and reco::writeSpecific().

{
    using namespace reco;

    typedef FFTAnyJet<T> OutputJet;
    typedef std::vector<OutputJet> OutputCollection;

    // Area of a single eta-phi cell for jet area calculations.
    // Set it to 0 in case the module configuration does not allow
    // us to calculate jet areas reliably.
    const double cellArea = useGriddedAlgorithm && 
                            recombinationDataCutoff < 0.0 ?
        energyFlow->etaBinWidth() * energyFlow->phiBinWidth() : 0.0;

    // allocate output jet collection
    std::auto_ptr<OutputCollection> jets(new OutputCollection());
    const unsigned nJets = recoJets.size();
    jets->reserve(nJets);

    for (unsigned ijet=0; ijet<nJets; ++ijet)
    {
        const RecoFFTJet& myjet(recoJets[ijet]);

        // Check if we should resum jet constituents
        VectorLike jet4vec(myjet.vec());
        if (resumConstituents)
        {
            VectorLike sum(0.0, 0.0, 0.0, 0.0);
            const unsigned nCon = constituents[ijet+1].size();
            const reco::CandidatePtr* cn = nCon ? &constituents[ijet+1][0] : 0;
            for (unsigned i=0; i<nCon; ++i)
                sum += cn[i]->p4();
            jet4vec = sum;
        }

        // Write the specifics to the jet (simultaneously sets 4-vector,
        // vertex, constituents). These are overridden functions that will
        // call the appropriate specific code.
        T jet;
        writeSpecific(jet, jet4vec, vertexUsed(),
                      constituents[ijet+1], iSetup);

        // calcuate the jet area
        jet.setJetArea(cellArea*myjet.ncells());

        // add jet to the list
        jets->push_back(OutputJet(jet, jetToStorable<float>(myjet)));
    }

    // put the collection into the event
    iEvent.put(jets, outputLabel);
}