#include <MultipleAlgoIterator.h>

Inheritance diagram for MultipleAlgoIterator:

Public Member Functions
void	calculateOrphanInput (std::vector< fastjet::PseudoJet > &orphanInput) override

void	calculatePedestal (std::vector< fastjet::PseudoJet > const &coll) override

double	getEt (const reco::CandidatePtr &in) const

double	getEta (const reco::CandidatePtr &in) const

	MultipleAlgoIterator (const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iC)

void	offsetCorrectJets () override

void	rescaleRMS (double s)

void	subtractPedestal (std::vector< fastjet::PseudoJet > &coll) override

Public Member Functions inherited from PileUpSubtractor
virtual double	getCone (double cone, double eta, double phi, double &et, double &pu)

virtual double	getMeanAtTower (const reco::CandidatePtr &in) const

int	getN (const reco::CandidatePtr &in) const

int	getNwithJets (const reco::CandidatePtr &in) const

virtual double	getPileUpAtTower (const reco::CandidatePtr &in) const

virtual double	getPileUpEnergy (int ijet) const

virtual double	getSigmaAtTower (const reco::CandidatePtr &in) const

int	ieta (const reco::CandidatePtr &in) const

int	iphi (const reco::CandidatePtr &in) const

	PileUpSubtractor (const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iC)

virtual void	reset (std::vector< edm::Ptr< reco::Candidate > > &input, std::vector< fastjet::PseudoJet > &towers, std::vector< fastjet::PseudoJet > &output)

virtual void	setDefinition (JetDefPtr const &jetDef)

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

virtual	~PileUpSubtractor ()

Private Attributes
bool	dropZeroTowers_

double	minimumTowersFraction_

bool	sumRecHits_

Additional Inherited Members
Public Types inherited from PileUpSubtractor
typedef std::shared_ptr < fastjet::GhostedAreaSpec >	ActiveAreaSpecPtr

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

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

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

Protected Attributes inherited from PileUpSubtractor
int	activeAreaRepeats

std::vector< HcalDetId >	allgeomid_

bool	doAreaFastjet_

bool	doRhoFastjet_

std::map< int, double >	emean_

std::map< int, double >	esigma_

ActiveAreaSpecPtr	fjActiveArea_

ClusterSequencePtr	fjClusterSeq_

std::vector< fastjet::PseudoJet > *	fjInputs_

JetDefPtr	fjJetDefinition_

std::vector< fastjet::PseudoJet > *	fjJets_

std::vector< fastjet::PseudoJet >	fjOriginalInputs_

CaloGeometry const *	geo_

std::map< int, int >	geomtowers_

edm::ESGetToken< CaloGeometry, CaloGeometryRecord >	geoToken_

double	ghostArea

double	ghostEtaMax

int	ietamax_

int	ietamin_

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

std::vector< double >	jetOffset_

double	jetPtMin_

double	nSigmaPU_

std::map< int, int >	ntowersWithJets_

double	puPtMin_

double	radiusPU_

bool	reRunAlgo_

Detailed Description

Definition at line 6 of file MultipleAlgoIterator.h.

Constructor & Destructor Documentation

MultipleAlgoIterator::MultipleAlgoIterator	(	const edm::ParameterSet &	iConfig,
		edm::ConsumesCollector &&	iC
	)

Definition at line 11 of file MultipleAlgoIterator.cc.

References LogDebug, and minimumTowersFraction_.

     : PileUpSubtractor(iConfig, std::move(iC)),
       minimumTowersFraction_(iConfig.getParameter<double>("minimumTowersFraction")),
       sumRecHits_(iConfig.getParameter<bool>("sumRecHits")),
       dropZeroTowers_(iConfig.getUntrackedParameter<bool>("dropZeroTowers", true)) {
   LogDebug("PileUpSubtractor") << "LIMITING THE MINIMUM TOWERS FRACTION TO : " << minimumTowersFraction_ << endl;
 }

Member Function Documentation

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

overridevirtual

Reimplemented from PileUpSubtractor.

Definition at line 172 of file MultipleAlgoIterator.cc.

References PileUpSubtractor::allgeomid_, reco::deltaR(), runTauDisplay::dr, spr::find(), PileUpSubtractor::fjInputs_, PileUpSubtractor::fjJets_, PileUpSubtractor::fjOriginalInputs_, PileUpSubtractor::geo_, PileUpSubtractor::geomtowers_, CaloGeometry::getPosition(), PileUpSubtractor::ieta(), PileUpSubtractor::inputs_, PileUpSubtractor::iphi(), LogDebug, minimumTowersFraction_, PileUpSubtractor::ntowersWithJets_, PileUpSubtractor::puPtMin_, and PileUpSubtractor::radiusPU_.

                                                                                      {
   LogDebug("PileUpSubtractor") << "The subtractor calculating orphan input...\n";
 
   (*fjInputs_) = fjOriginalInputs_;
 
   vector<int> jettowers;                   // vector of towers indexed by "user_index"
   vector<pair<int, int> > excludedTowers;  // vector of excluded ieta, iphi values
 
   for (auto const& pseudojetTMP : *fjJets_) {
     if (pseudojetTMP.perp() < puPtMin_)
       continue;
 
     // find towers within radiusPU_ of this jet
     for (auto const im : allgeomid_) {
       double dr = reco::deltaR(geo_->getPosition(im), pseudojetTMP);
       auto exclude = find(excludedTowers.begin(), excludedTowers.end(), pair<int, int>(im.ieta(), im.iphi()));
       if (dr < radiusPU_ && exclude == excludedTowers.end() &&
           (geomtowers_[im.ieta()] - ntowersWithJets_[im.ieta()]) > minimumTowersFraction_ * (geomtowers_[im.ieta()])) {
         ntowersWithJets_[im.ieta()]++;
 
         excludedTowers.push_back(pair<int, int>(im.ieta(), im.iphi()));
       }
     }
 
     for (auto const& it : *fjInputs_) {
       int index = it.user_index();
       int ie = ieta((*inputs_)[index]);
       int ip = iphi((*inputs_)[index]);
       auto exclude = find(excludedTowers.begin(), excludedTowers.end(), pair<int, int>(ie, ip));
       if (exclude != excludedTowers.end()) {
         jettowers.push_back(index);
       }
     }  // initial input collection
 
   }  // pseudojets
 
   //
   // Create a new collections from the towers not included in jets
   //
 
   for (auto const& it : *fjInputs_) {
     int index = it.user_index();
     vector<int>::const_iterator itjet = find(jettowers.begin(), jettowers.end(), index);
     if (itjet == jettowers.end()) {
       const reco::CandidatePtr& originalTower = (*inputs_)[index];
       fastjet::PseudoJet orphan(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy());
       orphan.set_user_index(index);
 
       orphanInput.push_back(orphan);
     }
   }
 }

void MultipleAlgoIterator::calculatePedestal ( std::vector< fastjet::PseudoJet > const & coll )

overridevirtual

Reimplemented from PileUpSubtractor.

Definition at line 108 of file MultipleAlgoIterator.cc.

References PileUpSubtractor::emean_, PileUpSubtractor::esigma_, PileUpSubtractor::geomtowers_, getEt(), submitPVValidationJobs::gt, mps_fire::i, PileUpSubtractor::ieta(), PileUpSubtractor::ietamax_, PileUpSubtractor::ietamin_, LogDebug, PileUpSubtractor::nSigmaPU_, nt, PileUpSubtractor::ntowersWithJets_, edm::second(), mathSSE::sqrt(), and sumRecHits_.

                                                                                  {
   LogDebug("PileUpSubtractor") << "The subtractor calculating pedestals...\n";
 
   map<int, double> emean2;
   map<int, int> ntowers;
 
   int ietaold = -10000;
   int ieta0 = -100;
 
   // Initial values for emean_, emean2, esigma_, ntowers
 
   for (int i = ietamin_; i < ietamax_ + 1; i++) {
     emean_[i] = 0.;
     emean2[i] = 0.;
     esigma_[i] = 0.;
     ntowers[i] = 0;
   }
 
   for (vector<fastjet::PseudoJet>::const_iterator input_object = coll.begin(), fjInputsEnd = coll.end();
        input_object != fjInputsEnd;
        ++input_object) {
     const reco::CandidatePtr& originalTower = (*inputs_)[input_object->user_index()];
     ieta0 = ieta(originalTower);
     double Original_Et = originalTower->et();
     if (sumRecHits_) {
       Original_Et = getEt(originalTower);
     }
 
     if (ieta0 - ietaold != 0) {
       emean_[ieta0] = emean_[ieta0] + Original_Et;
       emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et;
       ntowers[ieta0] = 1;
       ietaold = ieta0;
     } else {
       emean_[ieta0] = emean_[ieta0] + Original_Et;
       emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et;
       ntowers[ieta0]++;
     }
   }
 
   for (map<int, int>::const_iterator gt = geomtowers_.begin(); gt != geomtowers_.end(); gt++) {
     int it = (*gt).first;
 
     double e1 = (*(emean_.find(it))).second;
     double e2 = (*emean2.find(it)).second;
     int nt = (*gt).second - (*(ntowersWithJets_.find(it))).second;
 
     LogDebug("PileUpSubtractor") << " ieta : " << it << " number of towers : " << nt << " e1 : " << e1 << " e2 : " << e2
                                  << "\n";
 
     if (nt > 0) {
       emean_[it] = e1 / nt;
       double eee = e2 / nt - e1 * e1 / (nt * nt);
       if (eee < 0.)
         eee = 0.;
       esigma_[it] = nSigmaPU_ * sqrt(eee);
     } else {
       emean_[it] = 0.;
       esigma_[it] = 0.;
     }
     LogDebug("PileUpSubtractor") << " ieta : " << it << " Pedestals : " << emean_[it] << "  " << esigma_[it] << "\n";
   }
 }

double MultipleAlgoIterator::getEt ( const reco::CandidatePtr & in ) const

Definition at line 225 of file MultipleAlgoIterator.cc.

References CaloTower::emEnergy(), relval_parameters_module::energy, PileUpSubtractor::geo_, edm::Ptr< T >::get(), CaloGeometry::getPosition(), CaloTower::hadEnergy(), CaloTower::id(), funct::sin(), and PV3DBase< T, PVType, FrameType >::theta().

Referenced by calculatePedestal(), and subtractPedestal().

                                                                    {
   const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
   const GlobalPoint& pos = geo_->getPosition(ctc->id());
   double energy = ctc->emEnergy() + ctc->hadEnergy();
   double et = energy * sin(pos.theta());
   return et;
 }

double MultipleAlgoIterator::getEta ( const reco::CandidatePtr & in ) const

Definition at line 233 of file MultipleAlgoIterator.cc.

References PVValHelper::eta, PV3DBase< T, PVType, FrameType >::eta(), PileUpSubtractor::geo_, edm::Ptr< T >::get(), CaloGeometry::getPosition(), and CaloTower::id().

                                                                     {
   const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
   const GlobalPoint& pos = geo_->getPosition(ctc->id());
   double eta = pos.eta();
   return eta;
 }

void MultipleAlgoIterator::offsetCorrectJets ( )

overridevirtual

Reimplemented from PileUpSubtractor.

Definition at line 25 of file MultipleAlgoIterator.cc.

                                              {
   LogDebug("PileUpSubtractor") << "The subtractor correcting jets...\n";
   jetOffset_.clear();
 
   using namespace reco;
 
   (*fjInputs_) = fjOriginalInputs_;
   rescaleRMS(nSigmaPU_);
   subtractPedestal(*fjInputs_);
   const fastjet::JetDefinition& def = *fjJetDefinition_;
   if (!doAreaFastjet_ && !doRhoFastjet_) {
     fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(*fjInputs_, def);
   } else {
     fjClusterSeq_ = ClusterSequencePtr(new fastjet::ClusterSequenceArea(*fjInputs_, def, *fjActiveArea_));
   }
 
   (*fjJets_) = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
 
   jetOffset_.reserve(fjJets_->size());
 
   vector<fastjet::PseudoJet>::iterator pseudojetTMP = fjJets_->begin(), jetsEnd = fjJets_->end();
   for (; pseudojetTMP != jetsEnd; ++pseudojetTMP) {
     int ijet = pseudojetTMP - fjJets_->begin();
     jetOffset_[ijet] = 0;
 
     std::vector<fastjet::PseudoJet> towers = sorted_by_pt(pseudojetTMP->constituents());
 
     double newjetet = 0.;
     for (vector<fastjet::PseudoJet>::const_iterator ito = towers.begin(), towEnd = towers.end(); ito != towEnd; ++ito) {
       const reco::CandidatePtr& originalTower = (*inputs_)[ito->user_index()];
       int it = ieta(originalTower);
       double Original_Et = originalTower->et();
       double etnew = Original_Et - (*emean_.find(it)).second - (*esigma_.find(it)).second;
       if (etnew < 0.)
         etnew = 0;
       newjetet = newjetet + etnew;
       jetOffset_[ijet] += Original_Et - etnew;
     }
   }
 }

void MultipleAlgoIterator::rescaleRMS ( double s )

Definition at line 19 of file MultipleAlgoIterator.cc.

References PileUpSubtractor::esigma_.

Referenced by offsetCorrectJets().

                                               {
   for (std::map<int, double>::iterator iter = esigma_.begin(); iter != esigma_.end(); ++iter) {
     iter->second = s * (iter->second);
   }
 }

void MultipleAlgoIterator::subtractPedestal ( std::vector< fastjet::PseudoJet > & coll )

overridevirtual

Reimplemented from PileUpSubtractor.

Definition at line 66 of file MultipleAlgoIterator.cc.

References dropZeroTowers_, PileUpSubtractor::emean_, PileUpSubtractor::esigma_, getEt(), PileUpSubtractor::ieta(), PileUpSubtractor::iphi(), LogDebug, and sumRecHits_.

Referenced by offsetCorrectJets().

                                                                           {
   LogDebug("PileUpSubtractor") << "The subtractor subtracting pedestals...\n";
 
   int it = -100;
 
   vector<fastjet::PseudoJet> newcoll;
 
   for (vector<fastjet::PseudoJet>::iterator input_object = coll.begin(), fjInputsEnd = coll.end();
        input_object != fjInputsEnd;
        ++input_object) {
     reco::CandidatePtr const& itow = (*inputs_)[input_object->user_index()];
 
     it = ieta(itow);
     iphi(itow);
 
     double Original_Et = itow->et();
     if (sumRecHits_) {
       Original_Et = getEt(itow);
     }
 
     double etnew = Original_Et - (*(emean_.find(it))).second - (*(esigma_.find(it))).second;
     float mScale = etnew / input_object->Et();
     if (etnew < 0.)
       mScale = 0.;
 
     math::XYZTLorentzVectorD towP4(input_object->px() * mScale,
                                    input_object->py() * mScale,
                                    input_object->pz() * mScale,
                                    input_object->e() * mScale);
 
     int index = input_object->user_index();
     input_object->reset_momentum(towP4.px(), towP4.py(), towP4.pz(), towP4.energy());
     input_object->set_user_index(index);
 
     if (etnew > 0. && dropZeroTowers_)
       newcoll.push_back(*input_object);
   }
 
   if (dropZeroTowers_)
     coll = newcoll;
 }