#include <ParametrizedSubtractor.h>

Inheritance diagram for ParametrizedSubtractor:

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

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

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

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

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

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

double	getPU (int ieta, bool addMean, bool addSigma) const

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

virtual void	offsetCorrectJets ()

	ParametrizedSubtractor (const edm::ParameterSet &iConfig)

void	rescaleRMS (double s)

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

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

	~ParametrizedSubtractor ()

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

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

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

virtual double	getPileUpEnergy (int ijet) const

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

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

	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)

	~PileUpSubtractor ()

Public Attributes
int	bin_

const CentralityBins *	cbins_

double	centrality_

edm::InputTag	centTag_

bool	dropZeroTowers_

TF1 *	fMean

TF1 *	fPU

TF1 *	fRMS

TH1D *	hC

std::vector< TH1D * >	hEta

std::vector< TH1D * >	hEtaMean

std::vector< TH1D * >	hEtaRMS

bool	interpolate_

bool	sumRecHits_

Additional Inherited Members
Public Types inherited from PileUpSubtractor
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

Protected Attributes inherited from PileUpSubtractor
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 13 of file ParametrizedSubtractor.h.

Constructor & Destructor Documentation

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

Definition at line 22 of file ParametrizedSubtractor.cc.

References centTag_, data, fMean, fPU, fRMS, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hC, hEta, hEtaMean, hEtaRMS, i, EcalCondDB::inf, interpolate_, and sumRecHits_.

                                                                              : 
    PileUpSubtractor(iConfig),
    dropZeroTowers_(iConfig.getUntrackedParameter<bool>("dropZeroTowers",true)),
    cbins_(0)
 {
    centTag_ = iConfig.getUntrackedParameter<edm::InputTag>("centTag",edm::InputTag("hiCentrality","","RECO"));
 
    interpolate_ = iConfig.getParameter<bool>("interpolate");
    sumRecHits_ = iConfig.getParameter<bool>("sumRecHits");
 
    std::string ifname = "RecoHI/HiJetAlgos/data/PU_DATA.root";
    TFile* inf = new TFile(edm::FileInPath(ifname).fullPath().data());
    fPU = (TF1*)inf->Get("fPU");
    fMean = (TF1*)inf->Get("fMean");
    fRMS = (TF1*)inf->Get("fRMS");
    hC = (TH1D*)inf->Get("hC");
 
    for(int i = 0; i < 40; ++i){
       hEta.push_back((TH1D*)inf->Get(Form("hEta_%d",i)));
       hEtaMean.push_back((TH1D*)inf->Get(Form("hEtaMean_%d",i)));
       hEtaRMS.push_back((TH1D*)inf->Get(Form("hEtaRMS_%d",i)));
    }
 
 }

ParametrizedSubtractor::~ParametrizedSubtractor ( )

inline

Definition at line 29 of file ParametrizedSubtractor.h.

29 {;}

Member Function Documentation

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

virtual

Reimplemented from PileUpSubtractor.

Definition at line 147 of file ParametrizedSubtractor.cc.

References PileUpSubtractor::fjInputs_.

 {
    orphanInput = *fjInputs_;
 }

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

virtual

Reimplemented from PileUpSubtractor.

Definition at line 98 of file ParametrizedSubtractor.cc.

                                                                                      {
    return;
 }

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

Definition at line 223 of file ParametrizedSubtractor.cc.

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

Referenced by offsetCorrectJets(), and subtractPedestal().

                                                                       {
    const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
    const GlobalPoint& pos=geo_->getPosition(ctc->id());
    double energy = ctc->emEnergy() + ctc->hadEnergy();
 
    if(0){
       energy = 0;
       const std::vector<DetId>& hitids = ctc->constituents();
       for(unsigned int i = 0; i< hitids.size(); ++i){
 
       }
    }
 
    double et = energy*sin(pos.theta());
    return et;
 }

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

Definition at line 240 of file ParametrizedSubtractor.cc.

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

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

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

virtual

Reimplemented from PileUpSubtractor.

Definition at line 247 of file ParametrizedSubtractor.cc.

References getPU(), and PileUpSubtractor::ieta().

                                                                               {
    int it = ieta(in);
    return getPU(it,1,0);
 }

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

virtual

Reimplemented from PileUpSubtractor.

Definition at line 257 of file ParametrizedSubtractor.cc.

References getPU(), and PileUpSubtractor::ieta().

                                                                                  {
    int it = ieta(in);
    return getPU(it,1,1);
 }

double ParametrizedSubtractor::getPU	(	int	ieta,
		bool	addMean,
		bool	addSigma
	)		const

Definition at line 262 of file ParametrizedSubtractor.cc.

References bin_, centrality_, fMean, fRMS, hC, hEtaMean, hEtaRMS, interpolate_, n, and PileUpSubtractor::nSigmaPU_.

Referenced by getMeanAtTower(), getPileUpAtTower(), getSigmaAtTower(), offsetCorrectJets(), and subtractPedestal().

                                                                                {   
 
   //double e = hEta[bin_]->GetBinContent(hEta[bin_]->FindBin(ieta));
   //double c = fPU->Eval(centrality_);
 
    double em = hEtaMean[bin_]->GetBinContent(hEtaMean[bin_]->FindBin(ieta));
    double cm = fMean->Eval(centrality_);
 
    double er = hEtaRMS[bin_]->GetBinContent(hEtaRMS[bin_]->FindBin(ieta));
    double cr = fRMS->Eval(centrality_);
 
    if(interpolate_){
       double n = 0;
       int hbin = 40-bin_;
       double centerweight =  (centrality_ - hC->GetBinCenter(hbin));
       double lowerweight = (centrality_ - hC->GetBinLowEdge(hbin));
       double upperweight = (centrality_ - hC->GetBinLowEdge(hbin+1));
 
       em *= lowerweight*upperweight;
       er *= lowerweight*upperweight;
       n += lowerweight*upperweight;
 
       if(bin_ > 0){
      em += upperweight*centerweight*hEtaMean[bin_]->GetBinContent(hEtaMean[bin_-1]->FindBin(ieta));
      er += upperweight*centerweight*hEtaRMS[bin_]->GetBinContent(hEtaRMS[bin_-1]->FindBin(ieta));
      n += upperweight*centerweight;
       }
 
       if(bin_ < 39){
      em += lowerweight*centerweight*hEtaMean[bin_]->GetBinContent(hEtaMean[bin_+1]->FindBin(ieta));
          er += lowerweight*centerweight*hEtaRMS[bin_]->GetBinContent(hEtaRMS[bin_+1]->FindBin(ieta));
      n += lowerweight*centerweight;
       }
       em /= n;
       er /= n;
    }
 
    //   return e*c;
    return addMean*em*cm + addSigma*nSigmaPU_*er*cr;
 }

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

virtual

Reimplemented from PileUpSubtractor.

Definition at line 252 of file ParametrizedSubtractor.cc.

References getPU(), and PileUpSubtractor::ieta().

                                                                                 {
    int it = ieta(in);
    return getPU(it,0,1);
 }

void ParametrizedSubtractor::offsetCorrectJets ( )

virtual

Reimplemented from PileUpSubtractor.

Definition at line 154 of file ParametrizedSubtractor.cc.

References def, PileUpSubtractor::doAreaFastjet_, PileUpSubtractor::doRhoFastjet_, PileUpSubtractor::fjActiveArea_, PileUpSubtractor::fjClusterSeq_, PileUpSubtractor::fjInputs_, PileUpSubtractor::fjJets_, PileUpSubtractor::fjOriginalInputs_, getEt(), getPU(), PileUpSubtractor::ieta(), PileUpSubtractor::jetOffset_, PileUpSubtractor::jetPtMin_, LogDebug, PileUpSubtractor::nSigmaPU_, dt_dqm_sourceclient_common_cff::reco, rescaleRMS(), subtractPedestal(), and sumRecHits_.

 {
 
   LogDebug("PileUpSubtractor")<<"The subtractor correcting jets...\n";
   jetOffset_.clear();
 
   using namespace reco;
   
   (*fjInputs_) = fjOriginalInputs_;
   rescaleRMS(nSigmaPU_);
   subtractPedestal(*fjInputs_);
 
   if(0){
      const fastjet::JetDefinition& def = fjClusterSeq_->jet_def();
      if ( !doAreaFastjet_ && !doRhoFastjet_) {
     fastjet::ClusterSequence newseq( *fjInputs_, def );
     (*fjClusterSeq_) = newseq;
      } else {
     fastjet::ClusterSequenceArea newseq( *fjInputs_, def , *fjActiveArea_ );
     (*fjClusterSeq_) = newseq;
      }
      
      (*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(fjClusterSeq_->constituents(*pseudojetTMP));
     
     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();
       
     if(sumRecHits_){
        Original_Et = getEt(originalTower);
     }
 
     double etnew = Original_Et - getPU(it,1,1);
     if(etnew < 0.) etnew = 0;
     newjetet = newjetet + etnew;
     jetOffset_[ijet] += Original_Et - etnew;
       }
 
     if(sumRecHits_){       
        double mScale = newjetet/pseudojetTMP->Et();
        int cshist = pseudojetTMP->cluster_hist_index();
        pseudojetTMP->reset(pseudojetTMP->px()*mScale, pseudojetTMP->py()*mScale,
                pseudojetTMP->pz()*mScale, pseudojetTMP->e()*mScale);
        pseudojetTMP->set_cluster_hist_index(cshist);
     }
 
   }
 }

void ParametrizedSubtractor::rescaleRMS ( double s )

Definition at line 14 of file ParametrizedSubtractor.cc.

Referenced by offsetCorrectJets().

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

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

virtual

Reimplemented from PileUpSubtractor.

Definition at line 48 of file ParametrizedSubtractor.cc.

References PileUpSubtractor::allgeomid_, bin_, centrality_, centTag_, PileUpSubtractor::emean_, PileUpSubtractor::esigma_, PileUpSubtractor::geo_, PileUpSubtractor::geomtowers_, edm::EventSetup::get(), edm::Event::getByLabel(), CaloGeometry::getValidDetIds(), hC, DetId::Hcal, i, PileUpSubtractor::ieta(), PileUpSubtractor::ietamax_, PileUpSubtractor::ietamin_, LogDebug, PileUpSubtractor::ntowersWithJets_, and edm::ESHandle< class >::product().

                                                                                          {
    LogDebug("PileUpSubtractor")<<"The subtractor setting up geometry...\n";
 
    //   if(!cbins_) getCentralityBinsFromDB(iSetup);
 
    edm::Handle<reco::Centrality> cent;
    iEvent.getByLabel(centTag_,cent);
    
    centrality_ = cent->EtHFhitSum();
    bin_ = 40-hC->FindBin(centrality_);
    if(bin_ > 39) bin_ = 39;
    if(bin_ < 0) bin_ = 0;
 
    if(!geo_) {
       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++) {
       emean_[i] = 0.;
       esigma_[i] = 0.;
       ntowersWithJets_[i] = 0;
    }
 }

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

virtual

Reimplemented from PileUpSubtractor.

Definition at line 102 of file ParametrizedSubtractor.cc.

References dropZeroTowers_, getEt(), getPU(), PileUpSubtractor::ieta(), getHLTprescales::index, PileUpSubtractor::iphi(), LogDebug, and sumRecHits_.

Referenced by offsetCorrectJets().

 {
    if(0){
       return;
    }else{
       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 - getPU(it,1,1);
       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 ( 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;
 
    }
 }