ReflectedIterator.cc

Go to the documentation of this file.

#include "RecoHI/HiJetAlgos/interface/ReflectedIterator.h"
#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "DataFormats/Candidate/interface/Candidate.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
using namespace std;

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


void ReflectedIterator::offsetCorrectJets()
{

  LogDebug("PileUpSubtractor")<<"The subtractor correcting jets...\n";
  jetOffset_.clear();

  using namespace reco;
  
  (*fjInputs_) = fjOriginalInputs_;
  rescaleRMS(nSigmaPU_);
  subtractPedestal(*fjInputs_);
  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();
    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 ReflectedIterator::subtractPedestal(vector<fastjet::PseudoJet> & coll)
{

   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 ( 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;

}

void ReflectedIterator::calculatePedestal( vector<fastjet::PseudoJet> const & coll )
{
   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 ReflectedIterator::getEt(const reco::CandidatePtr & in) const {
   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 ReflectedIterator::getEta(const reco::CandidatePtr & in) const {
   const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
   const GlobalPoint& pos=geo_->getPosition(ctc->id());
   double eta = pos.eta();
   return eta;
}