Reference Manual

Go to the documentation of this file.
00001 /*
00002 Description: "Reflags" HB/HE hits based on their ECAL, HCAL, and tracking isolation.
00003 
00004 Original Author: John Paul Chou (Brown University)
00005                  Thursday, September 2, 2010
00006 */
00007 
00008 #include "HBHEIsolatedNoiseReflagger.h"
00009 
00010 #include "FWCore/Framework/interface/ESHandle.h"
00011 #include "FWCore/Framework/interface/EventSetup.h"
00012 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00013 #include "DataFormats/JetReco/interface/TrackExtrapolation.h"
00014 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00015 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
00016 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
00017 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"
00018 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00019 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00020 
00021 #include "RecoMET/METAlgorithms/interface/HcalHPDRBXMap.h"
00022 
00023 HBHEIsolatedNoiseReflagger::HBHEIsolatedNoiseReflagger(const edm::ParameterSet& iConfig) :
00024   hbheLabel_(iConfig.getParameter<edm::InputTag>("hbheInput")),
00025   ebLabel_(iConfig.getParameter<edm::InputTag>("ebInput")),
00026   eeLabel_(iConfig.getParameter<edm::InputTag>("eeInput")),
00027   trackExtrapolationLabel_(iConfig.getParameter<edm::InputTag>("trackExtrapolationInput")),
00028   
00029   LooseHcalIsol_(iConfig.getParameter<double>("LooseHcalIsol")),
00030   LooseEcalIsol_(iConfig.getParameter<double>("LooseEcalIsol")),
00031   LooseTrackIsol_(iConfig.getParameter<double>("LooseTrackIsol")),
00032   TightHcalIsol_(iConfig.getParameter<double>("TightHcalIsol")),
00033   TightEcalIsol_(iConfig.getParameter<double>("TightEcalIsol")),
00034   TightTrackIsol_(iConfig.getParameter<double>("TightTrackIsol")),
00035 
00036   LooseRBXEne1_(iConfig.getParameter<double>("LooseRBXEne1")),
00037   LooseRBXEne2_(iConfig.getParameter<double>("LooseRBXEne2")),
00038   LooseRBXHits1_(iConfig.getParameter<int>("LooseRBXHits1")),
00039   LooseRBXHits2_(iConfig.getParameter<int>("LooseRBXHits2")),
00040   TightRBXEne1_(iConfig.getParameter<double>("TightRBXEne1")),
00041   TightRBXEne2_(iConfig.getParameter<double>("TightRBXEne2")),
00042   TightRBXHits1_(iConfig.getParameter<int>("TightRBXHits1")),
00043   TightRBXHits2_(iConfig.getParameter<int>("TightRBXHits2")),
00044 
00045   LooseHPDEne1_(iConfig.getParameter<double>("LooseHPDEne1")),
00046   LooseHPDEne2_(iConfig.getParameter<double>("LooseHPDEne2")),
00047   LooseHPDHits1_(iConfig.getParameter<int>("LooseHPDHits1")),
00048   LooseHPDHits2_(iConfig.getParameter<int>("LooseHPDHits2")),
00049   TightHPDEne1_(iConfig.getParameter<double>("TightHPDEne1")),
00050   TightHPDEne2_(iConfig.getParameter<double>("TightHPDEne2")),
00051   TightHPDHits1_(iConfig.getParameter<int>("TightHPDHits1")),
00052   TightHPDHits2_(iConfig.getParameter<int>("TightHPDHits2")),
00053 
00054   LooseDiHitEne_(iConfig.getParameter<double>("LooseDiHitEne")),
00055   TightDiHitEne_(iConfig.getParameter<double>("TightDiHitEne")),
00056   LooseMonoHitEne_(iConfig.getParameter<double>("LooseMonoHitEne")),
00057   TightMonoHitEne_(iConfig.getParameter<double>("TightMonoHitEne")),
00058   
00059   debug_(iConfig.getUntrackedParameter<bool>("debug",true)),
00060   objvalidator_(iConfig)
00061 {
00062   produces<HBHERecHitCollection>();
00063 }
00064 
00065 HBHEIsolatedNoiseReflagger::~HBHEIsolatedNoiseReflagger()
00066 {
00067 }
00068 
00069 
00070 void
00071 HBHEIsolatedNoiseReflagger::produce(edm::Event& iEvent, const edm::EventSetup& evSetup)
00072 {
00073   // get the ECAL channel status map
00074   edm::ESHandle<EcalChannelStatus> ecalChStatus;
00075   evSetup.get<EcalChannelStatusRcd>().get( ecalChStatus );
00076   const EcalChannelStatus* dbEcalChStatus = ecalChStatus.product();
00077 
00078   // get the HCAL channel status map
00079   edm::ESHandle<HcalChannelQuality> hcalChStatus;    
00080   evSetup.get<HcalChannelQualityRcd>().get( hcalChStatus );
00081   const HcalChannelQuality* dbHcalChStatus = hcalChStatus.product();
00082 
00083   // get the severity level computers
00084   edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputerHndl;
00085   evSetup.get<HcalSeverityLevelComputerRcd>().get(hcalSevLvlComputerHndl);
00086   const HcalSeverityLevelComputer* hcalSevLvlComputer = hcalSevLvlComputerHndl.product();
00087 
00088   edm::ESHandle<EcalSeverityLevelAlgo> ecalSevLvlAlgoHndl;
00089   evSetup.get<EcalSeverityLevelAlgoRcd>().get(ecalSevLvlAlgoHndl);
00090   const EcalSeverityLevelAlgo* ecalSevLvlAlgo = ecalSevLvlAlgoHndl.product();
00091 
00092   // get the calotower mappings
00093   edm::ESHandle<CaloTowerConstituentsMap> ctcm;
00094   evSetup.get<IdealGeometryRecord>().get(ctcm);
00095   
00096   // get the HB/HE hits
00097   edm::Handle<HBHERecHitCollection> hbhehits_h;
00098   iEvent.getByLabel(hbheLabel_, hbhehits_h);
00099 
00100   // get the ECAL hits
00101   edm::Handle<EcalRecHitCollection> ebhits_h;
00102   iEvent.getByLabel(ebLabel_, ebhits_h);
00103   edm::Handle<EcalRecHitCollection> eehits_h;
00104   iEvent.getByLabel(eeLabel_, eehits_h);
00105 
00106   // get the tracks
00107   edm::Handle<std::vector<reco::TrackExtrapolation> > trackextraps_h;
00108   iEvent.getByLabel(trackExtrapolationLabel_, trackextraps_h);
00109 
00110   // set the status maps and severity level computers for the hit validator
00111   objvalidator_.setHcalChannelQuality(dbHcalChStatus);
00112   objvalidator_.setEcalChannelStatus(dbEcalChStatus);
00113   objvalidator_.setHcalSeverityLevelComputer(hcalSevLvlComputer);
00114   objvalidator_.setEcalSeverityLevelAlgo(ecalSevLvlAlgo);
00115   objvalidator_.setEBRecHitCollection(&(*ebhits_h));
00116   objvalidator_.setEERecHitCollection(&(*eehits_h));
00117 
00118   // organizer the hits
00119   PhysicsTowerOrganizer pto(iEvent, evSetup, hbhehits_h, ebhits_h, eehits_h, trackextraps_h, objvalidator_, *(ctcm.product()));
00120   HBHEHitMapOrganizer organizer(hbhehits_h, objvalidator_, pto);
00121 
00122   // get the rbxs, hpds, dihits, and monohits
00123   std::vector<HBHEHitMap> rbxs;
00124   std::vector<HBHEHitMap> hpds;
00125   std::vector<HBHEHitMap> dihits;
00126   std::vector<HBHEHitMap> monohits;
00127   organizer.getRBXs(rbxs, LooseRBXEne1_<TightRBXEne1_ ? LooseRBXEne1_ : TightRBXEne1_);
00128   organizer.getHPDs(hpds, LooseHPDEne1_<TightHPDEne1_ ? LooseHPDEne1_ : TightHPDEne1_);
00129   organizer.getDiHits(dihits, LooseDiHitEne_<TightDiHitEne_ ? LooseDiHitEne_ : TightDiHitEne_);
00130   organizer.getMonoHits(monohits, LooseMonoHitEne_<TightMonoHitEne_ ? LooseMonoHitEne_ : TightMonoHitEne_);
00131 
00132   if(debug_ && (rbxs.size()>0 || hpds.size()>0 || dihits.size()>0 || monohits.size()>0)) {
00133     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "RBXs:" << std::endl;
00134     DumpHBHEHitMap(rbxs);
00135     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "\nHPDs:" << std::endl;
00136     DumpHBHEHitMap(hpds);
00137     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "\nDiHits:" << std::endl;
00138     DumpHBHEHitMap(dihits);
00139     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "\nMonoHits:" << std::endl;
00140     DumpHBHEHitMap(monohits);
00141   }
00142 
00143   //  bool result=true;
00144 
00145   // determine which hits are noisy
00146   std::set<const HBHERecHit*> noisehits;
00147   for(int i=0; i<static_cast<int>(rbxs.size()); i++) {
00148     int nhits=rbxs[i].nHits();
00149     double ene=rbxs[i].hitEnergy();
00150     double trkfide=rbxs[i].hitEnergyTrackFiducial();
00151     double isolhcale=rbxs[i].hcalEnergySameTowers()+rbxs[i].hcalEnergyNeighborTowers();
00152     double isolecale=rbxs[i].ecalEnergySameTowers();
00153     double isoltrke=rbxs[i].trackEnergySameTowers()+rbxs[i].trackEnergyNeighborTowers();
00154     if((isolhcale/ene<LooseHcalIsol_ && isolecale/ene<LooseEcalIsol_ && isoltrke/ene<LooseTrackIsol_ && ((trkfide>LooseRBXEne1_ && nhits>=LooseRBXHits1_) || (trkfide>LooseRBXEne2_ && nhits>=LooseRBXHits2_))) ||
00155        (isolhcale/ene<TightHcalIsol_ && isolecale/ene<TightEcalIsol_ && isoltrke/ene<TightTrackIsol_ && ((trkfide>TightRBXEne1_ && nhits>=TightRBXHits1_) || (trkfide>TightRBXEne2_ && nhits>=TightRBXHits2_)))) {
00156       for(HBHEHitMap::hitmap_const_iterator it=rbxs[i].beginHits(); it!=rbxs[i].endHits(); ++it)
00157         noisehits.insert(it->first);
00158       //      result=false;
00159     }
00160   }
00161 
00162   for(int i=0; i<static_cast<int>(hpds.size()); i++) {
00163     int nhits=hpds[i].nHits();
00164     double ene=hpds[i].hitEnergy();
00165     double trkfide=hpds[i].hitEnergyTrackFiducial();
00166     double isolhcale=hpds[i].hcalEnergySameTowers()+hpds[i].hcalEnergyNeighborTowers();
00167     double isolecale=hpds[i].ecalEnergySameTowers();
00168     double isoltrke=hpds[i].trackEnergySameTowers()+hpds[i].trackEnergyNeighborTowers();
00169     if((isolhcale/ene<LooseHcalIsol_ && isolecale/ene<LooseEcalIsol_ && isoltrke/ene<LooseTrackIsol_ && ((trkfide>LooseHPDEne1_ && nhits>=LooseHPDHits1_) || (trkfide>LooseHPDEne2_ && nhits>=LooseHPDHits2_))) ||
00170        (isolhcale/ene<TightHcalIsol_ && isolecale/ene<TightEcalIsol_ && isoltrke/ene<TightTrackIsol_ && ((trkfide>TightHPDEne1_ && nhits>=TightHPDHits1_) || (trkfide>TightHPDEne2_ && nhits>=TightHPDHits2_)))) {
00171       for(HBHEHitMap::hitmap_const_iterator it=hpds[i].beginHits(); it!=hpds[i].endHits(); ++it)
00172         noisehits.insert(it->first);
00173       //      result=false;
00174     }
00175   }
00176 
00177   for(int i=0; i<static_cast<int>(dihits.size()); i++) {
00178     double ene=dihits[i].hitEnergy();
00179     double trkfide=dihits[i].hitEnergyTrackFiducial();
00180     double isolhcale=dihits[i].hcalEnergySameTowers()+dihits[i].hcalEnergyNeighborTowers();
00181     double isolecale=dihits[i].ecalEnergySameTowers();
00182     double isoltrke=dihits[i].trackEnergySameTowers()+dihits[i].trackEnergyNeighborTowers();
00183     if((isolhcale/ene<LooseHcalIsol_ && isolecale/ene<LooseEcalIsol_ && isoltrke/ene<LooseTrackIsol_ && trkfide>0.99*ene && trkfide>LooseDiHitEne_) ||
00184        (isolhcale/ene<TightHcalIsol_ && isolecale/ene<TightEcalIsol_ && isoltrke/ene<TightTrackIsol_ && ene>TightDiHitEne_)) {
00185       for(HBHEHitMap::hitmap_const_iterator it=dihits[i].beginHits(); it!=dihits[i].endHits(); ++it)
00186         noisehits.insert(it->first);
00187       //      result=false;
00188     }
00189   }
00190   
00191   for(int i=0; i<static_cast<int>(monohits.size()); i++) {
00192     double ene=monohits[i].hitEnergy();
00193     double trkfide=monohits[i].hitEnergyTrackFiducial();
00194     double isolhcale=monohits[i].hcalEnergySameTowers()+monohits[i].hcalEnergyNeighborTowers();
00195     double isolecale=monohits[i].ecalEnergySameTowers();
00196     double isoltrke=monohits[i].trackEnergySameTowers()+monohits[i].trackEnergyNeighborTowers();
00197     if((isolhcale/ene<LooseHcalIsol_ && isolecale/ene<LooseEcalIsol_ && isoltrke/ene<LooseTrackIsol_ && trkfide>0.99*ene && trkfide>LooseMonoHitEne_) ||
00198        (isolhcale/ene<TightHcalIsol_ && isolecale/ene<TightEcalIsol_ && isoltrke/ene<TightTrackIsol_ && ene>TightMonoHitEne_)) {
00199       for(HBHEHitMap::hitmap_const_iterator it=monohits[i].beginHits(); it!=monohits[i].endHits(); ++it)
00200         noisehits.insert(it->first);
00201       //      result=false;
00202     }
00203   }
00204 
00205   // prepare the output HBHE RecHit collection
00206   std::auto_ptr<HBHERecHitCollection> pOut(new HBHERecHitCollection());
00207   // loop over rechits, and set the new bit you wish to use
00208   for(HBHERecHitCollection::const_iterator it=hbhehits_h->begin(); it!=hbhehits_h->end(); ++it) {
00209     const HBHERecHit* hit=&(*it);
00210     HBHERecHit newhit(*hit);
00211     if(noisehits.end()!=noisehits.find(hit)) {
00212       newhit.setFlagField(1, HcalCaloFlagLabels::HBHEIsolatedNoise);
00213     }
00214     pOut->push_back(newhit);
00215   }
00216 
00217   iEvent.put(pOut);
00218 
00219   return;  
00220 }
00221 
00222 
00223 void HBHEIsolatedNoiseReflagger::DumpHBHEHitMap(std::vector<HBHEHitMap>& i) const
00224 {
00225   for(std::vector<HBHEHitMap>::const_iterator it=i.begin(); it!=i.end(); ++it) {
00226     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "hit energy=" << it->hitEnergy()
00227               << "; # hits=" << it->nHits()
00228               << "; hcal energy same=" << it->hcalEnergySameTowers()
00229                     << "; ecal energy same=" << it->ecalEnergySameTowers()
00230                   << "; track energy same=" << it->trackEnergySameTowers()
00231                   << "; neighbor hcal energy=" << it->hcalEnergyNeighborTowers() << std::endl;
00232         edm::LogInfo("HBHEIsolatedNoiseReflagger") << "hits:" << std::endl;
00233         for(HBHEHitMap::hitmap_const_iterator it2=it->beginHits(); it2!=it->endHits(); ++it2) {
00234           const HBHERecHit *hit=it2->first;
00235             edm::LogInfo("HBHEIsolatedNoiseReflagger") << "RBX #=" << HcalHPDRBXMap::indexRBX(hit->id())
00236                       << "; HPD #=" << HcalHPDRBXMap::indexHPD(hit->id())
00237                       << "; " << (*hit) << std::endl;
00238         }
00239   }
00240   return;
00241 }
00242 
00243 
00244 //define this as a plug-in
00245 DEFINE_FWK_MODULE(HBHEIsolatedNoiseReflagger);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoLocalCalo/HcalRecProducers/src/HBHEIsolatedNoiseReflagger.cc
