CMS 3D CMS Logo

/data/git/CMSSW_5_3_11_patch5/src/HLTrigger/JetMET/src/HLTHcalTowerNoiseCleaner.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Class:      HLTHcalTowerNoiseCleaner
00004 // 
00012 //
00013 // Original Author:  Alexander Mott
00014 //         Created:  Mon Nov 21 11:32:00 CEST 2011
00015 // 
00016 //
00017 //
00018 
00019 
00020 #include "DataFormats/Common/interface/Handle.h"
00021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00022 #include "FWCore/Framework/interface/ESHandle.h"
00023 #include "FWCore/Framework/interface/Event.h"
00024 #include "FWCore/Framework/interface/EventSetup.h"
00025 
00026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00027 
00028 #include "DataFormats/Candidate/interface/Candidate.h"
00029 #include "DataFormats/METReco/interface/HcalNoiseRBX.h"
00030 #include "DataFormats/JetReco/interface/Jet.h"
00031 #include "DataFormats/JetReco/interface/CaloJet.h"
00032 #include "DataFormats/CaloTowers/interface/CaloTower.h"
00033 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
00034 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00035 #include "DataFormats/Math/interface/LorentzVector.h"
00036 #include "DataFormats/Math/interface/Point3D.h"
00037 
00038 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
00039 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
00040 #include "FWCore/Utilities/interface/InputTag.h"
00041 
00042 #include <iostream>
00043 #include <string>
00044 #include <fstream>
00045 #include <TVector3.h>
00046 #include <TLorentzVector.h>
00047 #include <set>
00048 
00049 #include "HLTrigger/JetMET/interface/HLTHcalTowerNoiseCleaner.h"
00050 
00051 //#include <Point.h>
00052 
00053 HLTHcalTowerNoiseCleaner::HLTHcalTowerNoiseCleaner(const edm::ParameterSet& iConfig)
00054   : HcalNoiseRBXCollectionTag_(iConfig.getParameter<edm::InputTag>("HcalNoiseRBXCollection")),
00055     TowerCollectionTag_(iConfig.getParameter<edm::InputTag>("CaloTowerCollection")),
00056     severity_(iConfig.getParameter<int> ("severity")),
00057     maxNumRBXs_(iConfig.getParameter<int>("maxNumRBXs")),
00058     numRBXsToConsider_(iConfig.getParameter<int>("numRBXsToConsider")),
00059     needEMFCoincidence_(iConfig.getParameter<bool>("needEMFCoincidence")),
00060     minRBXEnergy_(iConfig.getParameter<double>("minRBXEnergy")),
00061     minRatio_(iConfig.getParameter<double>("minRatio")),
00062     maxRatio_(iConfig.getParameter<double>("maxRatio")),
00063     minHPDHits_(iConfig.getParameter<int>("minHPDHits")),
00064     minRBXHits_(iConfig.getParameter<int>("minRBXHits")),
00065     minHPDNoOtherHits_(iConfig.getParameter<int>("minHPDNoOtherHits")),
00066     minZeros_(iConfig.getParameter<int>("minZeros")),
00067     minHighEHitTime_(iConfig.getParameter<double>("minHighEHitTime")),
00068     maxHighEHitTime_(iConfig.getParameter<double>("maxHighEHitTime")),
00069     maxRBXEMF_(iConfig.getParameter<double>("maxRBXEMF")),
00070     minRecHitE_(iConfig.getParameter<double>("minRecHitE")),
00071     minLowHitE_(iConfig.getParameter<double>("minLowHitE")),
00072     minHighHitE_(iConfig.getParameter<double>("minHighHitE")),
00073     TS4TS5EnergyThreshold_(iConfig.getParameter<double>("TS4TS5EnergyThreshold"))
00074 {
00075 
00076   std::vector<double> TS4TS5UpperThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperThreshold");
00077   std::vector<double> TS4TS5UpperCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperCut");
00078   std::vector<double> TS4TS5LowerThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerThreshold");
00079   std::vector<double> TS4TS5LowerCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerCut");
00080 
00081   for(int i = 0; i < (int)TS4TS5UpperThresholdTemp.size() && i < (int)TS4TS5UpperCutTemp.size(); i++)
00082      TS4TS5UpperCut_.push_back(std::pair<double, double>(TS4TS5UpperThresholdTemp[i], TS4TS5UpperCutTemp[i]));
00083   sort(TS4TS5UpperCut_.begin(), TS4TS5UpperCut_.end());
00084 
00085   for(int i = 0; i < (int)TS4TS5LowerThresholdTemp.size() && i < (int)TS4TS5LowerCutTemp.size(); i++)
00086      TS4TS5LowerCut_.push_back(std::pair<double, double>(TS4TS5LowerThresholdTemp[i], TS4TS5LowerCutTemp[i]));
00087   sort(TS4TS5LowerCut_.begin(), TS4TS5LowerCut_.end());
00088 
00089   produces<CaloTowerCollection>();
00090 }
00091 
00092 
00093 HLTHcalTowerNoiseCleaner::~HLTHcalTowerNoiseCleaner(){}
00094 
00095 void
00096 HLTHcalTowerNoiseCleaner::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
00097   edm::ParameterSetDescription desc;
00098   desc.add<edm::InputTag>("HcalNoiseRBXCollection",edm::InputTag("hltHcalNoiseInfoProducer"));
00099   desc.add<edm::InputTag>("CaloTowerCollection",edm::InputTag("hltTowerMakerForAll"));
00100   desc.add<double>("maxTowerNoiseEnergyFraction",0.5);
00101   desc.add<int>("severity",1);
00102   desc.add<int>("maxNumRBXs",2);
00103   desc.add<int>("numRBXsToConsider",2);
00104   desc.add<bool>("needEMFCoincidence",true);
00105   desc.add<double>("minRBXEnergy",50.0);
00106   desc.add<double>("minRatio",-999.);
00107   desc.add<double>("maxRatio",999.);
00108   desc.add<int>("minHPDHits",17);
00109   desc.add<int>("minRBXHits",999);
00110   desc.add<int>("minHPDNoOtherHits",10);
00111   desc.add<int>("minZeros",10);
00112   desc.add<double>("minHighEHitTime",-9999.0);
00113   desc.add<double>("maxHighEHitTime",9999.0);
00114   desc.add<double>("maxRBXEMF",0.02);
00115   desc.add<double>("minRecHitE",1.5);
00116   desc.add<double>("minLowHitE",10.0);
00117   desc.add<double>("minHighHitE",25.0);
00118   desc.add<double>("TS4TS5EnergyThreshold",50.0);
00119 
00120   double TS4TS5UpperThresholdArray[5] = {70, 90, 100, 400, 4000 };
00121   double TS4TS5UpperCutArray[5] = {1, 0.8, 0.75, 0.72, 0.72};
00122   double TS4TS5LowerThresholdArray[7] = {100, 120, 150, 200, 300, 400, 500};
00123   double TS4TS5LowerCutArray[7] = {-1, -0.7, -0.4, -0.2, -0.08, 0, 0.1};
00124   std::vector<double> TS4TS5UpperThreshold(TS4TS5UpperThresholdArray, TS4TS5UpperThresholdArray+5);
00125   std::vector<double> TS4TS5UpperCut(TS4TS5UpperCutArray, TS4TS5UpperCutArray+5);
00126   std::vector<double> TS4TS5LowerThreshold(TS4TS5LowerThresholdArray, TS4TS5LowerThresholdArray+7);
00127   std::vector<double> TS4TS5LowerCut(TS4TS5LowerCutArray, TS4TS5LowerCutArray+7);
00128 
00129   desc.add<std::vector<double> >("TS4TS5UpperThreshold", TS4TS5UpperThreshold);
00130   desc.add<std::vector<double> >("TS4TS5UpperCut", TS4TS5UpperCut);
00131   desc.add<std::vector<double> >("TS4TS5LowerThreshold", TS4TS5LowerThreshold);
00132   desc.add<std::vector<double> >("TS4TS5LowerCut", TS4TS5LowerCut);
00133   descriptions.add("hltHcalTowerNoiseCleaner",desc);
00134 }
00135 
00136 //
00137 // member functions
00138 //
00139 
00140 void HLTHcalTowerNoiseCleaner::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00141 {
00142   using namespace reco;
00143 
00144   //get the calo MET / MHT
00145   edm::Handle<CaloTowerCollection> tower_h;
00146   iEvent.getByLabel(TowerCollectionTag_,tower_h);
00147   
00148   std::set<unsigned int> noisyTowers;
00149 
00150   if(not tower_h.isValid()){ //No towers MET, don't do anything and accept the event
00151     edm::LogError("HLTHcalTowerNoiseCleaner") << "Input Tower Collection is not Valid";
00152     return;
00153   }
00154 
00155   
00156   // get the RBXs produced by RecoMET/METProducers/HcalNoiseInfoProducer
00157   edm::Handle<HcalNoiseRBXCollection> rbxs_h;
00158   iEvent.getByLabel(HcalNoiseRBXCollectionTag_,rbxs_h);
00159   if(!rbxs_h.isValid()) {
00160     edm::LogWarning("HLTHcalTowerNoiseCleaner") << "Could not find HcalNoiseRBXCollection product named "
00161                                               << HcalNoiseRBXCollectionTag_ << "." << std::endl;
00162     severity_=0;
00163   }
00164                 
00165   // create a sorted set of the RBXs, ordered by energy
00166   noisedataset_t data;
00167   for(HcalNoiseRBXCollection::const_iterator it=rbxs_h->begin(); it!=rbxs_h->end(); ++it) {
00168     const HcalNoiseRBX &rbx=(*it);
00169     CommonHcalNoiseRBXData d(rbx, minRecHitE_, minLowHitE_, minHighHitE_, TS4TS5EnergyThreshold_,
00170                              TS4TS5UpperCut_, TS4TS5LowerCut_);
00171     data.insert(d);
00172   }
00173                 
00174   // data is now sorted by RBX energy
00175   // only consider top N=numRBXsToConsider_ energy RBXs
00176   if(severity_>0){
00177     for(noisedataset_t::const_iterator it=data.begin();
00178         it!=data.end();
00179         it++) {
00180       
00181       
00182       bool passFilter=true;
00183       bool passEMF=true;
00184       if(it->energy()>minRBXEnergy_) {
00185         if(it->validRatio() && it->ratio()<minRatio_)        passFilter=false;
00186         else if(it->validRatio() && it->ratio()>maxRatio_)   passFilter=false;
00187         else if(it->numHPDHits()>=minHPDHits_)               passFilter=false;
00188         else if(it->numRBXHits()>=minRBXHits_)               passFilter=false;
00189         else if(it->numHPDNoOtherHits()>=minHPDNoOtherHits_) passFilter=false;
00190         else if(it->numZeros()>=minZeros_)                   passFilter=false;
00191         else if(it->minHighEHitTime()<minHighEHitTime_)      passFilter=false;
00192         else if(it->maxHighEHitTime()>maxHighEHitTime_)      passFilter=false;
00193         else if(!it->PassTS4TS5())                           passFilter=false;
00194         
00195         if(it->RBXEMF()<maxRBXEMF_){
00196           passEMF=false;
00197         }
00198       }
00199       
00200       if((needEMFCoincidence_ && !passEMF && !passFilter) ||
00201          (!needEMFCoincidence_ && !passFilter)) { // check for noise
00202         LogDebug("") << "HLTHcalTowerNoiseCleaner debug: Found a noisy RBX: "
00203                      << "energy=" << it->energy() << "; "
00204                      << "ratio=" << it->ratio() << "; "
00205                      << "# RBX hits=" << it->numRBXHits() << "; "
00206                      << "# HPD hits=" << it->numHPDHits() << "; "
00207                      << "# Zeros=" << it->numZeros() << "; "
00208                      << "min time=" << it->minHighEHitTime() << "; "
00209                      << "max time=" << it->maxHighEHitTime() << "; "
00210                      << "passTS4TS5=" << it->PassTS4TS5() << "; "
00211                      << "RBX EMF=" << it->RBXEMF()
00212                      << std::endl;
00213         // add calotowers associated with this RBX to the noise list
00214         edm::RefVector<CaloTowerCollection> noiseTowers = it->rbxTowers();
00215         edm::RefVector<CaloTowerCollection>::const_iterator noiseTowersIt;
00216         //add these calotowers to the noisy list
00217         for( noiseTowersIt = noiseTowers.begin(); noiseTowersIt != noiseTowers.end(); noiseTowersIt++){
00218           edm::Ref<edm::SortedCollection<CaloTower> > tower_ref = *noiseTowersIt;
00219           CaloTowerDetId id = tower_ref->id();
00220           noisyTowers.insert( id.denseIndex() );
00221         }}
00222     } // done with noise loop
00223   }//if(severity_>0)
00224   
00225   //output collection
00226   std::auto_ptr<CaloTowerCollection> OutputTowers(new CaloTowerCollection() );
00227 
00228   CaloTowerCollection::const_iterator inTowersIt;
00229   
00230   for(inTowersIt = tower_h->begin(); inTowersIt != tower_h->end(); inTowersIt++){
00231     const CaloTower & tower = (*inTowersIt);
00232     CaloTowerDetId id = tower.id();
00233     if(noisyTowers.find( id.denseIndex() ) == noisyTowers.end()){ // the tower is not noisy
00234       OutputTowers->push_back(*inTowersIt);
00235     }
00236   }
00237   iEvent.put(OutputTowers);
00238 
00239 }
00240 
00241