#include <JetIDHelper.h>

Classes
struct	subtower
Public Member Functions
double	approximatefHPD () const
double	approximatefRBX () const
void	calculate (const edm::Event &event, const reco::CaloJet &jet, const int iDbg=0)
double	fEB () const
double	fEE () const
double	fHB () const
double	fHE () const
double	fHFOOT () const
double	fHO () const
double	fHPD () const
void	fillDescription (edm::ParameterSetDescription &iDesc)
double	fLong () const
double	fLSbad () const
double	fRBX () const
double	fShort () const
double	fSubDetector1 () const
double	fSubDetector2 () const
double	fSubDetector3 () const
double	fSubDetector4 () const
int	hitsInN90 () const
void	initValues ()
	JetIDHelper ()
	JetIDHelper (edm::ParameterSet const &pset)
int	n90Hits () const
int	nECALTowers () const
int	nHCALTowers () const
double	restrictedEMF () const
	~JetIDHelper ()
Private Types
enum	Region { unknown_region = -1, HFneg, HEneg, HBneg, HBpos, HEpos, HFpos }
Private Member Functions
void	classifyJetComponents (const edm::Event &event, const reco::CaloJet &jet, std::vector< double > &energies, std::vector< double > &subdet_energies, std::vector< double > &Ecal_energies, std::vector< double > &Hcal_energies, std::vector< double > &HO_energies, std::vector< double > &HPD_energies, std::vector< double > &RBX_energies, double &LS_bad_energy, double &HF_OOT_energy, const int iDbg=0)
void	classifyJetTowers (const edm::Event &event, const reco::CaloJet &jet, std::vector< subtower > &subtowers, std::vector< subtower > &Ecal_subtowers, std::vector< subtower > &Hcal_subtowers, std::vector< subtower > &HO_subtowers, std::vector< double > &HPD_energies, std::vector< double > &RBX_energies, const int iDbg=0)
int	HBHE_oddness (int iEta, int depth)
int	HBHE_oddness (int iEta)
Region	HBHE_region (int iEta, int depth)
unsigned int	hitsInNCarrying (double fraction, std::vector< subtower > descending_towers)
unsigned int	nCarrying (double fraction, std::vector< double > descending_energies)
Region	region (int iEta)
Private Attributes
double	approximatefHPD_
double	approximatefRBX_
double	approximatefSubDetector1_
double	approximatefSubDetector2_
double	approximatefSubDetector3_
double	approximatefSubDetector4_
edm::InputTag	ebRecHitsColl_
edm::InputTag	eeRecHitsColl_
double	fEB_
double	fEE_
double	fHB_
double	fHE_
double	fHFOOT_
double	fHO_
double	fHPD_
double	fLong_
double	fLS_
double	fRBX_
double	fShort_
double	fSubDetector1_
double	fSubDetector2_
double	fSubDetector3_
double	fSubDetector4_
edm::InputTag	hbheRecHitsColl_
edm::InputTag	hfRecHitsColl_
int	hitsInN90_
edm::InputTag	hoRecHitsColl_
int	n90Hits_
int	nECALTowers_
int	nHCALTowers_
double	restrictedEMF_
bool	useRecHits_
Static Private Attributes
static int	sanity_checks_left_ = 100

Detailed Description

Definition at line 15 of file JetIDHelper.h.

Member Enumeration Documentation

enum reco::helper::JetIDHelper::Region [private]

Enumerator:

unknown_region
HFneg
HEneg
HBneg
HBpos
HEpos
HFpos

Definition at line 90 of file JetIDHelper.h.

                 {
        unknown_region = -1,
        HFneg, HEneg, HBneg, HBpos, HEpos, HFpos };

Constructor & Destructor Documentation

reco::helper::JetIDHelper::JetIDHelper ( ) [inline]

Definition at line 19 of file JetIDHelper.h.

{}

reco::helper::JetIDHelper::JetIDHelper ( edm::ParameterSet const & pset )

Definition at line 40 of file JetIDHelper.cc.

References edm::ParameterSet::getParameter().

{
  useRecHits_ = pset.getParameter<bool>("useRecHits");
  if( useRecHits_ ) {
    hbheRecHitsColl_ = pset.getParameter<edm::InputTag>("hbheRecHitsColl");
    hoRecHitsColl_   = pset.getParameter<edm::InputTag>("hoRecHitsColl");
    hfRecHitsColl_   = pset.getParameter<edm::InputTag>("hfRecHitsColl");
    ebRecHitsColl_   = pset.getParameter<edm::InputTag>("ebRecHitsColl");
    eeRecHitsColl_   = pset.getParameter<edm::InputTag>("eeRecHitsColl");   
  }
  initValues();
}

reco::helper::JetIDHelper::~JetIDHelper ( ) [inline]

Definition at line 21 of file JetIDHelper.h.

{}

Member Function Documentation

double reco::helper::JetIDHelper::approximatefHPD ( ) const [inline]

Definition at line 54 of file JetIDHelper.h.

References approximatefHPD_.

Referenced by JetIDProducer::produce().

{ return    approximatefHPD_;}

double reco::helper::JetIDHelper::approximatefRBX ( ) const [inline]

Definition at line 55 of file JetIDHelper.h.

References approximatefRBX_.

Referenced by JetIDProducer::produce().

{ return    approximatefRBX_;}

void reco::helper::JetIDHelper::calculate	(	const edm::Event &	event,
		const reco::CaloJet &	jet,
		const int	iDbg = `0`
	)

Definition at line 86 of file JetIDHelper.cc.

References abs, gather_cfg::cout, reco::CaloJet::emEnergyInEB(), reco::CaloJet::emEnergyInEE(), reco::CaloJet::emEnergyInHF(), reco::LeafCandidate::energy(), reco::LeafCandidate::et(), reco::CaloJet::hadEnergyInHB(), reco::CaloJet::hadEnergyInHE(), reco::CaloJet::hadEnergyInHF(), reco::CaloJet::hadEnergyInHO(), reco::helper::hasNonPositiveE(), i, siStripFEDMonitor_P5_cff::Max, and reco::LeafCandidate::pt().

Referenced by JetMETHLTOfflineSource::analyze(), HLTJets::analyze(), HLTCaloJetIDProducer::produce(), JetIDProducer::produce(), BTagHLTOfflineSource::selectJets(), and FourVectorHLTOffline::selectJets().

{
  initValues();

  // --------------------------------------------------
  // 1) jet ID variable derived from existing fractions
  // --------------------------------------------------

  double E_EM = TMath::Max( float(0.), jet.emEnergyInHF() )
              + TMath::Max( float(0.), jet.emEnergyInEB() )
              + TMath::Max( float(0.), jet.emEnergyInEE() );
  double E_Had = TMath::Max( float(0.), jet.hadEnergyInHB() )
               + TMath::Max( float(0.), jet.hadEnergyInHE() )
               + TMath::Max( float(0.), jet.hadEnergyInHO() )
               + TMath::Max( float(0.), jet.hadEnergyInHF() );
  if( E_Had + E_EM > 0 ) restrictedEMF_ = E_EM / ( E_EM + E_Had );
  if( iDbg > 1 ) cout<<"jet pT: "<<jet.pt()<<", eT: "<<jet.et()<<", E: "
                     <<jet.energy()<<" rEMF: "<<restrictedEMF_<<endl;

  // ------------------------
  // 2) tower based variables
  // ------------------------
  vector<subtower> subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers;
  vector<double> HPD_energies, RBX_energies;

  classifyJetTowers( event, jet, 
                     subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers, 
                     HPD_energies, RBX_energies, iDbg );
  if( iDbg > 1 ) {
    cout<<"E:";
    for (unsigned int i=0; i<subtowers.size(); ++i) cout<<" "<<subtowers[i].E<<","<<subtowers[i].Nhit;
    cout<<"\nECal_E:";
    for (unsigned int i=0; i<Ecal_subtowers.size(); ++i) cout<<" "<<Ecal_subtowers[i].E<<","<<Ecal_subtowers[i].Nhit;
    cout<<"\nHCal_E:";
    for (unsigned int i=0; i<Hcal_subtowers.size(); ++i) cout<<" "<<Hcal_subtowers[i].E<<","<<Hcal_subtowers[i].Nhit;
    cout<<"\nHO_E:";
    for (unsigned int i=0; i<HO_subtowers.size(); ++i) cout<<" "<<HO_subtowers[i].E<<","<<HO_subtowers[i].Nhit;
    cout<<"\nHPD_E:";
    for (unsigned int i=0; i<HPD_energies.size(); ++i) cout<<" "<<HPD_energies[i];
    cout<<"\nRBX_E:";
    for (unsigned int i=0; i<RBX_energies.size(); ++i) cout<<" "<<RBX_energies[i];
    cout<<endl;
  }

  // counts
  hitsInN90_ = hitsInNCarrying( 0.9, subtowers );
  nHCALTowers_ = Hcal_subtowers.size();
  vector<subtower>::const_iterator it;
  it = find_if( Ecal_subtowers.begin(), Ecal_subtowers.end(), hasNonPositiveE );
  nECALTowers_ = it - Ecal_subtowers.begin(); // ignores negative energies from HF!

  // energy fractions
  if( jet.energy() > 0 ) {
    if( HPD_energies.size() > 0 ) approximatefHPD_ = HPD_energies.at( 0 ) / jet.energy();
    if( RBX_energies.size() > 0 ) approximatefRBX_ = RBX_energies.at( 0 ) / jet.energy();
  }

  // -----------------------
  // 3) cell based variables
  // -----------------------
  if( useRecHits_ ) {
    vector<double> energies, subdet_energies, Ecal_energies, Hcal_energies, HO_energies;
    double LS_bad_energy, HF_OOT_energy;
    classifyJetComponents( event, jet, 
                           energies, subdet_energies, Ecal_energies, Hcal_energies, HO_energies, 
                           HPD_energies, RBX_energies, LS_bad_energy, HF_OOT_energy, iDbg );

    // counts
    n90Hits_ = nCarrying( 0.9, energies );

    // energy fractions
    if( jet.energy() > 0 ) {
      if( HPD_energies.size() > 0 ) fHPD_ = HPD_energies.at( 0 ) / jet.energy();
      if( RBX_energies.size() > 0 ) fRBX_ = RBX_energies.at( 0 ) / jet.energy();
      if( subdet_energies.size() > 0 ) fSubDetector1_ = subdet_energies.at( 0 ) / jet.energy();
      if( subdet_energies.size() > 1 ) fSubDetector2_ = subdet_energies.at( 1 ) / jet.energy();
      if( subdet_energies.size() > 2 ) fSubDetector3_ = subdet_energies.at( 2 ) / jet.energy();
      if( subdet_energies.size() > 3 ) fSubDetector4_ = subdet_energies.at( 3 ) / jet.energy();
      fLS_ = LS_bad_energy / jet.energy();
      fHFOOT_ = HF_OOT_energy / jet.energy();

      if( sanity_checks_left_ > 0 ) {
        --sanity_checks_left_;
        double EH_sum = accumulate( Ecal_energies.begin(), Ecal_energies.end(), 0. );
        EH_sum = accumulate( Hcal_energies.begin(), Hcal_energies.end(), EH_sum );
        double EHO_sum = accumulate( HO_energies.begin(), HO_energies.end(), EH_sum );
        if( jet.energy() > 0.001 && abs( EH_sum/jet.energy() - 1 ) > 0.01 && abs( EHO_sum/jet.energy() - 1 ) > 0.01 )
          edm::LogWarning("BadInput")<<"Jet energy ("<<jet.energy()
                                     <<") does not match the total energy in the recHits ("<<EHO_sum
                                     <<", or without HO: "<<EH_sum<<") . Are these the right recHits? "
                                     <<"Were jet energy corrections mistakenly applied before jet ID? A bug?";
        if( iDbg > 1 ) cout<<"Sanity check - E: "<<jet.energy()<<" =? EH_sum: "<<EH_sum<<" / EHO_sum: "<<EHO_sum<<endl;
      }
    }

    if( iDbg > 1 ) {
      cout<<"DBG - fHPD: "<<fHPD_<<", fRBX: "<<fRBX_<<", nh90: "<<n90Hits_<<", fLS: "<<fLS_<<", fHFOOT: "<<fHFOOT_<<endl;
      cout<<"    -~fHPD: "<<approximatefHPD_<<", ~fRBX: "<<approximatefRBX_
          <<", hits in n90: "<<hitsInN90_<<endl;
      cout<<"    - nHCALTowers: "<<nHCALTowers_<<", nECALTowers: "<<nECALTowers_
          <<"; subdet fractions: "<<fSubDetector1_<<", "<<fSubDetector2_<<", "<<fSubDetector3_<<", "<<fSubDetector4_<<endl;
    }
  }
}

void reco::helper::JetIDHelper::classifyJetComponents	(	const edm::Event &	event,
		const reco::CaloJet &	jet,
		std::vector< double > &	energies,
		std::vector< double > &	subdet_energies,
		std::vector< double > &	Ecal_energies,
		std::vector< double > &	Hcal_energies,
		std::vector< double > &	HO_energies,
		std::vector< double > &	HPD_energies,
		std::vector< double > &	RBX_energies,
		double &	LS_bad_energy,
		double &	HF_OOT_energy,
		const int	iDbg = `0`
	)		`[private]`

Definition at line 225 of file JetIDHelper.cc.

References abs, gather_cfg::cout, HcalDetId::depth(), egHLT::errCodes::EBRecHits, DetId::Ecal, egHLT::errCodes::EERecHits, reco::LeafCandidate::energy(), flags, reco::CaloJet::getCaloConstituents(), egHLT::errCodes::HBHERecHits, DetId::Hcal, HcalForward, HcalOuter, HcalCaloFlagLabels::HFDigiTime, HcalCaloFlagLabels::HFLongShort, egHLT::errCodes::HFRecHits, HcalCaloFlagLabels::HFTimingTrustBits, HcalDetId::ieta(), reco::helper::select2nd(), python::multivaluedict::sort(), HcalDetId::subdet(), HcalCaloFlagLabels::TimingAddedBit, HcalCaloFlagLabels::TimingErrorBit, HcalCaloFlagLabels::TimingSubtractedBit, and create_public_pileup_plots::transform.

{
  energies.clear(); subdet_energies.clear(); Ecal_energies.clear(); Hcal_energies.clear(); HO_energies.clear();
  HPD_energies.clear(); RBX_energies.clear();
  LS_bad_energy = HF_OOT_energy = 0.;
  
  std::map< int, double > HPD_energy_map, RBX_energy_map;
  vector< double > EB_energies, EE_energies, HB_energies, HE_energies, short_energies, long_energies;
  edm::Handle<HBHERecHitCollection> HBHERecHits;
  edm::Handle<HORecHitCollection> HORecHits;
  edm::Handle<HFRecHitCollection> HFRecHits;
  edm::Handle<EBRecHitCollection> EBRecHits;
  edm::Handle<EERecHitCollection> EERecHits;
  // the jet only contains DetIds, so first read recHit collection
  event.getByLabel( hbheRecHitsColl_, HBHERecHits );
  event.getByLabel( hoRecHitsColl_, HORecHits );
  event.getByLabel( hfRecHitsColl_, HFRecHits );
  event.getByLabel( ebRecHitsColl_, EBRecHits );
  event.getByLabel( eeRecHitsColl_, EERecHits );
  if( iDbg > 2 ) cout<<"# of rechits found - HBHE: "<<HBHERecHits->size()
                     <<", HO: "<<HORecHits->size()<<", HF: "<<HFRecHits->size()
                     <<", EB: "<<EBRecHits->size()<<", EE: "<<EERecHits->size()<<endl;

  vector< CaloTowerPtr > towers = jet.getCaloConstituents ();
  int nTowers = towers.size();
  if( iDbg > 9 ) cout<<"In classifyJetComponents. # of towers found: "<<nTowers<<endl;

  for( int iTower = 0; iTower <nTowers ; iTower++ ) {

    CaloTowerPtr& tower = towers[iTower];

    int nCells = tower->constituentsSize();
    if( iDbg ) cout<<"tower #"<<iTower<<" has "<<nCells<<" cells. "
                   <<"It's at iEta: "<<tower->ieta()<<", iPhi: "<<tower->iphi()<<endl;

    const vector< DetId >& cellIDs = tower->constituents();  // cell == recHit
  
    for( int iCell = 0; iCell < nCells; ++iCell ) {
      DetId::Detector detNum = cellIDs[iCell].det();
      if( detNum == DetId::Hcal ) {
        HcalDetId HcalID = cellIDs[ iCell ];
        HcalSubdetector HcalNum = HcalID.subdet();
        double hitE = 0;
        if( HcalNum == HcalOuter ) {
          HORecHitCollection::const_iterator theRecHit=HORecHits->find(HcalID);
          if (theRecHit == HORecHits->end()) {
            edm::LogWarning("UnexpectedEventContents")<<"Can't find the HO recHit with ID: "<<HcalID;
            continue;
          }
          hitE = theRecHit->energy();
          HO_energies.push_back( hitE );
            
        } else if( HcalNum == HcalForward ) {
          
          HFRecHitCollection::const_iterator theRecHit=HFRecHits->find( HcalID );           
          if( theRecHit == HFRecHits->end() ) {
            edm::LogWarning("UnexpectedEventContents")<<"Can't find the HF recHit with ID: "<<HcalID;
            continue;
          }
          hitE = theRecHit->energy();
          if( iDbg>4 ) cout 
                         << "hit #"<<iCell<<" is  HF , E: "<<hitE<<" iEta: "<<theRecHit->id().ieta()
                         <<", depth: "<<theRecHit->id().depth()<<", iPhi: "
                         <<theRecHit->id().iphi();
          
          if( HcalID.depth() == 1 ) 
            long_energies.push_back( hitE );
          else
            short_energies.push_back( hitE );
          
          uint32_t flags = theRecHit->flags();
          if( flags & (1<<HcalCaloFlagLabels::HFLongShort) ) LS_bad_energy += hitE;
          if( flags & ( (1<<HcalCaloFlagLabels::HFDigiTime) | 
                        (3<<HcalCaloFlagLabels::HFTimingTrustBits) | 
                        (1<<HcalCaloFlagLabels::TimingSubtractedBit) |
                        (1<<HcalCaloFlagLabels::TimingAddedBit) |
                        (1<<HcalCaloFlagLabels::TimingErrorBit) ) ) HF_OOT_energy += hitE;
          if( iDbg>4 && flags ) cout<<"flags: "<<flags
                                    <<" -> LS_bad_energy: "<<LS_bad_energy<<", HF_OOT_energy: "<<HF_OOT_energy<<endl; 

        } else { // HBHE
          
          HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find( HcalID );     
          if( theRecHit == HBHERecHits->end() ) {
            edm::LogWarning("UnexpectedEventContents")<<"Can't find the HBHE recHit with ID: "<<HcalID; 
            continue;
          }
          hitE = theRecHit->energy();
          int iEta = theRecHit->id().ieta();
          int depth = theRecHit->id().depth();
          Region region = HBHE_region( iEta, depth );
          int hitIPhi = theRecHit->id().iphi();
          if( iDbg>3 ) cout<<"hit #"<<iCell<<" is HBHE, E: "<<hitE<<" iEta: "<<iEta
                           <<", depth: "<<depth<<", iPhi: "<<theRecHit->id().iphi()
                           <<" -> "<<region;
          int absIEta = TMath::Abs( theRecHit->id().ieta() );
          if( depth == 3 && (absIEta == 28 || absIEta == 29) ) {
            hitE /= 2; // Depth 3 at the HE forward edge is split over tower 28 & 29, and jet reco. assigns half each
          }
          int iHPD = 100 * region;
          int iRBX = 100 * region + ((hitIPhi + 1) % 72) / 4; // 71,72,1,2 are in the same RBX module
          
          if( std::abs( iEta ) >= 21 ) {
            if( (0x1 & hitIPhi) == 0 ) {
              edm::LogError("CodeAssumptionsViolated")<<"Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
              return;
            }
            bool oddnessIEta = HBHE_oddness( iEta, depth );
            bool upperIPhi = (( hitIPhi%4 ) == 1 || ( hitIPhi%4 ) == 2); // the upper iPhi indices in the HE wedge
            // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
            // 1) in the upper iPhis of the module, the even iEtas belong to the higher iPhi
            // 2) in the loewr iPhis of the module, the odd  iEtas belong to the higher iPhi
            if( upperIPhi != oddnessIEta ) ++hitIPhi; 
            // note that hitIPhi could not be 72 before, so it's still in the legal range [1,72]
          }
          iHPD += hitIPhi;
          
          // book the energies
          HPD_energy_map[ iHPD ] += hitE;
          RBX_energy_map[ iRBX ] += hitE;
          if( iDbg > 5 ) cout<<" --> H["<<iHPD<<"]="<<HPD_energy_map[iHPD]
                             <<", R["<<iRBX<<"]="<<RBX_energy_map[iRBX];
          if( iDbg > 1 ) cout<<endl;
          
          if( region == HBneg || region == HBpos )
            HB_energies.push_back( hitE );
          else
            HE_energies.push_back( hitE );
          
        } // if HBHE
        if( hitE == 0 ) edm::LogWarning("UnexpectedEventContents")<<"HCal hitE==0? (or unknown subdetector?)";
        
      } // if HCAL 
        
      else if (detNum == DetId::Ecal) {
        
        int EcalNum =  cellIDs[iCell].subdetId();
        double hitE = 0;
        if( EcalNum == 1 ){
          EBDetId EcalID = cellIDs[iCell];
          EBRecHitCollection::const_iterator theRecHit=EBRecHits->find(EcalID);
          if( theRecHit == EBRecHits->end() ) {edm::LogWarning("UnexpectedEventContents")
              <<"Can't find the EB recHit with ID: "<<EcalID; continue;}
          hitE = theRecHit->energy();
          EB_energies.push_back( hitE );
        } else if(  EcalNum == 2 ){
          EEDetId EcalID = cellIDs[iCell];
          EERecHitCollection::const_iterator theRecHit=EERecHits->find(EcalID);
          if( theRecHit == EERecHits->end() ) {edm::LogWarning("UnexpectedEventContents")
              <<"Can't find the EE recHit with ID: "<<EcalID; continue;}
          hitE = theRecHit->energy();
          EE_energies.push_back( hitE );
        }
        if( hitE == 0 ) edm::LogWarning("UnexpectedEventContents")<<"ECal hitE==0? (or unknown subdetector?)";
        if( iDbg > 6 ) cout<<"EcalNum: "<<EcalNum<<" hitE: "<<hitE<<endl;
      } // 
    } // loop on cells
  } // loop on towers

  /* Disabling check until HO is accounted for in EMF. Check was used in CMSSW_2, where HE was excluded.
  double expHcalE = jet.energy() * (1-jet.emEnergyFraction());
  if( totHcalE + expHcalE > 0 && 
      TMath::Abs( totHcalE - expHcalE ) > 0.01 && 
      ( totHcalE - expHcalE ) / ( totHcalE + expHcalE ) > 0.0001 ) {
    edm::LogWarning("CodeAssumptionsViolated")<<"failed to account for all Hcal energies"
                                              <<totHcalE<<"!="<<expHcalE;
                                              } */   
  // concatenate Hcal and Ecal lists
  Hcal_energies.insert( Hcal_energies.end(), HB_energies.begin(), HB_energies.end() );
  Hcal_energies.insert( Hcal_energies.end(), HE_energies.begin(), HE_energies.end() );
  Hcal_energies.insert( Hcal_energies.end(), short_energies.begin(), short_energies.end() );
  Hcal_energies.insert( Hcal_energies.end(), long_energies.begin(), long_energies.end() );
  Ecal_energies.insert( Ecal_energies.end(), EB_energies.begin(), EB_energies.end() );
  Ecal_energies.insert( Ecal_energies.end(), EE_energies.begin(), EE_energies.end() );

  // sort the energies
  std::sort( Hcal_energies.begin(), Hcal_energies.end(), greater<double>() );
  std::sort( Ecal_energies.begin(), Ecal_energies.end(), greater<double>() );

  // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
  std::transform( HPD_energy_map.begin(), HPD_energy_map.end(), 
                  std::inserter (HPD_energies, HPD_energies.end()), select2nd ); 
  //              std::select2nd<std::map<int,double>::value_type>());
  std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
  std::transform( RBX_energy_map.begin(), RBX_energy_map.end(), 
                  std::inserter (RBX_energies, RBX_energies.end()), select2nd );
  //              std::select2nd<std::map<int,double>::value_type>());
  std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );

  energies.insert( energies.end(), Hcal_energies.begin(), Hcal_energies.end() );
  energies.insert( energies.end(), Ecal_energies.begin(), Ecal_energies.end() );
  energies.insert( energies.end(), HO_energies.begin(), HO_energies.end() );
  std::sort( energies.begin(), energies.end(), greater<double>() );

  // prepare sub detector energies, then turn them into fractions
  fEB_ = std::accumulate( EB_energies.begin(), EB_energies.end(), 0. );
  fEE_ = std::accumulate( EE_energies.begin(), EE_energies.end(), 0. );
  fHB_ = std::accumulate( HB_energies.begin(), HB_energies.end(), 0. );
  fHE_ = std::accumulate( HE_energies.begin(), HE_energies.end(), 0. );
  fHO_ = std::accumulate( HO_energies.begin(), HO_energies.end(), 0. );
  fShort_ = std::accumulate( short_energies.begin(), short_energies.end(), 0. );
  fLong_ = std::accumulate( long_energies.begin(), long_energies.end(), 0. );
  subdet_energies.push_back( fEB_ );
  subdet_energies.push_back( fEE_ );
  subdet_energies.push_back( fHB_ );
  subdet_energies.push_back( fHE_ );
  subdet_energies.push_back( fHO_ );
  subdet_energies.push_back( fShort_ );
  subdet_energies.push_back( fLong_ );
  std::sort( subdet_energies.begin(), subdet_energies.end(), greater<double>() );
  if( jet.energy() > 0 ) {
    fEB_ /= jet.energy();
    fEE_ /= jet.energy();
    fHB_ /= jet.energy();
    fHE_ /= jet.energy();
    fHO_ /= jet.energy();
    fShort_ /= jet.energy();
    fLong_ /= jet.energy();
  } else {
    if( fEB_ > 0 || fEE_ > 0 || fHB_ > 0 || fHE_ > 0 || fHO_ > 0 || fShort_ > 0 || fLong_ > 0 )
      edm::LogError("UnexpectedEventContents")<<"Jet ID Helper found energy in subdetectors and jet E <= 0";
  }      
}

void reco::helper::JetIDHelper::classifyJetTowers	(	const edm::Event &	event,
		const reco::CaloJet &	jet,
		std::vector< subtower > &	subtowers,
		std::vector< subtower > &	Ecal_subtowers,
		std::vector< subtower > &	Hcal_subtowers,
		std::vector< subtower > &	HO_subtowers,
		std::vector< double > &	HPD_energies,
		std::vector< double > &	RBX_energies,
		const int	iDbg = `0`
	)		`[private]`

Definition at line 457 of file JetIDHelper.cc.

References abs, gather_cfg::cout, DetId::Ecal, reco::CaloJet::getCaloConstituents(), DetId::Hcal, HcalOuter, reco::helper::select2nd(), python::multivaluedict::sort(), HcalDetId::subdet(), reco::helper::subtower_has_greater_E(), and create_public_pileup_plots::transform.

{
  subtowers.clear(); Ecal_subtowers.clear(); Hcal_subtowers.clear(); HO_subtowers.clear();
  HPD_energies.clear(); RBX_energies.clear();

  std::map< int, double > HPD_energy_map, RBX_energy_map;

  vector< CaloTowerPtr > towers = jet.getCaloConstituents ();
  int nTowers = towers.size();
  if( iDbg > 9 ) cout<<"classifyJetTowers started. # of towers found: "<<nTowers<<endl;

  for( int iTower = 0; iTower <nTowers ; iTower++ ) {

    CaloTowerPtr& tower = towers[iTower];

    int nEM = 0, nHad = 0, nHO = 0;
    const vector< DetId >& cellIDs = tower->constituents();  // cell == recHit
    int nCells = cellIDs.size();
    if( iDbg ) cout<<"tower #"<<iTower<<" has "<<nCells<<" cells. "
                   <<"It's at iEta: "<<tower->ieta()<<", iPhi: "<<tower->iphi()<<endl;
  
    for( int iCell = 0; iCell < nCells; ++iCell ) {
      DetId::Detector detNum = cellIDs[iCell].det();
      if( detNum == DetId::Hcal ) {
        HcalDetId HcalID = cellIDs[ iCell ];
        HcalSubdetector HcalNum = HcalID.subdet();
        if( HcalNum == HcalOuter ) {
          ++nHO;
        } else {
          ++nHad;
        }
      } else if (detNum == DetId::Ecal) {
        ++nEM;
      }
    }

    double E_em = tower->emEnergy();
    if( E_em != 0 ) Ecal_subtowers.push_back( subtower( E_em, nEM ) );
      
    double E_HO = tower->outerEnergy();
    if( E_HO != 0 ) HO_subtowers.push_back( subtower( E_HO, nHO ) );
      
    double E_had = tower->hadEnergy();
    if( E_had != 0 ) {
      Hcal_subtowers.push_back( subtower( E_had, nHad ) );
      // totHcalE += E_had;
      
      int iEta = tower->ieta();
      Region reg = region( iEta );
      int iPhi = tower->iphi();
      if( iDbg>3 ) cout<<"tower has E_had: "<<E_had<<" iEta: "<<iEta
                       <<", iPhi: "<<iPhi<<" -> "<<reg;
      
      if( reg == HEneg || reg == HBneg || reg == HBpos || reg == HEpos ) {
        int oddnessIEta = HBHE_oddness( iEta );
        if( oddnessIEta < 0 ) break; // can't assign this tower to a single readout component
        
        int iHPD = 100 * reg;
        int iRBX = 100 * reg + ((iPhi + 1) % 72) / 4; // 71,72,1,2 are in the same RBX module
        
        if(( reg == HEneg || reg == HEpos ) && std::abs( iEta ) >= 21 ) { // at low-granularity edge of HE
          if( (0x1 & iPhi) == 0 ) {
            edm::LogError("CodeAssumptionsViolated")<<
              "Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
            return;
          }
          bool boolOddnessIEta = oddnessIEta;
          bool upperIPhi = (( iPhi%4 ) == 1 || ( iPhi%4 ) == 2); // the upper iPhi indices in the HE wedge
          // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
          // 1) in the upper iPhis of the module, the even IEtas belong to the higher iPhi
          // 2) in the loewr iPhis of the module, the odd  IEtas belong to the higher iPhi
          if( upperIPhi != boolOddnessIEta ) ++iPhi; 
          // note that iPhi could not be 72 before, so it's still in the legal range [1,72]
        } // if at low-granularity edge of HE
        iHPD += iPhi;
        
        // book the energies
        HPD_energy_map[ iHPD ] += E_had;
        RBX_energy_map[ iRBX ] += E_had;
        if( iDbg > 5 ) cout<<" --> H["<<iHPD<<"]="<<HPD_energy_map[iHPD]
                           <<", R["<<iRBX<<"]="<<RBX_energy_map[iRBX];
      } // HBHE
    } // E_had > 0
  } // loop on towers

  // sort the subtowers
  std::sort( Hcal_subtowers.begin(), Hcal_subtowers.end(), subtower_has_greater_E );
  std::sort( Ecal_subtowers.begin(), Ecal_subtowers.end(), subtower_has_greater_E );

  // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
  std::transform( HPD_energy_map.begin(), HPD_energy_map.end(), 
                  std::inserter (HPD_energies, HPD_energies.end()), select2nd ); 
  //              std::select2nd<std::map<int,double>::value_type>());
  std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
  std::transform( RBX_energy_map.begin(), RBX_energy_map.end(), 
                  std::inserter (RBX_energies, RBX_energies.end()), select2nd );
  //              std::select2nd<std::map<int,double>::value_type>());
  std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );

  subtowers.insert( subtowers.end(), Hcal_subtowers.begin(), Hcal_subtowers.end() );
  subtowers.insert( subtowers.end(), Ecal_subtowers.begin(), Ecal_subtowers.end() );
  subtowers.insert( subtowers.end(), HO_subtowers.begin(), HO_subtowers.end() );
  std::sort( subtowers.begin(), subtowers.end(), subtower_has_greater_E );
}

double reco::helper::JetIDHelper::fEB ( ) const [inline]

Definition at line 40 of file JetIDHelper.h.

References fEB_.

Referenced by JetIDProducer::produce().

{ return    fEB_;}

double reco::helper::JetIDHelper::fEE ( ) const [inline]

Definition at line 41 of file JetIDHelper.h.

References fEE_.

Referenced by JetIDProducer::produce().

{ return    fEE_;}

double reco::helper::JetIDHelper::fHB ( ) const [inline]

Definition at line 42 of file JetIDHelper.h.

References fHB_.

Referenced by JetIDProducer::produce().

{ return    fHB_;}

double reco::helper::JetIDHelper::fHE ( ) const [inline]

Definition at line 43 of file JetIDHelper.h.

References fHE_.

Referenced by JetIDProducer::produce().

{ return    fHE_;}

double reco::helper::JetIDHelper::fHFOOT ( ) const [inline]

Definition at line 48 of file JetIDHelper.h.

References fHFOOT_.

Referenced by JetIDProducer::produce().

{ return    fHFOOT_;}

double reco::helper::JetIDHelper::fHO ( ) const [inline]

Definition at line 44 of file JetIDHelper.h.

References fHO_.

Referenced by JetIDProducer::produce().

{ return    fHO_;}

double reco::helper::JetIDHelper::fHPD ( ) const [inline]

Definition at line 33 of file JetIDHelper.h.

References fHPD_.

Referenced by JetMETHLTOfflineSource::analyze(), JetIDProducer::produce(), BTagHLTOfflineSource::selectJets(), and FourVectorHLTOffline::selectJets().

{ return    fHPD_;}

void reco::helper::JetIDHelper::fillDescription ( edm::ParameterSetDescription & iDesc )

Definition at line 72 of file JetIDHelper.cc.

References edm::ParameterSetDescription::ifValue().

{
  iDesc.ifValue( edm::ParameterDescription<bool>("useRecHits", true, true),
                 true >> (edm::ParameterDescription<edm::InputTag>("hbheRecHitsColl", edm::InputTag(), true) and
                          edm::ParameterDescription<edm::InputTag>("hoRecHitsColl", edm::InputTag(), true) and
                          edm::ParameterDescription<edm::InputTag>("hfRecHitsColl", edm::InputTag(), true) and
                          edm::ParameterDescription<edm::InputTag>("ebRecHitsColl", edm::InputTag(), true) and
                          edm::ParameterDescription<edm::InputTag>("eeRecHitsColl", edm::InputTag(), true)
                          ) 
                 )->setComment("If using RecHits to calculate the precise jet ID variables that need them, "
                               "their sources need to be specified");
}

double reco::helper::JetIDHelper::fLong ( ) const [inline]

Definition at line 45 of file JetIDHelper.h.

References fLong_.

Referenced by JetIDProducer::produce().

{ return    fLong_;}

double reco::helper::JetIDHelper::fLSbad ( ) const [inline]

Definition at line 47 of file JetIDHelper.h.

References fLS_.

Referenced by JetIDProducer::produce().

{ return    fLS_;}

double reco::helper::JetIDHelper::fRBX ( ) const [inline]

Definition at line 34 of file JetIDHelper.h.

References fRBX_.

Referenced by JetIDProducer::produce().

{ return    fRBX_;}

double reco::helper::JetIDHelper::fShort ( ) const [inline]

Definition at line 46 of file JetIDHelper.h.

References fShort_.

Referenced by JetIDProducer::produce().

{ return    fShort_;}

double reco::helper::JetIDHelper::fSubDetector1 ( ) const [inline]

Definition at line 36 of file JetIDHelper.h.

References fSubDetector1_.

Referenced by JetIDProducer::produce().

{ return    fSubDetector1_;}

double reco::helper::JetIDHelper::fSubDetector2 ( ) const [inline]

Definition at line 37 of file JetIDHelper.h.

References fSubDetector2_.

Referenced by JetIDProducer::produce().

{ return    fSubDetector2_;}

double reco::helper::JetIDHelper::fSubDetector3 ( ) const [inline]

Definition at line 38 of file JetIDHelper.h.

References fSubDetector3_.

Referenced by JetIDProducer::produce().

{ return    fSubDetector3_;}

double reco::helper::JetIDHelper::fSubDetector4 ( ) const [inline]

Definition at line 39 of file JetIDHelper.h.

References fSubDetector4_.

Referenced by JetIDProducer::produce().

{ return    fSubDetector4_;}

int reco::helper::JetIDHelper::HBHE_oddness ( int iEta ) [private]

Definition at line 588 of file JetIDHelper.cc.

{
 int ae = TMath::Abs( iEta );
 if( ae == 29 ) return -1; // can't figure it out without RecHits
 return ae & 0x1;
}

int reco::helper::JetIDHelper::HBHE_oddness	(	int	iEta,
		int	depth
	)		`[private]`

Definition at line 572 of file JetIDHelper.cc.

{
 int ae = TMath::Abs( iEta );
 if( ae == 29 && depth == 1 ) ae += 1; // observed that: hits are at depths 1 & 2; 1 goes with the even pattern
 return ae & 0x1;
}

reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::HBHE_region	(	int	iEta,
		int	depth
	)		`[private]`

Definition at line 579 of file JetIDHelper.cc.

{
  // no error checking for HO indices (depth 2 & |ieta|<=14 or depth 3 & |ieta|=15)
  if( iEta <= -17 || ( depth == 3 && iEta == -16 ) ) return HEneg;
  if( iEta >=  17 || ( depth == 3 && iEta ==  16 ) ) return HEpos;
  if( iEta < 0 ) return HBneg;
  return HBpos;
}

int reco::helper::JetIDHelper::hitsInN90 ( ) const [inline]

Definition at line 56 of file JetIDHelper.h.

References hitsInN90_.

Referenced by JetIDProducer::produce().

{ return    hitsInN90_;}

unsigned int reco::helper::JetIDHelper::hitsInNCarrying	(	double	fraction,
		std::vector< subtower >	descending_towers
	)		`[private]`

Definition at line 209 of file JetIDHelper.cc.

References i.

{
  double totalE = 0;
  for( unsigned int i = 0; i < descending_towers.size(); ++i ) totalE += descending_towers[ i ].E;

  double runningE = 0;
  unsigned int NH = 0;
  
  // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
  for( unsigned int i = descending_towers.size(); i > 0; --i ) {
    runningE += descending_towers[ i-1 ].E;
    if (runningE >= ( 1-fraction ) * totalE) NH += descending_towers[ i-1 ].Nhit;
  }
  return NH;
}

void reco::helper::JetIDHelper::initValues ( )

Definition at line 53 of file JetIDHelper.cc.

{
  fHPD_= -1.0;
  fRBX_= -1.0;
  n90Hits_ = -1;
  fSubDetector1_= -1.0;
  fSubDetector2_= -1.0;
  fSubDetector3_= -1.0;
  fSubDetector4_= -1.0;
  restrictedEMF_= -1.0;
  nHCALTowers_ = -1;
  nECALTowers_ = -1;
  approximatefHPD_ = -1.0;
  approximatefRBX_ = -1.0;
  hitsInN90_ = -1;
  fEB_ = fEE_ = fHB_ = fHE_ = fHO_ = fLong_ = fShort_ = -1.0;
  fLS_ = fHFOOT_ = -1.0;
}

int reco::helper::JetIDHelper::n90Hits ( ) const [inline]

Definition at line 35 of file JetIDHelper.h.

References n90Hits_.

Referenced by HLTJets::analyze(), JetMETHLTOfflineSource::analyze(), HLTCaloJetIDProducer::produce(), and JetIDProducer::produce().

{ return    n90Hits_;}

unsigned int reco::helper::JetIDHelper::nCarrying	(	double	fraction,
		std::vector< double >	descending_energies
	)		`[private]`

Definition at line 192 of file JetIDHelper.cc.

References i.

{
  double totalE = 0;
  for( unsigned int i = 0; i < descending_energies.size(); ++i ) totalE += descending_energies[ i ];

  double runningE = 0;
  unsigned int NC = descending_energies.size();
  
  // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
  for( unsigned int i = descending_energies.size(); i > 0; --i ) {
    runningE += descending_energies[ i-1 ];
    if (runningE < ( 1-fraction ) * totalE) NC = i-1;
  }
  return NC;
}

int reco::helper::JetIDHelper::nECALTowers ( ) const [inline]

Definition at line 52 of file JetIDHelper.h.

References nECALTowers_.

Referenced by JetIDProducer::produce().

{ return    nECALTowers_;}

int reco::helper::JetIDHelper::nHCALTowers ( ) const [inline]

Definition at line 51 of file JetIDHelper.h.

References nHCALTowers_.

Referenced by JetIDProducer::produce().

{ return    nHCALTowers_;}

reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::region ( int iEta ) [private]

Definition at line 595 of file JetIDHelper.cc.

{
  if( iEta == 16 || iEta == -16 ) return unknown_region; // both HB and HE cells belong to these towers
  if( iEta == 29 || iEta == -29 ) return unknown_region; // both HE and HF cells belong to these towers
  if( iEta <= -30 ) return HFneg;
  if( iEta >=  30 ) return HFpos;
  if( iEta <= -17 ) return HEneg;
  if( iEta >=  17 ) return HEpos;
  if( iEta < 0 ) return HBneg;
  return HBpos;
}