#include <PileUpSubtractor.h>

Inheritance diagram for PileUpSubtractor:

Public Types
typedef boost::shared_ptr < fastjet::GhostedAreaSpec >	ActiveAreaSpecPtr
typedef boost::shared_ptr < fastjet::ClusterSequence >	ClusterSequencePtr
typedef boost::shared_ptr < fastjet::JetDefinition >	JetDefPtr
typedef boost::shared_ptr < fastjet::RangeDefinition >	RangeDefPtr
Public Member Functions
virtual void	calculateOrphanInput (std::vector< fastjet::PseudoJet > &orphanInput)
virtual void	calculatePedestal (std::vector< fastjet::PseudoJet > const &coll)
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
virtual void	offsetCorrectJets ()
	PileUpSubtractor (const edm::ParameterSet &iConfig)
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 void	subtractPedestal (std::vector< fastjet::PseudoJet > &coll)
virtual	~PileUpSubtractor ()
Protected Attributes
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_
RangeDefPtr	fjRangeDef_
CaloGeometry const *	geo_
std::map< int, int >	geomtowers_
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 19 of file PileUpSubtractor.h.

Member Typedef Documentation

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

Definition at line 24 of file PileUpSubtractor.h.

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

Definition at line 23 of file PileUpSubtractor.h.

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

Definition at line 26 of file PileUpSubtractor.h.

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

Definition at line 25 of file PileUpSubtractor.h.

Constructor & Destructor Documentation

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

Definition at line 18 of file PileUpSubtractor.cc.

References doAreaFastjet_, doRhoFastjet_, edm::ParameterSet::exists(), fjActiveArea_, fjRangeDef_, edm::ParameterSet::getParameter(), and puPtMin_.

                                                                 :
  reRunAlgo_ (iConfig.getUntrackedParameter<bool>("reRunAlgo",false)),
  doAreaFastjet_ (iConfig.getParameter<bool>         ("doAreaFastjet")),
  doRhoFastjet_  (iConfig.getParameter<bool>         ("doRhoFastjet")),
  jetPtMin_(iConfig.getParameter<double>       ("jetPtMin")),
  nSigmaPU_(iConfig.getParameter<double>("nSigmaPU")),
  radiusPU_(iConfig.getParameter<double>("radiusPU")),
  geo_(0)
{
  if (iConfig.exists("puPtMin"))
    puPtMin_=iConfig.getParameter<double>       ("puPtMin");
  else{
    puPtMin_=10;
    edm::LogWarning("MisConfiguration")<<"the parameter puPtMin is now necessary for PU substraction. setting it to "<<puPtMin_;
  }
   if ( doAreaFastjet_ || doRhoFastjet_ ) {
      // 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");
      fjActiveArea_ =  ActiveAreaSpecPtr(new fastjet::ActiveAreaSpec(ghostEtaMax,
                                                                     activeAreaRepeats,
                                                                     ghostArea));
      fjRangeDef_ = RangeDefPtr( new fastjet::RangeDefinition(ghostEtaMax) );
   } 
}

virtual PileUpSubtractor::~PileUpSubtractor ( ) [inline, virtual]

Definition at line 29 of file PileUpSubtractor.h.

{;}

Member Function Documentation

void PileUpSubtractor::calculateOrphanInput ( std::vector< fastjet::PseudoJet > & orphanInput ) [virtual]

Reimplemented in ParametrizedSubtractor.

Definition at line 147 of file ParametrizedSubtractor.cc.

References fjInputs_.

{
   orphanInput = *fjInputs_;
}

void PileUpSubtractor::calculatePedestal ( std::vector< fastjet::PseudoJet > const & coll ) [virtual]

Reimplemented in MultipleAlgoIterator, ParametrizedSubtractor, and ReflectedIterator.

Definition at line 109 of file MultipleAlgoIterator.cc.

References gt, i, LogDebug, nt, edm::second(), and mathSSE::sqrt().

{
   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 PileUpSubtractor::getCone	(	double	cone,
		double	eta,
		double	phi,
		double &	et,
		double &	pu
	)		`[virtual]`

Definition at line 308 of file PileUpSubtractor.cc.

References allgeomid_, deltaR(), emean_, esigma_, eta(), PV3DBase< T, PVType, FrameType >::eta(), geo_, CaloGeometry::getPosition(), PV3DBase< T, PVType, FrameType >::phi(), phi, and point.

                                                                                           {
  pu = 0;
  
  for(vector<HcalDetId>::const_iterator im = allgeomid_.begin(); im != allgeomid_.end(); im++){
     if( im->depth() != 1 ) continue;    
    const GlobalPoint& point = geo_->getPosition((DetId)(*im));
    double dr = reco::deltaR(point.eta(),point.phi(),eta,phi);      
    if( dr < cone){
      pu += (*emean_.find(im->ieta())).second+(*esigma_.find(im->ieta())).second;
    }
  }
  
  return pu;
}

double PileUpSubtractor::getMeanAtTower ( const reco::CandidatePtr & in ) const [virtual]

Reimplemented in ParametrizedSubtractor.

Definition at line 323 of file PileUpSubtractor.cc.

References emean_, ieta(), and edm::second().

                                                                        {
  int it = ieta(in);
  return (*emean_.find(it)).second;
}

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

Definition at line 338 of file PileUpSubtractor.cc.

References geomtowers_, ieta(), n, and ntowersWithJets_.

                                                            {
   int it = ieta(in);
   
   int n = (*(geomtowers_.find(it))).second - (*(ntowersWithJets_.find(it))).second;
   return n;

}

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

Definition at line 346 of file PileUpSubtractor.cc.

References ieta(), n, and ntowersWithJets_.

                                                                    {
   int it = ieta(in);
   int n = (*(ntowersWithJets_.find(it))).second;
   return n;

}

double PileUpSubtractor::getPileUpAtTower ( const reco::CandidatePtr & in ) const [virtual]

Reimplemented in ParametrizedSubtractor.

Definition at line 333 of file PileUpSubtractor.cc.

References emean_, esigma_, ieta(), and edm::second().

                                                                           {
  int it = ieta(in);
  return (*emean_.find(it)).second + (*esigma_.find(it)).second;
}

virtual double PileUpSubtractor::getPileUpEnergy ( int ijet ) const [inline, virtual]

Definition at line 43 of file PileUpSubtractor.h.

References jetOffset_.

{return jetOffset_[ijet];}

double PileUpSubtractor::getSigmaAtTower ( const reco::CandidatePtr & in ) const [virtual]

Reimplemented in ParametrizedSubtractor.

Definition at line 328 of file PileUpSubtractor.cc.

References esigma_, ieta(), and edm::second().

                                                                          {
   int it = ieta(in);
   return (*esigma_.find(it)).second;
}

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

Definition at line 354 of file PileUpSubtractor.cc.

References Exception, edm::Ptr< T >::get(), CaloTower::id(), and CaloTowerDetId::ieta().

Referenced by ParametrizedSubtractor::getMeanAtTower(), getMeanAtTower(), getN(), getNwithJets(), getPileUpAtTower(), ParametrizedSubtractor::getPileUpAtTower(), ParametrizedSubtractor::getSigmaAtTower(), getSigmaAtTower(), offsetCorrectJets(), ParametrizedSubtractor::offsetCorrectJets(), setupGeometryMap(), and ParametrizedSubtractor::setupGeometryMap().

                                                            {
  int it = 0;
  const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
  if(ctc){
    it = ctc->id().ieta();
  } else
    {
      throw cms::Exception("Invalid Constituent") << "CaloJet constituent is not of CaloTower type";
    }
  return it;
}

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

Definition at line 366 of file PileUpSubtractor.cc.

References Exception, edm::Ptr< T >::get(), CaloTower::id(), and CaloTowerDetId::iphi().

                                                            {
  int it = 0;
  const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
  if(ctc){
    it = ctc->id().iphi();
  } else
    {
      throw cms::Exception("Invalid Constituent") << "CaloJet constituent is not of CaloTower type";
    }
  return it;
}

void PileUpSubtractor::offsetCorrectJets ( ) [virtual]

Reimplemented in MultipleAlgoIterator, ParametrizedSubtractor, and ReflectedIterator.

Definition at line 262 of file PileUpSubtractor.cc.

References emean_, esigma_, fjJets_, ieta(), jetOffset_, LogDebug, L1Trigger_dataformats::reco, and edm::second().

{
  LogDebug("PileUpSubtractor")<<"The subtractor correcting jets...\n";
  jetOffset_.clear();
  using namespace reco;
  
  //    
  // Reestimate energy of jet (energy of jet with initial map)
  //
  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 =
      fastjet::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;
      }
    double mScale = newjetet/pseudojetTMP->Et();
    LogDebug("PileUpSubtractor")<<"pseudojetTMP->Et() : "<<pseudojetTMP->Et()<<"\n";
    LogDebug("PileUpSubtractor")<<"newjetet : "<<newjetet<<"\n";
    LogDebug("PileUpSubtractor")<<"jetOffset_[ijet] : "<<jetOffset_[ijet]<<"\n";
    LogDebug("PileUpSubtractor")<<"pseudojetTMP->Et() - jetOffset_[ijet] : "<<pseudojetTMP->Et() - jetOffset_[ijet]<<"\n";
    LogDebug("PileUpSubtractor")<<"Scale is : "<<mScale<<"\n";
    int cshist = pseudojetTMP->cluster_hist_index();
    pseudojetTMP->reset_momentum(pseudojetTMP->px()*mScale, pseudojetTMP->py()*mScale,
                                 pseudojetTMP->pz()*mScale, pseudojetTMP->e()*mScale);
    pseudojetTMP->set_cluster_hist_index(cshist);
    
  }
}

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

Definition at line 47 of file PileUpSubtractor.cc.

References fjInputs_, fjJets_, fjOriginalInputs_, i, LaserDQM_cfg::input, inputs_, and convertSQLitetoXML_cfg::output.

                                                                 {
  
  inputs_ = &input;
  fjInputs_ = &towers;
  fjJets_ = &output;
  fjOriginalInputs_ = (*fjInputs_);
  for (unsigned int i = 0; i < fjInputs_->size(); ++i){
    fjOriginalInputs_[i].set_user_index((*fjInputs_)[i].user_index());
  }
  
}

void PileUpSubtractor::setDefinition ( JetDefPtr const & jetDef ) [virtual]

Definition at line 61 of file PileUpSubtractor.cc.

References fjJetDefinition_.

                                                            {
  fjJetDefinition_ = JetDefPtr( new fastjet::JetDefinition( *jetDef ) );
}

void PileUpSubtractor::setupGeometryMap	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)		`[virtual]`

Reimplemented in ParametrizedSubtractor.

Definition at line 65 of file PileUpSubtractor.cc.

References allgeomid_, geo_, geomtowers_, edm::EventSetup::get(), CaloGeometry::getValidDetIds(), DetId::Hcal, i, ieta(), ietamax_, ietamin_, LogDebug, ntowersWithJets_, and edm::ESHandle< T >::product().

{

  LogDebug("PileUpSubtractor")<<"The subtractor setting up geometry...\n";

  if(geo_ == 0) {
    edm::ESHandle<CaloGeometry> pG;
    iSetup.get<CaloGeometryRecord>().get(pG);
    geo_ = pG.product();
    std::vector<DetId> alldid =  geo_->getValidDetIds();

    int ietaold = -10000;
    ietamax_ = -10000;
    ietamin_ = 10000;
    for(std::vector<DetId>::const_iterator did=alldid.begin(); did != alldid.end(); did++){
      if( (*did).det() == DetId::Hcal ){
        HcalDetId hid = HcalDetId(*did);
        if( (hid).depth() == 1 ) {
          allgeomid_.push_back(*did);

          if((hid).ieta() != ietaold){
            ietaold = (hid).ieta();
            geomtowers_[(hid).ieta()] = 1;
            if((hid).ieta() > ietamax_) ietamax_ = (hid).ieta();
            if((hid).ieta() < ietamin_) ietamin_ = (hid).ieta();
          }
          else{
            geomtowers_[(hid).ieta()]++;
          }
        }
      }
    }
  }

  for (int i = ietamin_; i<ietamax_+1; i++) {
    ntowersWithJets_[i] = 0;
  }
}

void PileUpSubtractor::subtractPedestal ( std::vector< fastjet::PseudoJet > & coll ) [virtual]

Reimplemented in MultipleAlgoIterator, ParametrizedSubtractor, and ReflectedIterator.

Definition at line 65 of file MultipleAlgoIterator.cc.

References getHLTprescales::index, and LogDebug.

{

   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;

}