Reference Manual

Go to the documentation of this file.
00001 #include "RecoJets/JetProducers/interface/JetIDHelper.h"
00002 
00003 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00004 #include "DataFormats/HcalRecHit/interface/HcalRecHitFwd.h"
00005 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00006 #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
00007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00008 
00009 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"
00010 // JetIDHelper needs a much more detailed description that the one in HcalTopology, 
00011 // so to be consistent, all needed constants are hardwired in JetIDHelper.cc itself
00012 // #include "Geometry/CaloTopology/interface/HcalTopology.h"
00013 
00014 #include "TMath.h"
00015 #include <vector>
00016 #include <numeric>
00017 #include <iostream>
00018 
00019 using namespace std;
00020 
00021 // defining stuff useful only for this class (not needed in header)
00022 namespace reco {
00023 
00024   namespace helper {
00025 
00026     // select2nd exists only in some std and boost implementations, so let's control our own fate
00027     // and it can't be a non-static member function.
00028     static double select2nd (std::map<int,double>::value_type const &pair) {return pair.second;}
00029     
00030     bool hasNonPositiveE( reco::helper::JetIDHelper::subtower x ) {
00031       return x.E <= 0;
00032     }
00033     
00034     bool subtower_has_greater_E( reco::helper::JetIDHelper::subtower i, 
00035                                  reco::helper::JetIDHelper::subtower j ) { return i.E > j.E; }
00036     int JetIDHelper::sanity_checks_left_ = 100;
00037   }
00038 }
00039 
00040 reco::helper::JetIDHelper::JetIDHelper( edm::ParameterSet const & pset )
00041 {
00042   useRecHits_ = pset.getParameter<bool>("useRecHits");
00043   if( useRecHits_ ) {
00044     hbheRecHitsColl_ = pset.getParameter<edm::InputTag>("hbheRecHitsColl");
00045     hoRecHitsColl_   = pset.getParameter<edm::InputTag>("hoRecHitsColl");
00046     hfRecHitsColl_   = pset.getParameter<edm::InputTag>("hfRecHitsColl");
00047     ebRecHitsColl_   = pset.getParameter<edm::InputTag>("ebRecHitsColl");
00048     eeRecHitsColl_   = pset.getParameter<edm::InputTag>("eeRecHitsColl");   
00049   }
00050   initValues();
00051 }
00052   
00053 void reco::helper::JetIDHelper::initValues()
00054 {
00055   fHPD_= -1.0;
00056   fRBX_= -1.0;
00057   n90Hits_ = -1;
00058   fSubDetector1_= -1.0;
00059   fSubDetector2_= -1.0;
00060   fSubDetector3_= -1.0;
00061   fSubDetector4_= -1.0;
00062   restrictedEMF_= -1.0;
00063   nHCALTowers_ = -1;
00064   nECALTowers_ = -1;
00065   approximatefHPD_ = -1.0;
00066   approximatefRBX_ = -1.0;
00067   hitsInN90_ = -1;
00068   fEB_ = fEE_ = fHB_ = fHE_ = fHO_ = fLong_ = fShort_ = -1.0;
00069   fLS_ = fHFOOT_ = -1.0;
00070 }
00071 
00072 void reco::helper::JetIDHelper::fillDescription(edm::ParameterSetDescription& iDesc)
00073 {
00074   iDesc.ifValue( edm::ParameterDescription<bool>("useRecHits", true, true),
00075                  true >> (edm::ParameterDescription<edm::InputTag>("hbheRecHitsColl", edm::InputTag(), true) and
00076                           edm::ParameterDescription<edm::InputTag>("hoRecHitsColl", edm::InputTag(), true) and
00077                           edm::ParameterDescription<edm::InputTag>("hfRecHitsColl", edm::InputTag(), true) and
00078                           edm::ParameterDescription<edm::InputTag>("ebRecHitsColl", edm::InputTag(), true) and
00079                           edm::ParameterDescription<edm::InputTag>("eeRecHitsColl", edm::InputTag(), true)
00080                           ) 
00081                  )->setComment("If using RecHits to calculate the precise jet ID variables that need them, "
00082                                "their sources need to be specified");
00083 }
00084 
00085 
00086 void reco::helper::JetIDHelper::calculate( const edm::Event& event, const reco::CaloJet &jet, const int iDbg )
00087 {
00088   initValues();
00089 
00090   // --------------------------------------------------
00091   // 1) jet ID variable derived from existing fractions
00092   // --------------------------------------------------
00093 
00094   double E_EM = TMath::Max( float(0.), jet.emEnergyInHF() )
00095               + TMath::Max( float(0.), jet.emEnergyInEB() )
00096               + TMath::Max( float(0.), jet.emEnergyInEE() );
00097   double E_Had = TMath::Max( float(0.), jet.hadEnergyInHB() )
00098                + TMath::Max( float(0.), jet.hadEnergyInHE() )
00099                + TMath::Max( float(0.), jet.hadEnergyInHO() )
00100                + TMath::Max( float(0.), jet.hadEnergyInHF() );
00101   if( E_Had + E_EM > 0 ) restrictedEMF_ = E_EM / ( E_EM + E_Had );
00102   if( iDbg > 1 ) cout<<"jet pT: "<<jet.pt()<<", eT: "<<jet.et()<<", E: "
00103                      <<jet.energy()<<" rEMF: "<<restrictedEMF_<<endl;
00104 
00105   // ------------------------
00106   // 2) tower based variables
00107   // ------------------------
00108   vector<subtower> subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers;
00109   vector<double> HPD_energies, RBX_energies;
00110 
00111   classifyJetTowers( event, jet, 
00112                      subtowers, Ecal_subtowers, Hcal_subtowers, HO_subtowers, 
00113                      HPD_energies, RBX_energies, iDbg );
00114   if( iDbg > 1 ) {
00115     cout<<"E:";
00116     for (unsigned int i=0; i<subtowers.size(); ++i) cout<<" "<<subtowers[i].E<<","<<subtowers[i].Nhit;
00117     cout<<"\nECal_E:";
00118     for (unsigned int i=0; i<Ecal_subtowers.size(); ++i) cout<<" "<<Ecal_subtowers[i].E<<","<<Ecal_subtowers[i].Nhit;
00119     cout<<"\nHCal_E:";
00120     for (unsigned int i=0; i<Hcal_subtowers.size(); ++i) cout<<" "<<Hcal_subtowers[i].E<<","<<Hcal_subtowers[i].Nhit;
00121     cout<<"\nHO_E:";
00122     for (unsigned int i=0; i<HO_subtowers.size(); ++i) cout<<" "<<HO_subtowers[i].E<<","<<HO_subtowers[i].Nhit;
00123     cout<<"\nHPD_E:";
00124     for (unsigned int i=0; i<HPD_energies.size(); ++i) cout<<" "<<HPD_energies[i];
00125     cout<<"\nRBX_E:";
00126     for (unsigned int i=0; i<RBX_energies.size(); ++i) cout<<" "<<RBX_energies[i];
00127     cout<<endl;
00128   }
00129 
00130   // counts
00131   hitsInN90_ = hitsInNCarrying( 0.9, subtowers );
00132   nHCALTowers_ = Hcal_subtowers.size();
00133   vector<subtower>::const_iterator it;
00134   it = find_if( Ecal_subtowers.begin(), Ecal_subtowers.end(), hasNonPositiveE );
00135   nECALTowers_ = it - Ecal_subtowers.begin(); // ignores negative energies from HF!
00136 
00137   // energy fractions
00138   if( jet.energy() > 0 ) {
00139     if( HPD_energies.size() > 0 ) approximatefHPD_ = HPD_energies.at( 0 ) / jet.energy();
00140     if( RBX_energies.size() > 0 ) approximatefRBX_ = RBX_energies.at( 0 ) / jet.energy();
00141   }
00142 
00143   // -----------------------
00144   // 3) cell based variables
00145   // -----------------------
00146   if( useRecHits_ ) {
00147     vector<double> energies, subdet_energies, Ecal_energies, Hcal_energies, HO_energies;
00148     double LS_bad_energy, HF_OOT_energy;
00149     classifyJetComponents( event, jet, 
00150                            energies, subdet_energies, Ecal_energies, Hcal_energies, HO_energies, 
00151                            HPD_energies, RBX_energies, LS_bad_energy, HF_OOT_energy, iDbg );
00152 
00153     // counts
00154     n90Hits_ = nCarrying( 0.9, energies );
00155 
00156     // energy fractions
00157     if( jet.energy() > 0 ) {
00158       if( HPD_energies.size() > 0 ) fHPD_ = HPD_energies.at( 0 ) / jet.energy();
00159       if( RBX_energies.size() > 0 ) fRBX_ = RBX_energies.at( 0 ) / jet.energy();
00160       if( subdet_energies.size() > 0 ) fSubDetector1_ = subdet_energies.at( 0 ) / jet.energy();
00161       if( subdet_energies.size() > 1 ) fSubDetector2_ = subdet_energies.at( 1 ) / jet.energy();
00162       if( subdet_energies.size() > 2 ) fSubDetector3_ = subdet_energies.at( 2 ) / jet.energy();
00163       if( subdet_energies.size() > 3 ) fSubDetector4_ = subdet_energies.at( 3 ) / jet.energy();
00164       fLS_ = LS_bad_energy / jet.energy();
00165       fHFOOT_ = HF_OOT_energy / jet.energy();
00166 
00167       if( sanity_checks_left_ > 0 ) {
00168         --sanity_checks_left_;
00169         double EH_sum = accumulate( Ecal_energies.begin(), Ecal_energies.end(), 0. );
00170         EH_sum = accumulate( Hcal_energies.begin(), Hcal_energies.end(), EH_sum );
00171         double EHO_sum = accumulate( HO_energies.begin(), HO_energies.end(), EH_sum );
00172         if( jet.energy() > 0.001 && abs( EH_sum/jet.energy() - 1 ) > 0.01 && abs( EHO_sum/jet.energy() - 1 ) > 0.01 )
00173           edm::LogWarning("BadInput")<<"Jet energy ("<<jet.energy()
00174                                      <<") does not match the total energy in the recHits ("<<EHO_sum
00175                                      <<", or without HO: "<<EH_sum<<") . Are these the right recHits? "
00176                                      <<"Were jet energy corrections mistakenly applied before jet ID? A bug?";
00177         if( iDbg > 1 ) cout<<"Sanity check - E: "<<jet.energy()<<" =? EH_sum: "<<EH_sum<<" / EHO_sum: "<<EHO_sum<<endl;
00178       }
00179     }
00180 
00181     if( iDbg > 1 ) {
00182       cout<<"DBG - fHPD: "<<fHPD_<<", fRBX: "<<fRBX_<<", nh90: "<<n90Hits_<<", fLS: "<<fLS_<<", fHFOOT: "<<fHFOOT_<<endl;
00183       cout<<"    -~fHPD: "<<approximatefHPD_<<", ~fRBX: "<<approximatefRBX_
00184           <<", hits in n90: "<<hitsInN90_<<endl;
00185       cout<<"    - nHCALTowers: "<<nHCALTowers_<<", nECALTowers: "<<nECALTowers_
00186           <<"; subdet fractions: "<<fSubDetector1_<<", "<<fSubDetector2_<<", "<<fSubDetector3_<<", "<<fSubDetector4_<<endl;
00187     }
00188   }
00189 }
00190 
00191 
00192 unsigned int reco::helper::JetIDHelper::nCarrying( double fraction, vector< double > descending_energies )
00193 {
00194   double totalE = 0;
00195   for( unsigned int i = 0; i < descending_energies.size(); ++i ) totalE += descending_energies[ i ];
00196 
00197   double runningE = 0;
00198   unsigned int NC = descending_energies.size();
00199   
00200   // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
00201   for( unsigned int i = descending_energies.size(); i > 0; --i ) {
00202     runningE += descending_energies[ i-1 ];
00203     if (runningE < ( 1-fraction ) * totalE) NC = i-1;
00204   }
00205   return NC;
00206 }
00207 
00208 
00209 unsigned int reco::helper::JetIDHelper::hitsInNCarrying( double fraction, vector< subtower > descending_towers )
00210 {
00211   double totalE = 0;
00212   for( unsigned int i = 0; i < descending_towers.size(); ++i ) totalE += descending_towers[ i ].E;
00213 
00214   double runningE = 0;
00215   unsigned int NH = 0;
00216   
00217   // slightly odd loop structure avoids round-off problems when runningE never catches up with totalE
00218   for( unsigned int i = descending_towers.size(); i > 0; --i ) {
00219     runningE += descending_towers[ i-1 ].E;
00220     if (runningE >= ( 1-fraction ) * totalE) NH += descending_towers[ i-1 ].Nhit;
00221   }
00222   return NH;
00223 }
00224 
00225 void reco::helper::JetIDHelper::classifyJetComponents( const edm::Event& event, const reco::CaloJet &jet, 
00226                                                        vector< double > &energies,      
00227                                                        vector< double > &subdet_energies,      
00228                                                        vector< double > &Ecal_energies, 
00229                                                        vector< double > &Hcal_energies, 
00230                                                        vector< double > &HO_energies,
00231                                                        vector< double > &HPD_energies,  
00232                                                        vector< double > &RBX_energies,
00233                                                        double& LS_bad_energy, double& HF_OOT_energy, const int iDbg )
00234 {
00235   energies.clear(); subdet_energies.clear(); Ecal_energies.clear(); Hcal_energies.clear(); HO_energies.clear();
00236   HPD_energies.clear(); RBX_energies.clear();
00237   LS_bad_energy = HF_OOT_energy = 0.;
00238   
00239   std::map< int, double > HPD_energy_map, RBX_energy_map;
00240   vector< double > EB_energies, EE_energies, HB_energies, HE_energies, short_energies, long_energies;
00241   edm::Handle<HBHERecHitCollection> HBHERecHits;
00242   edm::Handle<HORecHitCollection> HORecHits;
00243   edm::Handle<HFRecHitCollection> HFRecHits;
00244   edm::Handle<EBRecHitCollection> EBRecHits;
00245   edm::Handle<EERecHitCollection> EERecHits;
00246   // the jet only contains DetIds, so first read recHit collection
00247   event.getByLabel( hbheRecHitsColl_, HBHERecHits );
00248   event.getByLabel( hoRecHitsColl_, HORecHits );
00249   event.getByLabel( hfRecHitsColl_, HFRecHits );
00250   event.getByLabel( ebRecHitsColl_, EBRecHits );
00251   event.getByLabel( eeRecHitsColl_, EERecHits );
00252   if( iDbg > 2 ) cout<<"# of rechits found - HBHE: "<<HBHERecHits->size()
00253                      <<", HO: "<<HORecHits->size()<<", HF: "<<HFRecHits->size()
00254                      <<", EB: "<<EBRecHits->size()<<", EE: "<<EERecHits->size()<<endl;
00255 
00256   vector< CaloTowerPtr > towers = jet.getCaloConstituents ();
00257   int nTowers = towers.size();
00258   if( iDbg > 9 ) cout<<"In classifyJetComponents. # of towers found: "<<nTowers<<endl;
00259 
00260   for( int iTower = 0; iTower <nTowers ; iTower++ ) {
00261 
00262     CaloTowerPtr& tower = towers[iTower];
00263 
00264     int nCells = tower->constituentsSize();
00265     if( iDbg ) cout<<"tower #"<<iTower<<" has "<<nCells<<" cells. "
00266                    <<"It's at iEta: "<<tower->ieta()<<", iPhi: "<<tower->iphi()<<endl;
00267 
00268     const vector< DetId >& cellIDs = tower->constituents();  // cell == recHit
00269   
00270     for( int iCell = 0; iCell < nCells; ++iCell ) {
00271       DetId::Detector detNum = cellIDs[iCell].det();
00272       if( detNum == DetId::Hcal ) {
00273         HcalDetId HcalID = cellIDs[ iCell ];
00274         HcalSubdetector HcalNum = HcalID.subdet();
00275         double hitE = 0;
00276         if( HcalNum == HcalOuter ) {
00277           HORecHitCollection::const_iterator theRecHit=HORecHits->find(HcalID);
00278           if (theRecHit == HORecHits->end()) {
00279             edm::LogWarning("UnexpectedEventContents")<<"Can't find the HO recHit with ID: "<<HcalID;
00280             continue;
00281           }
00282           hitE = theRecHit->energy();
00283           HO_energies.push_back( hitE );
00284             
00285         } else if( HcalNum == HcalForward ) {
00286           
00287           HFRecHitCollection::const_iterator theRecHit=HFRecHits->find( HcalID );           
00288           if( theRecHit == HFRecHits->end() ) {
00289             edm::LogWarning("UnexpectedEventContents")<<"Can't find the HF recHit with ID: "<<HcalID;
00290             continue;
00291           }
00292           hitE = theRecHit->energy();
00293           if( iDbg>4 ) cout 
00294                          << "hit #"<<iCell<<" is  HF , E: "<<hitE<<" iEta: "<<theRecHit->id().ieta()
00295                          <<", depth: "<<theRecHit->id().depth()<<", iPhi: "
00296                          <<theRecHit->id().iphi();
00297           
00298           if( HcalID.depth() == 1 ) 
00299             long_energies.push_back( hitE );
00300           else
00301             short_energies.push_back( hitE );
00302           
00303           uint32_t flags = theRecHit->flags();
00304           if( flags & (1<<HcalCaloFlagLabels::HFLongShort) ) LS_bad_energy += hitE;
00305           if( flags & ( (1<<HcalCaloFlagLabels::HFDigiTime) | 
00306                         (3<<HcalCaloFlagLabels::HFTimingTrustBits) | 
00307                         (1<<HcalCaloFlagLabels::TimingSubtractedBit) |
00308                         (1<<HcalCaloFlagLabels::TimingAddedBit) |
00309                         (1<<HcalCaloFlagLabels::TimingErrorBit) ) ) HF_OOT_energy += hitE;
00310           if( iDbg>4 && flags ) cout<<"flags: "<<flags
00311                                     <<" -> LS_bad_energy: "<<LS_bad_energy<<", HF_OOT_energy: "<<HF_OOT_energy<<endl; 
00312 
00313         } else { // HBHE
00314           
00315           HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find( HcalID );     
00316           if( theRecHit == HBHERecHits->end() ) {
00317             edm::LogWarning("UnexpectedEventContents")<<"Can't find the HBHE recHit with ID: "<<HcalID; 
00318             continue;
00319           }
00320           hitE = theRecHit->energy();
00321           int iEta = theRecHit->id().ieta();
00322           int depth = theRecHit->id().depth();
00323           Region region = HBHE_region( iEta, depth );
00324           int hitIPhi = theRecHit->id().iphi();
00325           if( iDbg>3 ) cout<<"hit #"<<iCell<<" is HBHE, E: "<<hitE<<" iEta: "<<iEta
00326                            <<", depth: "<<depth<<", iPhi: "<<theRecHit->id().iphi()
00327                            <<" -> "<<region;
00328           int absIEta = TMath::Abs( theRecHit->id().ieta() );
00329           if( depth == 3 && (absIEta == 28 || absIEta == 29) ) {
00330             hitE /= 2; // Depth 3 at the HE forward edge is split over tower 28 & 29, and jet reco. assigns half each
00331           }
00332           int iHPD = 100 * region;
00333           int iRBX = 100 * region + ((hitIPhi + 1) % 72) / 4; // 71,72,1,2 are in the same RBX module
00334           
00335           if( std::abs( iEta ) >= 21 ) {
00336             if( (0x1 & hitIPhi) == 0 ) {
00337               edm::LogError("CodeAssumptionsViolated")<<"Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
00338               return;
00339             }
00340             bool oddnessIEta = HBHE_oddness( iEta, depth );
00341             bool upperIPhi = (( hitIPhi%4 ) == 1 || ( hitIPhi%4 ) == 2); // the upper iPhi indices in the HE wedge
00342             // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
00343             // 1) in the upper iPhis of the module, the even iEtas belong to the higher iPhi
00344             // 2) in the loewr iPhis of the module, the odd  iEtas belong to the higher iPhi
00345             if( upperIPhi != oddnessIEta ) ++hitIPhi; 
00346             // note that hitIPhi could not be 72 before, so it's still in the legal range [1,72]
00347           }
00348           iHPD += hitIPhi;
00349           
00350           // book the energies
00351           HPD_energy_map[ iHPD ] += hitE;
00352           RBX_energy_map[ iRBX ] += hitE;
00353           if( iDbg > 5 ) cout<<" --> H["<<iHPD<<"]="<<HPD_energy_map[iHPD]
00354                              <<", R["<<iRBX<<"]="<<RBX_energy_map[iRBX];
00355           if( iDbg > 1 ) cout<<endl;
00356           
00357           if( region == HBneg || region == HBpos )
00358             HB_energies.push_back( hitE );
00359           else
00360             HE_energies.push_back( hitE );
00361           
00362         } // if HBHE
00363         if( hitE == 0 ) edm::LogWarning("UnexpectedEventContents")<<"HCal hitE==0? (or unknown subdetector?)";
00364         
00365       } // if HCAL 
00366         
00367       else if (detNum == DetId::Ecal) {
00368         
00369         int EcalNum =  cellIDs[iCell].subdetId();
00370         double hitE = 0;
00371         if( EcalNum == 1 ){
00372           EBDetId EcalID = cellIDs[iCell];
00373           EBRecHitCollection::const_iterator theRecHit=EBRecHits->find(EcalID);
00374           if( theRecHit == EBRecHits->end() ) {edm::LogWarning("UnexpectedEventContents")
00375               <<"Can't find the EB recHit with ID: "<<EcalID; continue;}
00376           hitE = theRecHit->energy();
00377           EB_energies.push_back( hitE );
00378         } else if(  EcalNum == 2 ){
00379           EEDetId EcalID = cellIDs[iCell];
00380           EERecHitCollection::const_iterator theRecHit=EERecHits->find(EcalID);
00381           if( theRecHit == EERecHits->end() ) {edm::LogWarning("UnexpectedEventContents")
00382               <<"Can't find the EE recHit with ID: "<<EcalID; continue;}
00383           hitE = theRecHit->energy();
00384           EE_energies.push_back( hitE );
00385         }
00386         if( hitE == 0 ) edm::LogWarning("UnexpectedEventContents")<<"ECal hitE==0? (or unknown subdetector?)";
00387         if( iDbg > 6 ) cout<<"EcalNum: "<<EcalNum<<" hitE: "<<hitE<<endl;
00388       } // 
00389     } // loop on cells
00390   } // loop on towers
00391 
00392   /* Disabling check until HO is accounted for in EMF. Check was used in CMSSW_2, where HE was excluded.
00393   double expHcalE = jet.energy() * (1-jet.emEnergyFraction());
00394   if( totHcalE + expHcalE > 0 && 
00395       TMath::Abs( totHcalE - expHcalE ) > 0.01 && 
00396       ( totHcalE - expHcalE ) / ( totHcalE + expHcalE ) > 0.0001 ) {
00397     edm::LogWarning("CodeAssumptionsViolated")<<"failed to account for all Hcal energies"
00398                                               <<totHcalE<<"!="<<expHcalE;
00399                                               } */   
00400   // concatenate Hcal and Ecal lists
00401   Hcal_energies.insert( Hcal_energies.end(), HB_energies.begin(), HB_energies.end() );
00402   Hcal_energies.insert( Hcal_energies.end(), HE_energies.begin(), HE_energies.end() );
00403   Hcal_energies.insert( Hcal_energies.end(), short_energies.begin(), short_energies.end() );
00404   Hcal_energies.insert( Hcal_energies.end(), long_energies.begin(), long_energies.end() );
00405   Ecal_energies.insert( Ecal_energies.end(), EB_energies.begin(), EB_energies.end() );
00406   Ecal_energies.insert( Ecal_energies.end(), EE_energies.begin(), EE_energies.end() );
00407 
00408   // sort the energies
00409   std::sort( Hcal_energies.begin(), Hcal_energies.end(), greater<double>() );
00410   std::sort( Ecal_energies.begin(), Ecal_energies.end(), greater<double>() );
00411 
00412   // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
00413   std::transform( HPD_energy_map.begin(), HPD_energy_map.end(), 
00414                   std::inserter (HPD_energies, HPD_energies.end()), select2nd ); 
00415   //              std::select2nd<std::map<int,double>::value_type>());
00416   std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
00417   std::transform( RBX_energy_map.begin(), RBX_energy_map.end(), 
00418                   std::inserter (RBX_energies, RBX_energies.end()), select2nd );
00419   //              std::select2nd<std::map<int,double>::value_type>());
00420   std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );
00421 
00422   energies.insert( energies.end(), Hcal_energies.begin(), Hcal_energies.end() );
00423   energies.insert( energies.end(), Ecal_energies.begin(), Ecal_energies.end() );
00424   energies.insert( energies.end(), HO_energies.begin(), HO_energies.end() );
00425   std::sort( energies.begin(), energies.end(), greater<double>() );
00426 
00427   // prepare sub detector energies, then turn them into fractions
00428   fEB_ = std::accumulate( EB_energies.begin(), EB_energies.end(), 0. );
00429   fEE_ = std::accumulate( EE_energies.begin(), EE_energies.end(), 0. );
00430   fHB_ = std::accumulate( HB_energies.begin(), HB_energies.end(), 0. );
00431   fHE_ = std::accumulate( HE_energies.begin(), HE_energies.end(), 0. );
00432   fHO_ = std::accumulate( HO_energies.begin(), HO_energies.end(), 0. );
00433   fShort_ = std::accumulate( short_energies.begin(), short_energies.end(), 0. );
00434   fLong_ = std::accumulate( long_energies.begin(), long_energies.end(), 0. );
00435   subdet_energies.push_back( fEB_ );
00436   subdet_energies.push_back( fEE_ );
00437   subdet_energies.push_back( fHB_ );
00438   subdet_energies.push_back( fHE_ );
00439   subdet_energies.push_back( fHO_ );
00440   subdet_energies.push_back( fShort_ );
00441   subdet_energies.push_back( fLong_ );
00442   std::sort( subdet_energies.begin(), subdet_energies.end(), greater<double>() );
00443   if( jet.energy() > 0 ) {
00444     fEB_ /= jet.energy();
00445     fEE_ /= jet.energy();
00446     fHB_ /= jet.energy();
00447     fHE_ /= jet.energy();
00448     fHO_ /= jet.energy();
00449     fShort_ /= jet.energy();
00450     fLong_ /= jet.energy();
00451   } else {
00452     if( fEB_ > 0 || fEE_ > 0 || fHB_ > 0 || fHE_ > 0 || fHO_ > 0 || fShort_ > 0 || fLong_ > 0 )
00453       edm::LogError("UnexpectedEventContents")<<"Jet ID Helper found energy in subdetectors and jet E <= 0";
00454   }      
00455 }
00456 
00457 void reco::helper::JetIDHelper::classifyJetTowers( const edm::Event& event, const reco::CaloJet &jet, 
00458                                                    vector< subtower > &subtowers,      
00459                                                    vector< subtower > &Ecal_subtowers, 
00460                                                    vector< subtower > &Hcal_subtowers, 
00461                                                    vector< subtower > &HO_subtowers,
00462                                                    vector< double > &HPD_energies,  
00463                                                    vector< double > &RBX_energies,
00464                                                    const int iDbg )
00465 {
00466   subtowers.clear(); Ecal_subtowers.clear(); Hcal_subtowers.clear(); HO_subtowers.clear();
00467   HPD_energies.clear(); RBX_energies.clear();
00468 
00469   std::map< int, double > HPD_energy_map, RBX_energy_map;
00470 
00471   vector< CaloTowerPtr > towers = jet.getCaloConstituents ();
00472   int nTowers = towers.size();
00473   if( iDbg > 9 ) cout<<"classifyJetTowers started. # of towers found: "<<nTowers<<endl;
00474 
00475   for( int iTower = 0; iTower <nTowers ; iTower++ ) {
00476 
00477     CaloTowerPtr& tower = towers[iTower];
00478 
00479     int nEM = 0, nHad = 0, nHO = 0;
00480     const vector< DetId >& cellIDs = tower->constituents();  // cell == recHit
00481     int nCells = cellIDs.size();
00482     if( iDbg ) cout<<"tower #"<<iTower<<" has "<<nCells<<" cells. "
00483                    <<"It's at iEta: "<<tower->ieta()<<", iPhi: "<<tower->iphi()<<endl;
00484   
00485     for( int iCell = 0; iCell < nCells; ++iCell ) {
00486       DetId::Detector detNum = cellIDs[iCell].det();
00487       if( detNum == DetId::Hcal ) {
00488         HcalDetId HcalID = cellIDs[ iCell ];
00489         HcalSubdetector HcalNum = HcalID.subdet();
00490         if( HcalNum == HcalOuter ) {
00491           ++nHO;
00492         } else {
00493           ++nHad;
00494         }
00495       } else if (detNum == DetId::Ecal) {
00496         ++nEM;
00497       }
00498     }
00499 
00500     double E_em = tower->emEnergy();
00501     if( E_em != 0 ) Ecal_subtowers.push_back( subtower( E_em, nEM ) );
00502       
00503     double E_HO = tower->outerEnergy();
00504     if( E_HO != 0 ) HO_subtowers.push_back( subtower( E_HO, nHO ) );
00505       
00506     double E_had = tower->hadEnergy();
00507     if( E_had != 0 ) {
00508       Hcal_subtowers.push_back( subtower( E_had, nHad ) );
00509       // totHcalE += E_had;
00510       
00511       int iEta = tower->ieta();
00512       Region reg = region( iEta );
00513       int iPhi = tower->iphi();
00514       if( iDbg>3 ) cout<<"tower has E_had: "<<E_had<<" iEta: "<<iEta
00515                        <<", iPhi: "<<iPhi<<" -> "<<reg;
00516       
00517       if( reg == HEneg || reg == HBneg || reg == HBpos || reg == HEpos ) {
00518         int oddnessIEta = HBHE_oddness( iEta );
00519         if( oddnessIEta < 0 ) break; // can't assign this tower to a single readout component
00520         
00521         int iHPD = 100 * reg;
00522         int iRBX = 100 * reg + ((iPhi + 1) % 72) / 4; // 71,72,1,2 are in the same RBX module
00523         
00524         if(( reg == HEneg || reg == HEpos ) && std::abs( iEta ) >= 21 ) { // at low-granularity edge of HE
00525           if( (0x1 & iPhi) == 0 ) {
00526             edm::LogError("CodeAssumptionsViolated")<<
00527               "Bug?! Jet ID code assumes no even iPhi recHits at HE edges";
00528             return;
00529           }
00530           bool boolOddnessIEta = oddnessIEta;
00531           bool upperIPhi = (( iPhi%4 ) == 1 || ( iPhi%4 ) == 2); // the upper iPhi indices in the HE wedge
00532           // remap the iPhi so it fits the one in the inner HE regions, change in needed in two cases:
00533           // 1) in the upper iPhis of the module, the even IEtas belong to the higher iPhi
00534           // 2) in the loewr iPhis of the module, the odd  IEtas belong to the higher iPhi
00535           if( upperIPhi != boolOddnessIEta ) ++iPhi; 
00536           // note that iPhi could not be 72 before, so it's still in the legal range [1,72]
00537         } // if at low-granularity edge of HE
00538         iHPD += iPhi;
00539         
00540         // book the energies
00541         HPD_energy_map[ iHPD ] += E_had;
00542         RBX_energy_map[ iRBX ] += E_had;
00543         if( iDbg > 5 ) cout<<" --> H["<<iHPD<<"]="<<HPD_energy_map[iHPD]
00544                            <<", R["<<iRBX<<"]="<<RBX_energy_map[iRBX];
00545       } // HBHE
00546     } // E_had > 0
00547   } // loop on towers
00548 
00549   // sort the subtowers
00550   std::sort( Hcal_subtowers.begin(), Hcal_subtowers.end(), subtower_has_greater_E );
00551   std::sort( Ecal_subtowers.begin(), Ecal_subtowers.end(), subtower_has_greater_E );
00552 
00553   // put the energy sums (the 2nd entry in each pair of the maps) into the output vectors and sort them
00554   std::transform( HPD_energy_map.begin(), HPD_energy_map.end(), 
00555                   std::inserter (HPD_energies, HPD_energies.end()), select2nd ); 
00556   //              std::select2nd<std::map<int,double>::value_type>());
00557   std::sort( HPD_energies.begin(), HPD_energies.end(), greater<double>() );
00558   std::transform( RBX_energy_map.begin(), RBX_energy_map.end(), 
00559                   std::inserter (RBX_energies, RBX_energies.end()), select2nd );
00560   //              std::select2nd<std::map<int,double>::value_type>());
00561   std::sort( RBX_energies.begin(), RBX_energies.end(), greater<double>() );
00562 
00563   subtowers.insert( subtowers.end(), Hcal_subtowers.begin(), Hcal_subtowers.end() );
00564   subtowers.insert( subtowers.end(), Ecal_subtowers.begin(), Ecal_subtowers.end() );
00565   subtowers.insert( subtowers.end(), HO_subtowers.begin(), HO_subtowers.end() );
00566   std::sort( subtowers.begin(), subtowers.end(), subtower_has_greater_E );
00567 }
00568 
00569 // ---------------------------------------------------------------------------------------------------------
00570 // private helper functions to figure out detector & readout geometry as needed
00571 
00572 int reco::helper::JetIDHelper::HBHE_oddness( int iEta, int depth )
00573 {
00574  int ae = TMath::Abs( iEta );
00575  if( ae == 29 && depth == 1 ) ae += 1; // observed that: hits are at depths 1 & 2; 1 goes with the even pattern
00576  return ae & 0x1;
00577 }
00578 
00579 reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::HBHE_region( int iEta, int depth )
00580 {
00581   // no error checking for HO indices (depth 2 & |ieta|<=14 or depth 3 & |ieta|=15)
00582   if( iEta <= -17 || ( depth == 3 && iEta == -16 ) ) return HEneg;
00583   if( iEta >=  17 || ( depth == 3 && iEta ==  16 ) ) return HEpos;
00584   if( iEta < 0 ) return HBneg;
00585   return HBpos;
00586 }
00587 
00588 int reco::helper::JetIDHelper::HBHE_oddness( int iEta )
00589 {
00590  int ae = TMath::Abs( iEta );
00591  if( ae == 29 ) return -1; // can't figure it out without RecHits
00592  return ae & 0x1;
00593 }
00594 
00595 reco::helper::JetIDHelper::Region reco::helper::JetIDHelper::region( int iEta )
00596 {
00597   if( iEta == 16 || iEta == -16 ) return unknown_region; // both HB and HE cells belong to these towers
00598   if( iEta == 29 || iEta == -29 ) return unknown_region; // both HE and HF cells belong to these towers
00599   if( iEta <= -30 ) return HFneg;
00600   if( iEta >=  30 ) return HFpos;
00601   if( iEta <= -17 ) return HEneg;
00602   if( iEta >=  17 ) return HEpos;
00603   if( iEta < 0 ) return HBneg;
00604   return HBpos;
00605 }
CMSSW_4_4_3_patch1/src/RecoJets/JetProducers/src/JetIDHelper.cc
