#include <VirtualJetProducer.h>

Inheritance diagram for VirtualJetProducer:

Classes
struct	JetType

Public Types
typedef boost::shared_ptr < fastjet::GhostedAreaSpec >	ActiveAreaSpecPtr

typedef boost::shared_ptr < fastjet::AreaDefinition >	AreaDefinitionPtr

typedef boost::shared_ptr < fastjet::ClusterSequence >	ClusterSequencePtr

typedef boost::shared_ptr < fastjet::JetDefinition >	JetDefPtr

typedef boost::shared_ptr < fastjet::JetDefinition::Plugin >	PluginPtr

typedef boost::shared_ptr < fastjet::RangeDefinition >	RangeDefPtr

Public Types inherited from edm::stream::EDProducer<>
typedef CacheContexts< T...>	CacheTypes

typedef CacheTypes::GlobalCache	GlobalCache

typedef AbilityChecker< T...>	HasAbility

typedef CacheTypes::LuminosityBlockCache	LuminosityBlockCache

typedef LuminosityBlockContextT < LuminosityBlockCache, RunCache, GlobalCache >	LuminosityBlockContext

typedef CacheTypes::LuminosityBlockSummaryCache	LuminosityBlockSummaryCache

typedef CacheTypes::RunCache	RunCache

typedef RunContextT< RunCache, GlobalCache >	RunContext

typedef CacheTypes::RunSummaryCache	RunSummaryCache

Public Types inherited from edm::stream::EDProducerBase
typedef EDProducerAdaptorBase	ModuleType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Public Member Functions
std::string	jetType () const

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

	VirtualJetProducer (const edm::ParameterSet &iConfig)

virtual	~VirtualJetProducer ()

Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default

Public Member Functions inherited from edm::stream::EDProducerBase
	EDProducerBase ()

ModuleDescription const &	moduleDescription () const

virtual	~EDProducerBase ()

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

	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
std::vector< ConsumesInfo >	consumesInfo () const

	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

void	modulesDependentUpon (std::string const &iProcessName, std::string const &iModuleLabel, bool iPrint, std::vector< char const * > &oModuleLabels) const

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

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

Protected Member Functions
virtual void	copyConstituents (const std::vector< fastjet::PseudoJet > &fjConstituents, reco::Jet *jet)

virtual std::vector < reco::CandidatePtr >	getConstituents (const std::vector< fastjet::PseudoJet > &fjConstituents)

virtual void	inputTowers ()

virtual bool	isAnomalousTower (reco::CandidatePtr input)

bool	makeBasicJet (const JetType::Type &fTag)

bool	makeCaloJet (const JetType::Type &fTag)

bool	makeGenJet (const JetType::Type &fTag)

bool	makePFClusterJet (const JetType::Type &fTag)

bool	makePFJet (const JetType::Type &fTag)

virtual void	makeProduces (std::string s, std::string tag="")

bool	makeTrackJet (const JetType::Type &fTag)

void	offsetCorrectJets (std::vector< fastjet::PseudoJet > &orphanInput)

virtual void	output (edm::Event &iEvent, edm::EventSetup const &iSetup)

virtual void	runAlgorithm (edm::Event &iEvent, const edm::EventSetup &iSetup)=0

template<typename T >
void	writeCompoundJets (edm::Event &iEvent, edm::EventSetup const &iSetup)
	function template to write out the outputs More...

template<typename T >
void	writeJets (edm::Event &iEvent, edm::EventSetup const &iSetup)

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 Attributes
bool	doAreaDiskApprox_

bool	doAreaFastjet_

bool	doFastJetNonUniform_

bool	doPUOffsetCorr_

bool	doPVCorrection_

bool	doRhoFastjet_

ActiveAreaSpecPtr	fjActiveArea_

AreaDefinitionPtr	fjAreaDefinition_

ClusterSequencePtr	fjClusterSeq_

std::vector< fastjet::PseudoJet >	fjInputs_

JetDefPtr	fjJetDefinition_

std::vector< fastjet::PseudoJet >	fjJets_

PluginPtr	fjPlugin_

RangeDefPtr	fjRangeDef_

edm::EDGetTokenT < reco::VertexCollection >	input_vertex_token_

double	inputEMin_

double	inputEtMin_

std::vector< edm::Ptr < reco::Candidate > >	inputs_

std::string	jetAlgorithm_

std::string	jetCollInstanceName_

double	jetPtMin_

std::string	jetType_

JetType::Type	jetTypeE

unsigned int	maxInputs_

unsigned int	minSeed_

std::string	moduleLabel_

unsigned int	nExclude_

std::vector< double >	puCenters_

std::string	puSubtractorName_

double	puWidth_

bool	restrictInputs_

double	rParam_

edm::InputTag	src_

edm::InputTag	srcPVs_

boost::shared_ptr < PileUpSubtractor >	subtractor_

bool	useDeterministicSeed_

bool	useExplicitGhosts_

int	verbosity_

reco::Particle::Point	vertex_

double	voronoiRfact_

bool	writeCompound_

Private Attributes
std::auto_ptr< AnomalousTower >	anomalousTowerDef_

edm::EDGetTokenT< std::vector < edm::FwdPtr < reco::PFCandidate > > >	input_candidatefwdptr_token_

edm::EDGetTokenT < reco::CandidateView >	input_candidateview_token_

edm::EDGetTokenT< std::vector < edm::FwdPtr < pat::PackedCandidate > > >	input_packedcandidatefwdptr_token_

Additional Inherited Members
Static Public Member Functions inherited from edm::stream::EDProducerBase
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Detailed Description

Definition at line 31 of file VirtualJetProducer.h.

Member Typedef Documentation

typedef boost::shared_ptr<fastjet::GhostedAreaSpec> VirtualJetProducer::ActiveAreaSpecPtr

Definition at line 84 of file VirtualJetProducer.h.

typedef boost::shared_ptr<fastjet::AreaDefinition> VirtualJetProducer::AreaDefinitionPtr

Definition at line 85 of file VirtualJetProducer.h.

typedef boost::shared_ptr<fastjet::ClusterSequence> VirtualJetProducer::ClusterSequencePtr

Definition at line 81 of file VirtualJetProducer.h.

typedef boost::shared_ptr<fastjet::JetDefinition> VirtualJetProducer::JetDefPtr

Definition at line 83 of file VirtualJetProducer.h.

typedef boost::shared_ptr<fastjet::JetDefinition::Plugin> VirtualJetProducer::PluginPtr

Definition at line 82 of file VirtualJetProducer.h.

typedef boost::shared_ptr<fastjet::RangeDefinition> VirtualJetProducer::RangeDefPtr

Definition at line 86 of file VirtualJetProducer.h.

Constructor & Destructor Documentation

VirtualJetProducer::VirtualJetProducer ( const edm::ParameterSet & iConfig )

explicit

Definition at line 115 of file VirtualJetProducer.cc.

References HLT_25ns10e33_v2_cff::alias, anomalousTowerDef_, VirtualJetProducer::JetType::byName(), edm::EDConsumerBase::consumesCollector(), doAreaDiskApprox_, doAreaFastjet_, doFastJetNonUniform_, doPUOffsetCorr_, doRhoFastjet_, Exception, edm::ParameterSet::exists(), fjActiveArea_, fjAreaDefinition_, fjJetDefinition_, fjPlugin_, fjRangeDef_, reco::get(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), input_candidatefwdptr_token_, input_candidateview_token_, input_packedcandidatefwdptr_token_, input_vertex_token_, jetAlgorithm_, jetCollInstanceName_, jetType_, jetTypeE, edm::InputTag::label(), makeProduces(), maxInputs_, minSeed_, moduleLabel_, nExclude_, puCenters_, puSubtractorName_, puWidth_, restrictInputs_, rParam_, src_, srcPVs_, AlCaHLTBitMon_QueryRunRegistry::string, subtractor_, useDeterministicSeed_, useExplicitGhosts_, verbosity_, voronoiRfact_, and writeCompound_.

   : moduleLabel_   (iConfig.getParameter<string>       ("@module_label"))
   , src_           (iConfig.getParameter<edm::InputTag>("src"))
   , srcPVs_        (iConfig.getParameter<edm::InputTag>("srcPVs"))
   , jetType_       (iConfig.getParameter<string>       ("jetType"))
   , jetAlgorithm_  (iConfig.getParameter<string>       ("jetAlgorithm"))
   , rParam_        (iConfig.getParameter<double>       ("rParam"))
   , inputEtMin_    (iConfig.getParameter<double>       ("inputEtMin"))
   , inputEMin_     (iConfig.getParameter<double>       ("inputEMin"))
   , jetPtMin_      (iConfig.getParameter<double>       ("jetPtMin"))
   , doPVCorrection_(iConfig.getParameter<bool>         ("doPVCorrection"))
   , restrictInputs_(false)
   , maxInputs_(99999999)
   , doAreaFastjet_ (iConfig.getParameter<bool>         ("doAreaFastjet"))
   , useExplicitGhosts_(false)
   , doAreaDiskApprox_       (false)
   , doRhoFastjet_  (iConfig.getParameter<bool>         ("doRhoFastjet"))
   , voronoiRfact_           (-9)
   , doPUOffsetCorr_(iConfig.getParameter<bool>         ("doPUOffsetCorr"))
   , puWidth_(0)
   , nExclude_(0)
   , jetCollInstanceName_ ("")
   , writeCompound_ ( false )
   , verbosity_(0)
 {
   anomalousTowerDef_ = std::auto_ptr<AnomalousTower>(new AnomalousTower(iConfig));
 
   //
   // additional parameters to think about:
   // - overlap threshold (set to 0.75 for the time being)
   // - p parameter for generalized kT (set to -2 for the time being)
   // - fastjet PU subtraction parameters (not yet considered)
   //
   if (jetAlgorithm_=="SISCone") {
     fjPlugin_ = PluginPtr( new fastjet::SISConePlugin(rParam_,0.75,0,0.0,false,
                                                       fastjet::SISConePlugin::SM_pttilde) );
     fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(&*fjPlugin_) );
   }
   else if (jetAlgorithm_=="IterativeCone") {
     fjPlugin_ = PluginPtr(new fastjet::CMSIterativeConePlugin(rParam_,1.0));
     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(&*fjPlugin_));
   }
   else if (jetAlgorithm_=="CDFMidPoint") {
     fjPlugin_ = PluginPtr(new fastjet::CDFMidPointPlugin(rParam_,0.75));
     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(&*fjPlugin_));
   }
   else if (jetAlgorithm_=="ATLASCone") {
     fjPlugin_ = PluginPtr(new fastjet::ATLASConePlugin(rParam_));
     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(&*fjPlugin_));
   }
   else if (jetAlgorithm_=="Kt")
     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(fastjet::kt_algorithm,rParam_));
   else if (jetAlgorithm_=="CambridgeAachen")
     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(fastjet::cambridge_algorithm,
                                                            rParam_) );
   else if (jetAlgorithm_=="AntiKt")
     fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::antikt_algorithm,rParam_) );
   else if (jetAlgorithm_=="GeneralizedKt")
     fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::genkt_algorithm,
                                                             rParam_,-2) );
   else
     throw cms::Exception("Invalid jetAlgorithm")
       <<"Jet algorithm for VirtualJetProducer is invalid, Abort!\n";
   
   jetTypeE=JetType::byName(jetType_);
 
   if ( iConfig.exists("jetCollInstanceName") ) {
     jetCollInstanceName_ = iConfig.getParameter<string>("jetCollInstanceName");
   }
 
   if ( doPUOffsetCorr_ ) {
 
      if(iConfig.exists("subtractorName")) puSubtractorName_  =  iConfig.getParameter<string> ("subtractorName");
      else puSubtractorName_ = std::string();
      
      if(puSubtractorName_.empty()){
        edm::LogWarning("VirtualJetProducer") << "Pile Up correction on; however, pile up type is not specified. Using default... \n";
        subtractor_ =  boost::shared_ptr<PileUpSubtractor>(new PileUpSubtractor(iConfig, consumesCollector()));
      }else{
        subtractor_ =  boost::shared_ptr<PileUpSubtractor>(PileUpSubtractorFactory::get()->create( puSubtractorName_, iConfig, consumesCollector()));
      }
   }
 
   // use explicit ghosts in the fastjet clustering sequence?
   if ( iConfig.exists("useExplicitGhosts") ) {
     useExplicitGhosts_ = iConfig.getParameter<bool>("useExplicitGhosts");
   }
 
   // do approximate disk-based area calculation => warn if conflicting request
   if (iConfig.exists("doAreaDiskApprox")) {
     doAreaDiskApprox_ = iConfig.getParameter<bool>("doAreaDiskApprox");
     if (doAreaDiskApprox_ && doAreaFastjet_)
       throw cms::Exception("Conflicting area calculations") << "Both the calculation of jet area via fastjet and via an analytical disk approximation have been requested. Please decide on one.\n";
 
   }
   // turn off jet collection output for speed
   // Voronoi-based area calculation allows for an empirical scale factor
   if (iConfig.exists("voronoiRfact"))
     voronoiRfact_     = iConfig.getParameter<double>("voronoiRfact");
 
 
   // do fasjet area / rho calcluation? => accept corresponding parameters
   if ( doAreaFastjet_ || doRhoFastjet_ ) {
     // Eta range of jets to be considered for Rho calculation
     // Should be at most (jet acceptance - jet radius)
     double rhoEtaMax=iConfig.getParameter<double>("Rho_EtaMax");
     // default Ghost_EtaMax should be 5
     double ghostEtaMax = iConfig.getParameter<double>("Ghost_EtaMax");
     // default Active_Area_Repeats 1
     int    activeAreaRepeats = iConfig.getParameter<int> ("Active_Area_Repeats");
     // default GhostArea 0.01
     double ghostArea = iConfig.getParameter<double> ("GhostArea");
     if (voronoiRfact_ <= 0) {
       fjActiveArea_     = ActiveAreaSpecPtr(new fastjet::GhostedAreaSpec(ghostEtaMax,activeAreaRepeats,ghostArea));
       fjActiveArea_->set_fj2_placement(true);
       if ( ! useExplicitGhosts_ ) {
     fjAreaDefinition_ = AreaDefinitionPtr( new fastjet::AreaDefinition(fastjet::active_area, *fjActiveArea_ ) );
       } else {
     fjAreaDefinition_ = AreaDefinitionPtr( new fastjet::AreaDefinition(fastjet::active_area_explicit_ghosts, *fjActiveArea_ ) );
       }
     }
     fjRangeDef_ = RangeDefPtr( new fastjet::RangeDefinition(rhoEtaMax) );
   } 
 
   // restrict inputs to first "maxInputs" towers?
   if ( iConfig.exists("restrictInputs") ) {
     restrictInputs_ = iConfig.getParameter<bool>("restrictInputs");
     maxInputs_      = iConfig.getParameter<unsigned int>("maxInputs");
   }
  
 
   string alias=iConfig.getUntrackedParameter<string>("alias",moduleLabel_);
 
 
   // Check to see if we are writing compound jets for substructure
   // and jet grooming
   if ( iConfig.exists("writeCompound") ) {
     writeCompound_ = iConfig.getParameter<bool>("writeCompound");
   }
 
   // make the "produces" statements
   makeProduces( alias, jetCollInstanceName_ );
 
   doFastJetNonUniform_ = false;
   if(iConfig.exists("doFastJetNonUniform")) doFastJetNonUniform_ = iConfig.getParameter<bool>   ("doFastJetNonUniform");
   if(doFastJetNonUniform_){
     puCenters_ = iConfig.getParameter<std::vector<double> >("puCenters");
     puWidth_ = iConfig.getParameter<double>("puWidth");
     nExclude_ = iConfig.getParameter<unsigned int>("nExclude");
   }
 
   useDeterministicSeed_ = false;
   minSeed_ = 0;
   if ( iConfig.exists("useDeterministicSeed") ) {
     useDeterministicSeed_ = iConfig.getParameter<bool>("useDeterministicSeed");
     minSeed_ =              iConfig.getParameter<unsigned int>("minSeed");
   }
   
   if ( iConfig.exists("verbosity" ) ) {
     verbosity_ = iConfig.getParameter<int>("verbosity");
   }
   
   produces<std::vector<double> >("rhos");
   produces<std::vector<double> >("sigmas");
   produces<double>("rho");
   produces<double>("sigma");
 
   if (!srcPVs_.label().empty()) input_vertex_token_ = consumes<reco::VertexCollection>(srcPVs_);
   input_candidateview_token_ = consumes<reco::CandidateView>(src_);
   input_candidatefwdptr_token_ = consumes<std::vector<edm::FwdPtr<reco::PFCandidate> > >(src_);
   input_packedcandidatefwdptr_token_ = consumes<std::vector<edm::FwdPtr<pat::PackedCandidate> > >(src_);
   
 }

VirtualJetProducer::~VirtualJetProducer ( )

virtual

Definition at line 290 of file VirtualJetProducer.cc.

291 {

292 }

Member Function Documentation

void VirtualJetProducer::copyConstituents	(	const std::vector< fastjet::PseudoJet > &	fjConstituents,
		reco::Jet *	jet
	)

protectedvirtual

Definition at line 504 of file VirtualJetProducer.cc.

References reco::CompositePtrCandidate::addDaughter(), i, cmsHarvester::index, and inputs_.

 {
   for (unsigned int i=0;i<fjConstituents.size();++i) { 
     int index = fjConstituents[i].user_index();
     if ( index >= 0 && static_cast<unsigned int>(index) < inputs_.size() )
       jet->addDaughter(inputs_[index]);
   }
 }

vector< reco::CandidatePtr > VirtualJetProducer::getConstituents ( const std::vector< fastjet::PseudoJet > & fjConstituents )

protectedvirtual

Definition at line 517 of file VirtualJetProducer.cc.

References i, cmsHarvester::index, inputs_, and mps_fire::result.

Referenced by FastjetJetProducer::produceTrackJets(), cms::SubEventGenJetProducer::runAlgorithm(), writeCompoundJets(), and writeJets().

 {
   vector<reco::CandidatePtr> result; result.reserve(fjConstituents.size()/2);
   for (unsigned int i=0;i<fjConstituents.size();i++) {
     auto index = fjConstituents[i].user_index();
     if ( index >= 0 && static_cast<unsigned int>(index) < inputs_.size() ) {
       result.emplace_back(inputs_[index]);
     }
   }
   return result;
 }

void VirtualJetProducer::inputTowers ( )

protectedvirtual

Reimplemented in cms::CompoundJetProducer, SubjetFilterJetProducer, and cms::SubEventGenJetProducer.

Definition at line 445 of file VirtualJetProducer.cc.

References doPVCorrection_, geometryDiff::epsilon, fjInputs_, i, input, inputEMin_, inputEtMin_, inputs_, isAnomalousTower(), edm::isNotFinite(), jetTypeE, makeCaloJet(), maxInputs_, CaloTower::p4(), restrictInputs_, and vertex_.

Referenced by SubjetFilterJetProducer::inputTowers(), cms::CompoundJetProducer::inputTowers(), produce(), and FastjetJetProducer::produceTrackJets().

 {
   auto inBegin = inputs_.begin(),
     inEnd = inputs_.end(), i = inBegin;
   for (; i != inEnd; ++i ) {
     auto const & input = **i;
     // std::cout << "CaloTowerVI jets " << input->pt() << " " << input->et() << ' '<< input->energy() << ' ' << (isAnomalousTower(input) ? " bad" : " ok") << std::endl; 
     if (edm::isNotFinite(input.pt()))           continue;
     if (input.et()    <inputEtMin_)  continue;
     if (input.energy()<inputEMin_)   continue;
     if (isAnomalousTower(*i))      continue;
     // Change by SRR : this is no longer an error nor warning, this can happen with PU mitigation algos.
     // Also switch to something more numerically safe. (VI: 10^-42GeV????)
     if (input.pt() < 100 * std::numeric_limits<double>::epsilon() ) { 
       continue;
     }
     if (makeCaloJet(jetTypeE)&&doPVCorrection_) {
       const CaloTower & tower = dynamic_cast<const CaloTower &>(input);
       auto const &  ct = tower.p4(vertex_);  // very expensive as computed in eta/phi
       fjInputs_.emplace_back(ct.px(),ct.py(),ct.pz(),ct.energy());
       //std::cout << "tower:" << *tower << '\n';
     }
     else {
       /*
       if(makePFJet(jetTypeE)) {
     reco::PFCandidate& pfc = (reco::PFCandidate&)input;
     std::cout << "PF cand:" << pfc << '\n';
       }
       */
       fjInputs_.emplace_back(input.px(),input.py(),input.pz(),
                          input.energy());
     }
     fjInputs_.back().set_user_index(i - inBegin);
   }
 
   if ( restrictInputs_ && fjInputs_.size() > maxInputs_ ) {
     reco::helper::GreaterByPtPseudoJet   pTComparator;
     std::sort(fjInputs_.begin(), fjInputs_.end(), pTComparator);
     fjInputs_.resize(maxInputs_);
     edm::LogWarning("JetRecoTooManyEntries") << "Too many inputs in the event, limiting to first " << maxInputs_ << ". Output is suspect.";
   }
 }

bool VirtualJetProducer::isAnomalousTower ( reco::CandidatePtr input )

protectedvirtual

Definition at line 489 of file VirtualJetProducer.cc.

References anomalousTowerDef_, jetTypeE, and makeCaloJet().

Referenced by cms::SubEventGenJetProducer::inputTowers(), and inputTowers().

 {
   if (!makeCaloJet(jetTypeE)) 
       return false;
   else
       return (*anomalousTowerDef_)(*input);
 }

std::string VirtualJetProducer::jetType ( ) const

inline

Definition at line 93 of file VirtualJetProducer.h.

93 { return jetType_; }

VirtualJetProducer::jetType_

std::string jetType_

Definition: VirtualJetProducer.h:154

bool VirtualJetProducer::makeBasicJet ( const JetType::Type & fTag )

inlineprotected

Definition at line 68 of file VirtualJetProducer.h.

                                                     {
     return fTag == JetType::BasicJet;
   }

bool VirtualJetProducer::makeCaloJet ( const JetType::Type & fTag )

inlineprotected

Definition at line 53 of file VirtualJetProducer.h.

Referenced by inputTowers(), isAnomalousTower(), and produce().

                                                    {
     return fTag == JetType::CaloJet;
   }

bool VirtualJetProducer::makeGenJet ( const JetType::Type & fTag )

inlineprotected

Definition at line 59 of file VirtualJetProducer.h.

                                                   {
     return fTag == JetType::GenJet;
   }

bool VirtualJetProducer::makePFClusterJet ( const JetType::Type & fTag )

inlineprotected

Definition at line 65 of file VirtualJetProducer.h.

                                                         {
     return fTag == JetType::PFClusterJet;
   }

bool VirtualJetProducer::makePFJet ( const JetType::Type & fTag )

inlineprotected

Definition at line 56 of file VirtualJetProducer.h.

Referenced by produce().

                                                  {
     return fTag == JetType::PFJet;
   }

void VirtualJetProducer::makeProduces	(	std::string	s,
		std::string	tag = `""`
	)

protectedvirtual

Definition at line 82 of file VirtualJetProducer.cc.

References GlobalPosition_Frontier_DevDB_cff::tag.

Referenced by SubjetFilterJetProducer::SubjetFilterJetProducer(), and VirtualJetProducer().

 {
 
 
   if ( writeCompound_ ) {
     produces<reco::BasicJetCollection>();
   }
 
   if (makeCaloJet(jetTypeE)) {
     produces<reco::CaloJetCollection>(tag).setBranchAlias(alias);
   }
   else if (makePFJet(jetTypeE)) {
     produces<reco::PFJetCollection>(tag).setBranchAlias(alias);
   }
   else if (makeGenJet(jetTypeE)) {
     produces<reco::GenJetCollection>(tag).setBranchAlias(alias);
   }
   else if (makeTrackJet(jetTypeE)) {
     produces<reco::TrackJetCollection>(tag).setBranchAlias(alias);
   }
   else if (makePFClusterJet(jetTypeE)) {
     produces<reco::PFClusterJetCollection>(tag).setBranchAlias(alias);
   }
   else if (makeBasicJet(jetTypeE)) {
     produces<reco::BasicJetCollection>(tag).setBranchAlias(alias);
   }
 }

bool VirtualJetProducer::makeTrackJet ( const JetType::Type & fTag )

inlineprotected

Definition at line 62 of file VirtualJetProducer.h.

Referenced by FastjetJetProducer::produce().

                                                     {
     return fTag == JetType::TrackJet;
   }

void VirtualJetProducer::offsetCorrectJets ( std::vector< fastjet::PseudoJet > & orphanInput )

protected

void VirtualJetProducer::output	(	edm::Event &	iEvent,
		edm::EventSetup const &	iSetup
	)

protectedvirtual

Reimplemented in cms::CompoundJetProducer, and SubjetFilterJetProducer.

Definition at line 532 of file VirtualJetProducer.cc.

References VirtualJetProducer::JetType::BasicJet, VirtualJetProducer::JetType::CaloJet, Exception, VirtualJetProducer::JetType::GenJet, iEvent, jetTypeE, VirtualJetProducer::JetType::PFClusterJet, VirtualJetProducer::JetType::PFJet, VirtualJetProducer::JetType::TrackJet, and writeCompound_.

Referenced by produce().

 {
   // Write jets and constitutents. Will use fjJets_, inputs_
   // and fjClusterSeq_
 
   if ( !writeCompound_ ) {
     switch( jetTypeE ) {
     case JetType::CaloJet :
       writeJets<reco::CaloJet>( iEvent, iSetup);
       break;
     case JetType::PFJet :
       writeJets<reco::PFJet>( iEvent, iSetup);
       break;
     case JetType::GenJet :
       writeJets<reco::GenJet>( iEvent, iSetup);
       break;
     case JetType::TrackJet :
       writeJets<reco::TrackJet>( iEvent, iSetup);
       break;
     case JetType::PFClusterJet :
       writeJets<reco::PFClusterJet>( iEvent, iSetup);
       break;
     case JetType::BasicJet :
       writeJets<reco::BasicJet>( iEvent, iSetup);
       break;
     default:
       throw cms::Exception("InvalidInput") << "invalid jet type in VirtualJetProducer\n";
       break;
     };
   } else {
     // Write jets and constitutents.
     switch( jetTypeE ) {
     case JetType::CaloJet :
       writeCompoundJets<reco::CaloJet>( iEvent, iSetup );
       break;
     case JetType::PFJet :
       writeCompoundJets<reco::PFJet>( iEvent, iSetup );
       break;
     case JetType::GenJet :
       writeCompoundJets<reco::GenJet>( iEvent, iSetup );
       break;
     case JetType::BasicJet :
       writeCompoundJets<reco::BasicJet>( iEvent, iSetup );
       break;
     default:
       throw cms::Exception("InvalidInput") << "invalid jet type in CompoundJetProducer\n";
       break;
     };
   }
 }

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

overridevirtual

Implements edm::stream::EDProducerBase.

Reimplemented in cms::CATopJetProducer, FastjetJetProducer, SubjetFilterJetProducer, CMSInsideOutJetProducer, cms::SubEventGenJetProducer, and cms::SubJetProducer.

Definition at line 300 of file VirtualJetProducer.cc.

References doPUOffsetCorr_, doPVCorrection_, edm::EventID::event(), fjInputs_, fjJetDefinition_, fjJets_, edm::Event::getByToken(), i, edm::EventBase::id(), input_candidatefwdptr_token_, input_candidateview_token_, input_packedcandidatefwdptr_token_, input_vertex_token_, inputs_, inputTowers(), jetTypeE, LogDebug, makeCaloJet(), makePFJet(), bookConverter::max, minSeed_, output(), edm::EventID::run(), runAlgorithm(), subtractor_, swap(), useDeterministicSeed_, and vertex_.

Referenced by cms::SubJetProducer::produce(), CMSInsideOutJetProducer::produce(), SubjetFilterJetProducer::produce(), and FastjetJetProducer::produce().

 {
 
   // If requested, set the fastjet random seed to a deterministic function
   // of the run/lumi/event. 
   // NOTE!!! The fastjet random number sequence is a global singleton.
   // Thus, we have to create an object and get access to the global singleton
   // in order to change it. 
   if ( useDeterministicSeed_ ) {
     fastjet::GhostedAreaSpec gas;
     std::vector<int> seeds(2);
     unsigned int runNum_uint = static_cast <unsigned int> (iEvent.id().run());
     unsigned int evNum_uint = static_cast <unsigned int> (iEvent.id().event()); 
     seeds[0] = std::max(runNum_uint,minSeed_ + 3) + 3 * evNum_uint;
     seeds[1] = std::max(runNum_uint,minSeed_ + 5) + 5 * evNum_uint;
     gas.set_random_status(seeds);
   }
 
   LogDebug("VirtualJetProducer") << "Entered produce\n";
   //determine signal vertex2
   vertex_=reco::Jet::Point(0,0,0);
   if ( (makeCaloJet(jetTypeE) || makePFJet(jetTypeE)) &&doPVCorrection_) {
     LogDebug("VirtualJetProducer") << "Adding PV info\n";
     edm::Handle<reco::VertexCollection> pvCollection;
     iEvent.getByToken(input_vertex_token_ , pvCollection);
     if (pvCollection->size()>0) vertex_=pvCollection->begin()->position();
   }
 
   // For Pileup subtraction using offset correction:
   // set up geometry map
   if ( doPUOffsetCorr_ ) {
      subtractor_->setupGeometryMap(iEvent, iSetup);
   }
 
   // clear data
   LogDebug("VirtualJetProducer") << "Clear data\n";
   fjInputs_.clear();
   fjJets_.clear();
   inputs_.clear();  
   
   // get inputs and convert them to the fastjet format (fastjet::PeudoJet)
   edm::Handle<reco::CandidateView> inputsHandle;
   
   edm::Handle< std::vector<edm::FwdPtr<reco::PFCandidate> > > pfinputsHandleAsFwdPtr; 
   edm::Handle< std::vector<edm::FwdPtr<pat::PackedCandidate> > > packedinputsHandleAsFwdPtr; 
   
   bool isView = iEvent.getByToken(input_candidateview_token_, inputsHandle);
   if ( isView ) {
     if ( inputsHandle->size() == 0) {
       output( iEvent, iSetup );
       return;
     }
     for (size_t i = 0; i < inputsHandle->size(); ++i) {
       inputs_.push_back(inputsHandle->ptrAt(i));
     }
   } else {
     bool isPF = iEvent.getByToken(input_candidatefwdptr_token_, pfinputsHandleAsFwdPtr);
     if ( isPF ) {
       if ( pfinputsHandleAsFwdPtr->size() == 0) {
     output( iEvent, iSetup );
     return;
       }
       for (size_t i = 0; i < pfinputsHandleAsFwdPtr->size(); ++i) {
     if ( (*pfinputsHandleAsFwdPtr)[i].ptr().isAvailable() ) {
       inputs_.push_back( (*pfinputsHandleAsFwdPtr)[i].ptr() );
     }
     else if ( (*pfinputsHandleAsFwdPtr)[i].backPtr().isAvailable() ) {
       inputs_.push_back( (*pfinputsHandleAsFwdPtr)[i].backPtr() );
     }
       }
     } else {
       iEvent.getByToken(input_packedcandidatefwdptr_token_, packedinputsHandleAsFwdPtr);
       if ( packedinputsHandleAsFwdPtr->size() == 0) {
     output( iEvent, iSetup );
     return;
       }
       for (size_t i = 0; i < packedinputsHandleAsFwdPtr->size(); ++i) {
     if ( (*packedinputsHandleAsFwdPtr)[i].ptr().isAvailable() ) {
       inputs_.push_back( (*packedinputsHandleAsFwdPtr)[i].ptr() );
     }
     else if ( (*packedinputsHandleAsFwdPtr)[i].backPtr().isAvailable() ) {
       inputs_.push_back( (*packedinputsHandleAsFwdPtr)[i].backPtr() );
     }
       }
     }
   }
   LogDebug("VirtualJetProducer") << "Got inputs\n";
   
   // Convert candidates to fastjet::PseudoJets.
   // Also correct to Primary Vertex. Will modify fjInputs_
   // and use inputs_
   fjInputs_.reserve(inputs_.size());
   inputTowers();
   LogDebug("VirtualJetProducer") << "Inputted towers\n";
 
   // For Pileup subtraction using offset correction:
   // Subtract pedestal. 
   if ( doPUOffsetCorr_ ) {
      subtractor_->setDefinition(fjJetDefinition_);
      subtractor_->reset(inputs_,fjInputs_,fjJets_);
      subtractor_->calculatePedestal(fjInputs_); 
      subtractor_->subtractPedestal(fjInputs_);    
      LogDebug("VirtualJetProducer") << "Subtracted pedestal\n";
   }
   // Run algorithm. Will modify fjJets_ and allocate fjClusterSeq_. 
   // This will use fjInputs_
   runAlgorithm( iEvent, iSetup );
 
   // if ( doPUOffsetCorr_ ) {
   //    subtractor_->setAlgorithm(fjClusterSeq_);
   // }
 
   LogDebug("VirtualJetProducer") << "Ran algorithm\n";
   // For Pileup subtraction using offset correction:
   // Now we find jets and need to recalculate their energy,
   // mark towers participated in jet,
   // remove occupied towers from the list and recalculate mean and sigma
   // put the initial towers collection to the jet,   
   // and subtract from initial towers in jet recalculated mean and sigma of towers 
   if ( doPUOffsetCorr_ ) {
     LogDebug("VirtualJetProducer") << "Do PUOffsetCorr\n";
     vector<fastjet::PseudoJet> orphanInput;
     subtractor_->calculateOrphanInput(orphanInput);
     subtractor_->calculatePedestal(orphanInput);
     subtractor_->offsetCorrectJets();
   }
   // Write the output jets.
   // This will (by default) call the member function template
   // "writeJets", but can be overridden. 
   // this will use inputs_
   output( iEvent, iSetup );
   LogDebug("VirtualJetProducer") << "Wrote jets\n";
   
   // Clear the work vectors so that memory is free for other modules.
   // Use the trick of swapping with an empty vector so that the memory
   // is actually given back rather than silently kept.
   decltype(fjInputs_)().swap(fjInputs_);
   decltype(fjJets_)().swap(fjJets_);
   decltype(inputs_)().swap(inputs_);  
 
   return;
 }

virtual void VirtualJetProducer::runAlgorithm	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

protectedpure virtual

Implemented in cms::CATopJetProducer, FastjetJetProducer, SubjetFilterJetProducer, CMSInsideOutJetProducer, cms::SubEventGenJetProducer, and cms::SubJetProducer.

Referenced by produce().

template<class T >

void VirtualJetProducer::writeCompoundJets	(	edm::Event &	iEvent,
		edm::EventSetup const &	iSetup
	)

protected

function template to write out the outputs

Definition at line 755 of file VirtualJetProducer.cc.

References gather_cfg::cout, doAreaFastjet_, fjJets_, getConstituents(), metsig::jet, jetCollInstanceName_, moduleLabel_, point, edm::Event::put(), reco::Jet::setJetArea(), verbosity_, and reco::writeSpecific().

 {
   if ( verbosity_ >= 1 ) { 
     std::cout << "<VirtualJetProducer::writeCompoundJets (moduleLabel = " << moduleLabel_ << ")>:" << std::endl;
   }
 
   // get a list of output jets
   std::auto_ptr<reco::BasicJetCollection>  jetCollection( new reco::BasicJetCollection() );
   // get a list of output subjets
   std::auto_ptr<std::vector<T> >  subjetCollection( new std::vector<T>() );
 
   // This will store the handle for the subjets after we write them
   edm::OrphanHandle< std::vector<T> > subjetHandleAfterPut;
   // this is the mapping of subjet to hard jet
   std::vector< std::vector<int> > indices;
   // this is the list of hardjet 4-momenta
   std::vector<math::XYZTLorentzVector> p4_hardJets;
   // this is the hardjet areas
   std::vector<double> area_hardJets;
 
   // Loop over the hard jets
   std::vector<fastjet::PseudoJet>::const_iterator it = fjJets_.begin(),
     iEnd = fjJets_.end(),
     iBegin = fjJets_.begin();
   indices.resize( fjJets_.size() );
   for ( ; it != iEnd; ++it ) {
     fastjet::PseudoJet const & localJet = *it;
     unsigned int jetIndex = it - iBegin;
     // Get the 4-vector for the hard jet
     p4_hardJets.push_back( math::XYZTLorentzVector(localJet.px(), localJet.py(), localJet.pz(), localJet.e() ));
     double localJetArea = 0.0;
     if ( doAreaFastjet_ && localJet.has_area() ) {
       localJetArea = localJet.area();
     }
     area_hardJets.push_back( localJetArea );
 
     // create the subjet list
     std::vector<fastjet::PseudoJet> constituents;
     if ( it->has_pieces() ) {
       constituents = it->pieces();
     } else if ( it->has_constituents() ) {
       constituents = it->constituents();
     }
 
     std::vector<fastjet::PseudoJet>::const_iterator itSubJetBegin = constituents.begin(),
       itSubJet = itSubJetBegin, itSubJetEnd = constituents.end();
     for (; itSubJet != itSubJetEnd; ++itSubJet ){
 
       fastjet::PseudoJet const & subjet = *itSubJet;      
       if ( verbosity_ >= 1 ) {
     std::cout << "subjet #" << (itSubJet - itSubJetBegin) << ": Pt = " << subjet.pt() << ", eta = " << subjet.eta() << ", phi = " << subjet.phi() << ", mass = " << subjet.m() 
           << " (#constituents = " << subjet.constituents().size() << ")" << std::endl;
     std::vector<fastjet::PseudoJet> subjet_constituents = subjet.constituents();
     int idx_constituent = 0;
     for ( std::vector<fastjet::PseudoJet>::const_iterator constituent = subjet_constituents.begin();
           constituent != subjet_constituents.end(); ++constituent ) {
       if ( constituent->pt() < 1.e-3 ) continue; // CV: skip ghosts
       std::cout << "  constituent #" << idx_constituent << ": Pt = " << constituent->pt() << ", eta = " << constituent->eta() << ", phi = " << constituent->phi() << "," 
             << " mass = " << constituent->m() << std::endl;
       ++idx_constituent;
     }
       }
 
       if ( verbosity_ >= 1 ) {
     std::cout << "subjet #" << (itSubJet - itSubJetBegin) << ": Pt = " << subjet.pt() << ", eta = " << subjet.eta() << ", phi = " << subjet.phi() << ", mass = " << subjet.m() 
           << " (#constituents = " << subjet.constituents().size() << ")" << std::endl;
     std::vector<fastjet::PseudoJet> subjet_constituents = subjet.constituents();
     int idx_constituent = 0;
     for ( std::vector<fastjet::PseudoJet>::const_iterator constituent = subjet_constituents.begin();
           constituent != subjet_constituents.end(); ++constituent ) {
       if ( constituent->pt() < 1.e-3 ) continue; // CV: skip ghosts
       std::cout << " constituent #" << idx_constituent << ": Pt = " << constituent->pt() << ", eta = " << constituent->eta() << ", phi = " << constituent->phi() << "," 
             << " mass = " << constituent->m() << std::endl;
       ++idx_constituent;
     }
       }
 
       math::XYZTLorentzVector p4Subjet(subjet.px(), subjet.py(), subjet.pz(), subjet.e() );
       reco::Particle::Point point(0,0,0);
 
       // This will hold ptr's to the subjets
       std::vector<reco::CandidatePtr> subjetConstituents;
 
       // Get the transient subjet constituents from fastjet
       std::vector<fastjet::PseudoJet> subjetFastjetConstituents = subjet.constituents();
       std::vector<reco::CandidatePtr> constituents =
     getConstituents(subjetFastjetConstituents );    
 
       indices[jetIndex].push_back( subjetCollection->size() );
 
       // Add the concrete subjet type to the subjet list to write to event record
       T jet;
       reco::writeSpecific( jet, p4Subjet, point, constituents, iSetup);
       double subjetArea = 0.0;
       if ( doAreaFastjet_ && itSubJet->has_area() ){
     subjetArea = itSubJet->area();
       }
       jet.setJetArea( subjetArea );
       subjetCollection->push_back( jet );
     }
   }
 
   // put subjets into event record
   subjetHandleAfterPut = iEvent.put( subjetCollection, jetCollInstanceName_ );
   
   // Now create the hard jets with ptr's to the subjets as constituents
   std::vector<math::XYZTLorentzVector>::const_iterator ip4 = p4_hardJets.begin(),
     ip4Begin = p4_hardJets.begin(),
     ip4End = p4_hardJets.end();
 
   for ( ; ip4 != ip4End; ++ip4 ) {
     int p4_index = ip4 - ip4Begin;
     std::vector<int> & ind = indices[p4_index];
     std::vector<reco::CandidatePtr> i_hardJetConstituents;
     // Add the subjets to the hard jet
     for( std::vector<int>::const_iterator isub = ind.begin();
      isub != ind.end(); ++isub ) {
       reco::CandidatePtr candPtr( subjetHandleAfterPut, *isub, false );
       i_hardJetConstituents.push_back( candPtr );
     }   
     reco::Particle::Point point(0,0,0);
     reco::BasicJet toput( *ip4, point, i_hardJetConstituents);
     toput.setJetArea( area_hardJets[ip4 - ip4Begin] );
     jetCollection->push_back( toput );
   }
 
   // put hard jets into event record
   iEvent.put( jetCollection);
 }

template<typename T >

void VirtualJetProducer::writeJets	(	edm::Event &	iEvent,
		edm::EventSetup const &	iSetup
	)

protected

Definition at line 595 of file VirtualJetProducer.cc.

References reco::deltaR2(), HLT_25ns10e33_v2_cff::distance, doAreaDiskApprox_, doAreaFastjet_, doFastJetNonUniform_, doPUOffsetCorr_, doRhoFastjet_, eta, Exception, plotBeamSpotDB::first, fjClusterSeq_, fjJets_, fjRangeDef_, getConstituents(), reco::helper::VirtualJetProducerHelper::intersection(), edm::isNotFinite(), metsig::jet, jetCollInstanceName_, fwrapper::jets, relval_2017::k, M_PI, eostools::move(), HLT_25ns10e33_v2_cff::nEta, nExclude_, puCenters_, edm::Event::put(), puWidth_, dt_dqm_sourceclient_common_cff::reco, rho, fastjet::BackgroundEstimator::rho(), rParam_, edm::second(), fastjet::BackgroundEstimator::set_excluded_jets(), fastjet::BackgroundEstimator::sigma(), mathSSE::sqrt(), subtractor_, vertex_, and reco::writeSpecific().

 {
   // std::cout << "writeJets " << typeid(T).name() 
   //          << (doRhoFastjet_ ? " doRhoFastjet " : "")
   //          << (doAreaFastjet_ ? " doAreaFastjet " : "")
   //          << (doAreaDiskApprox_ ? " doAreaDiskApprox " : "")
   //          << std::endl;
 
   if (doRhoFastjet_) {
     // declare jet collection without the two jets, 
     // for unbiased background estimation.
     std::vector<fastjet::PseudoJet> fjexcluded_jets;
     fjexcluded_jets=fjJets_;
     
     if(fjexcluded_jets.size()>2) fjexcluded_jets.resize(nExclude_);
     
     if(doFastJetNonUniform_){
       std::auto_ptr<std::vector<double> > rhos(new std::vector<double>);
       std::auto_ptr<std::vector<double> > sigmas(new std::vector<double>);
       int nEta = puCenters_.size();
       rhos->reserve(nEta);
       sigmas->reserve(nEta);
       fastjet::ClusterSequenceAreaBase const* clusterSequenceWithArea =
         dynamic_cast<fastjet::ClusterSequenceAreaBase const *> ( &*fjClusterSeq_ );
 
       if (clusterSequenceWithArea ==nullptr ){
     if (fjJets_.size() > 0) {
       throw cms::Exception("LogicError")<<"fjClusterSeq is not initialized while inputs are present\n ";
     }
       } else {
     for(int ie = 0; ie < nEta; ++ie){
       double eta = puCenters_[ie];
       double etamin=eta-puWidth_;
       double etamax=eta+puWidth_;
       fastjet::RangeDefinition range_rho(etamin,etamax);
       fastjet::BackgroundEstimator bkgestim(*clusterSequenceWithArea,range_rho);
       bkgestim.set_excluded_jets(fjexcluded_jets);
       rhos->push_back(bkgestim.rho());
       sigmas->push_back(bkgestim.sigma());
     }
       }
       iEvent.put(rhos,"rhos");
       iEvent.put(sigmas,"sigmas");
     }else{
       std::auto_ptr<double> rho(new double(0.0));
       std::auto_ptr<double> sigma(new double(0.0));
       double mean_area = 0;
       
       fastjet::ClusterSequenceAreaBase const* clusterSequenceWithArea =
         dynamic_cast<fastjet::ClusterSequenceAreaBase const *> ( &*fjClusterSeq_ );
       /*
     const double nemptyjets = clusterSequenceWithArea->n_empty_jets(*fjRangeDef_);
     if(( nemptyjets  < -15 ) || ( nemptyjets > fjRangeDef_->area()+ 15)) {
     edm::LogWarning("StrangeNEmtpyJets") << "n_empty_jets is : " << clusterSequenceWithArea->n_empty_jets(*fjRangeDef_) << " with range " << fjRangeDef_->description() << ".";
     }
       */
       if (clusterSequenceWithArea ==nullptr ){
     if (fjJets_.size() > 0) {
       throw cms::Exception("LogicError")<<"fjClusterSeq is not initialized while inputs are present\n ";
     }
       } else {
     clusterSequenceWithArea->get_median_rho_and_sigma(*fjRangeDef_,false,*rho,*sigma,mean_area);
     if((*rho < 0)|| (edm::isNotFinite(*rho))) {
       edm::LogError("BadRho") << "rho value is " << *rho << " area:" << mean_area << " and n_empty_jets: " << clusterSequenceWithArea->n_empty_jets(*fjRangeDef_) << " with range " << fjRangeDef_->description()
                   <<". Setting rho to rezo.";
       *rho = 0;
     }
       }
       iEvent.put(rho,"rho");
       iEvent.put(sigma,"sigma");
     }
   } // doRhoFastjet_
   
   // produce output jet collection
   
   using namespace reco;
 
   // allocate fjJets_.size() Ts in vector  
   auto jets = std::make_unique<std::vector<T>>(fjJets_.size());
   
   // Distance between jet centers and overlap area -- for disk-based area calculation
   using RIJ = std::pair<double,double>; 
   std::vector<RIJ>   rijStorage(fjJets_.size()*(fjJets_.size()/2));
   RIJ * rij[fjJets_.size()];
   unsigned int k=0;
   for (unsigned int ijet=0;ijet<fjJets_.size();++ijet) {
      rij[ijet] = &rijStorage[k]; k+=ijet;
   }
 
   float etaJ[fjJets_.size()],  phiJ[fjJets_.size()];
 
   auto orParam_ = 1./rParam_;
   // fill jets 
   for (unsigned int ijet=0;ijet<fjJets_.size();++ijet) {
     auto & jet = (*jets)[ijet];
     // get the fastjet jet
     const fastjet::PseudoJet& fjJet = fjJets_[ijet];
     // get the constituents from fastjet
     std::vector<fastjet::PseudoJet> const & fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
     // convert them to CandidatePtr vector
     std::vector<CandidatePtr> const & constituents = getConstituents(fjConstituents);
 
     // write the specifics to the jet (simultaneously sets 4-vector, vertex).
     // These are overridden functions that will call the appropriate
     // specific allocator.
     writeSpecific(jet,
                   Particle::LorentzVector(fjJet.px(),
                                           fjJet.py(),
                                           fjJet.pz(),
                                           fjJet.E()),
                   vertex_,
                   constituents, iSetup);
     phiJ[ijet] = jet.phi();
     etaJ[ijet] = jet.eta();
   }
 
    // calcuate the jet area
   for (unsigned int ijet=0;ijet<fjJets_.size();++ijet) {
     // calcuate the jet area
     double jetArea=0.0;
     // get the fastjet jet
     const auto & fjJet = fjJets_[ijet];
     if ( doAreaFastjet_ && fjJet.has_area() ) {
       jetArea = fjJet.area();
     }
     else if ( doAreaDiskApprox_ ) {
       // Here it is assumed that fjJets_ is in decreasing order of pT, 
       // which should happen in FastjetJetProducer::runAlgorithm() 
       jetArea   = M_PI;
         RIJ *  distance  = rij[ijet];
         for (unsigned jJet = 0; jJet < ijet; ++jJet) {
           distance[jJet].first      = std::sqrt(reco::deltaR2(etaJ[ijet],phiJ[ijet], etaJ[jJet],phiJ[jJet]))*orParam_;
           distance[jJet].second = reco::helper::VirtualJetProducerHelper::intersection(distance[jJet].first);
           jetArea            -=distance[jJet].second;
           for (unsigned kJet = 0; kJet < jJet; ++kJet) {
             jetArea          += reco::helper::VirtualJetProducerHelper::intersection(distance[jJet].first, distance[kJet].first, rij[jJet][kJet].first, 
                                                                                      distance[jJet].second, distance[kJet].second, rij[jJet][kJet].second);
           } // end loop over harder jets
         } // end loop over harder jets
       jetArea  *= (rParam_*rParam_);
     } 
     auto & jet = (*jets)[ijet]; 
     jet.setJetArea (jetArea);
     
     if(doPUOffsetCorr_){
       jet.setPileup(subtractor_->getPileUpEnergy(ijet));
     }else{
       jet.setPileup (0.0);
     }
         
    // std::cout << "area " << ijet << " " << jetArea << " " << Area<T>::get(jet) << std::endl;
    // std::cout << "JetVI " << ijet << ' '<< jet.pt() << " " << jet.et() << ' '<< jet.energy() << ' '<< jet.mass() << std::endl;
 
   }
   // put the jets in the collection
   iEvent.put(std::move(jets),jetCollInstanceName_);
 }