#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 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_

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_, runTheMatrix::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.

{;}

Member Function Documentation

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

Reimplemented from PileUpSubtractor.

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

Reimplemented from PileUpSubtractor.

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

Definition at line 224 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().

                                                                      {
   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 241 of file ParametrizedSubtractor.cc.

References eta(), PV3DBase< T, PVType, FrameType >::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 248 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 258 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 263 of file ParametrizedSubtractor.cc.

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

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

                                                                               {   

  //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 253 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 155 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< T >::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;
   }
}