Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Class:      HLTHcalMETNoiseCleaner
00004 // 
00012 //
00013 // Original Author:  Alexander Mott
00014 //         Created:  Mon Nov 21 11:32:00 CEST 2011
00015 // 
00016 //
00017 //
00018 
00019 #include "HLTrigger/JetMET/interface/HLTHcalMETNoiseCleaner.h"
00020 
00021 #include "DataFormats/Common/interface/Handle.h"
00022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00023 #include "FWCore/Framework/interface/ESHandle.h"
00024 #include "FWCore/Framework/interface/Event.h"
00025 #include "FWCore/Framework/interface/EventSetup.h"
00026 
00027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00028 
00029 #include "DataFormats/Candidate/interface/Candidate.h"
00030 #include "DataFormats/METReco/interface/HcalNoiseRBX.h"
00031 #include "DataFormats/METReco/interface/CaloMET.h"
00032 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00033 #include "DataFormats/METReco/interface/SpecificCaloMETData.h"
00034 #include "DataFormats/CaloTowers/interface/CaloTower.h"
00035 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00036 #include "DataFormats/Math/interface/LorentzVector.h"
00037 #include "DataFormats/Math/interface/Point3D.h"
00038 
00039 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
00040 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
00041 #include "FWCore/Utilities/interface/InputTag.h"
00042 
00043 #include <iostream>
00044 #include <string>
00045 #include <fstream>
00046 #include <TVector3.h>
00047 #include <TLorentzVector.h>
00048 //#include <Point.h>
00049 
00050 HLTHcalMETNoiseCleaner::HLTHcalMETNoiseCleaner(const edm::ParameterSet& iConfig)
00051   : HcalNoiseRBXCollectionTag_(iConfig.getParameter<edm::InputTag>("HcalNoiseRBXCollection")),
00052     CaloMetCollectionTag_(iConfig.getParameter<edm::InputTag>("CaloMetCollection")),
00053     CaloMetCut_(iConfig.getParameter<double>("CaloMetCut")),
00054     severity_(iConfig.getParameter<int> ("severity")),
00055     maxNumRBXs_(iConfig.getParameter<int>("maxNumRBXs")),
00056     numRBXsToConsider_(iConfig.getParameter<int>("numRBXsToConsider")),
00057     accept2NoiseRBXEvents_(iConfig.getParameter<bool>("accept2NoiseRBXEvents")),
00058     needEMFCoincidence_(iConfig.getParameter<bool>("needEMFCoincidence")),
00059     minRBXEnergy_(iConfig.getParameter<double>("minRBXEnergy")),
00060     minRatio_(iConfig.getParameter<double>("minRatio")),
00061     maxRatio_(iConfig.getParameter<double>("maxRatio")),
00062     minHPDHits_(iConfig.getParameter<int>("minHPDHits")),
00063     minRBXHits_(iConfig.getParameter<int>("minRBXHits")),
00064     minHPDNoOtherHits_(iConfig.getParameter<int>("minHPDNoOtherHits")),
00065     minZeros_(iConfig.getParameter<int>("minZeros")),
00066     minHighEHitTime_(iConfig.getParameter<double>("minHighEHitTime")),
00067     maxHighEHitTime_(iConfig.getParameter<double>("maxHighEHitTime")),
00068     maxRBXEMF_(iConfig.getParameter<double>("maxRBXEMF")),
00069     minRecHitE_(iConfig.getParameter<double>("minRecHitE")),
00070     minLowHitE_(iConfig.getParameter<double>("minLowHitE")),
00071     minHighHitE_(iConfig.getParameter<double>("minHighHitE")),
00072     TS4TS5EnergyThreshold_(iConfig.getParameter<double>("TS4TS5EnergyThreshold"))
00073 {
00074 
00075   std::vector<double> TS4TS5UpperThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperThreshold");
00076   std::vector<double> TS4TS5UpperCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperCut");
00077   std::vector<double> TS4TS5LowerThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerThreshold");
00078   std::vector<double> TS4TS5LowerCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerCut");
00079 
00080   for(int i = 0; i < (int)TS4TS5UpperThresholdTemp.size() && i < (int)TS4TS5UpperCutTemp.size(); i++)
00081      TS4TS5UpperCut_.push_back(std::pair<double, double>(TS4TS5UpperThresholdTemp[i], TS4TS5UpperCutTemp[i]));
00082   sort(TS4TS5UpperCut_.begin(), TS4TS5UpperCut_.end());
00083 
00084   for(int i = 0; i < (int)TS4TS5LowerThresholdTemp.size() && i < (int)TS4TS5LowerCutTemp.size(); i++)
00085      TS4TS5LowerCut_.push_back(std::pair<double, double>(TS4TS5LowerThresholdTemp[i], TS4TS5LowerCutTemp[i]));
00086   sort(TS4TS5LowerCut_.begin(), TS4TS5LowerCut_.end());
00087 
00088   produces<reco::CaloMETCollection>();
00089 }
00090 
00091 
00092 HLTHcalMETNoiseCleaner::~HLTHcalMETNoiseCleaner(){}
00093 
00094 void
00095 HLTHcalMETNoiseCleaner::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
00096   edm::ParameterSetDescription desc;
00097   desc.add<edm::InputTag>("HcalNoiseRBXCollection",edm::InputTag("hltHcalNoiseInfoProducer"));
00098   desc.add<edm::InputTag>("CaloMetCollection",edm::InputTag("hltMet"));
00099   desc.add<double>("CaloMetCut",0.0);
00100   desc.add<int>("severity",1);
00101   desc.add<int>("maxNumRBXs",2);
00102   desc.add<int>("numRBXsToConsider",2);
00103   desc.add<bool>("accept2NoiseRBXEvents",true);
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("hltHcalMETNoiseCleaner",desc);
00134 }
00135 
00136 //
00137 // member functions
00138 //
00139 
00140 bool HLTHcalMETNoiseCleaner::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
00141 {
00142   using namespace reco;
00143 
00144   //output collection
00145   std::auto_ptr<CaloMETCollection> CleanedMET(new CaloMETCollection);
00146 
00147   //get the calo MET / MHT
00148   edm::Handle<CaloMETCollection> met_h;
00149   iEvent.getByLabel(CaloMetCollectionTag_,met_h);
00150   
00151   if(not met_h.isValid() or met_h->size()==0 or met_h->front().pt()<0){ //No Valid MET, don't do anything and accept the event
00152     return true;  // we shouldn't get here, but lets not crash
00153   }
00154   
00155   reco::CaloMET inCaloMet = met_h->front();
00156 
00157 
00158   // in this case, do not filter anything
00159   if(severity_==0){
00160     CleanedMET->push_back(inCaloMet);
00161     iEvent.put(CleanedMET);
00162     return true;
00163   }
00164 
00165   // get the RBXs produced by RecoMET/METProducers/HcalNoiseInfoProducer
00166   edm::Handle<HcalNoiseRBXCollection> rbxs_h;
00167   iEvent.getByLabel(HcalNoiseRBXCollectionTag_,rbxs_h);
00168   if(!rbxs_h.isValid()) {
00169     edm::LogError("DataNotFound") << "HLTHcalMETNoiseCleaner: Could not find HcalNoiseRBXCollection product named "
00170                                   << HcalNoiseRBXCollectionTag_ << "." << std::endl;
00171     CleanedMET->push_back(inCaloMet);
00172     iEvent.put(CleanedMET);    
00173     return true; // no valid RBXs
00174   }
00175 
00176   // create a sorted set of the RBXs, ordered by energy
00177   noisedataset_t data;
00178   for(HcalNoiseRBXCollection::const_iterator it=rbxs_h->begin(); it!=rbxs_h->end(); ++it) {
00179     const HcalNoiseRBX &rbx=(*it);
00180     CommonHcalNoiseRBXData d(rbx, minRecHitE_, minLowHitE_, minHighHitE_, TS4TS5EnergyThreshold_,
00181                              TS4TS5UpperCut_, TS4TS5LowerCut_);
00182     data.insert(d);
00183   }
00184   //if 0 RBXs are in the list, just accept
00185   if(data.size()<1){
00186     CleanedMET->push_back(inCaloMet);
00187     iEvent.put(CleanedMET);
00188     return true;
00189   }
00190   // data is now sorted by RBX energy
00191   // only consider top N=numRBXsToConsider_ energy RBXs
00192   int cntr=0;
00193   int nNoise=0;
00194 
00195   TVector3 metVec;
00196   metVec.SetPtEtaPhi(met_h->front().pt(), 0, met_h->front().phi() );
00197 
00198   TVector3 noiseHPDVector(0,0,0);
00199   TVector3 secondHPDVector(0,0,0);
00200   for(noisedataset_t::const_iterator it=data.begin();
00201       it!=data.end() && cntr<numRBXsToConsider_;
00202       it++, cntr++) {
00203     bool isNoise=false;
00204     bool passFilter=true;
00205     bool passEMF=true;
00206     if(it->energy()>minRBXEnergy_) {
00207       if(it->validRatio() && it->ratio()<minRatio_)        passFilter=false;
00208       else if(it->validRatio() && it->ratio()>maxRatio_)   passFilter=false;
00209       else if(it->numHPDHits()>=minHPDHits_)               passFilter=false;
00210       else if(it->numRBXHits()>=minRBXHits_)               passFilter=false;
00211       else if(it->numHPDNoOtherHits()>=minHPDNoOtherHits_) passFilter=false;
00212       else if(it->numZeros()>=minZeros_)                   passFilter=false;
00213       else if(it->minHighEHitTime()<minHighEHitTime_)      passFilter=false;
00214       else if(it->maxHighEHitTime()>maxHighEHitTime_)      passFilter=false;
00215       else if(!it->PassTS4TS5())                           passFilter=false;
00216 
00217       if(it->RBXEMF()<maxRBXEMF_){
00218         passEMF=false;
00219       }
00220     }
00221     
00222     if((needEMFCoincidence_ && !passEMF && !passFilter) ||
00223        (!needEMFCoincidence_ && !passFilter)) { // check for noise
00224       LogDebug("") << "HLTHcalMETNoiseCleaner debug: Found a noisy RBX: "
00225                    << "energy=" << it->energy() << "; "
00226                    << "ratio=" << it->ratio() << "; "
00227                    << "# RBX hits=" << it->numRBXHits() << "; "
00228                    << "# HPD hits=" << it->numHPDHits() << "; "
00229                    << "# Zeros=" << it->numZeros() << "; "
00230                    << "min time=" << it->minHighEHitTime() << "; "
00231                    << "max time=" << it->maxHighEHitTime() << "; "
00232                    << "passTS4TS5=" << it->PassTS4TS5() << "; "
00233                    << "RBX EMF=" << it->RBXEMF()
00234                    << std::endl;
00235       nNoise++;
00236       isNoise=true;
00237     }// OK, checked for noise
00238 
00239     //------------First Noisy RBX-----------------------
00240     if(isNoise && nNoise==1){
00241       edm::RefVector<CaloTowerCollection> noiseTowers = it->rbxTowers();
00242       edm::RefVector<CaloTowerCollection>::const_iterator noiseTowersIt;
00243       // get the energy vector for this RBX from the calotowers
00244       for( noiseTowersIt = noiseTowers.begin(); noiseTowersIt != noiseTowers.end(); noiseTowersIt++){
00245         TVector3 towerVec;
00246         towerVec.SetPtEtaPhi((*noiseTowersIt)->pt(),(*noiseTowersIt)->eta(),(*noiseTowersIt)->phi());
00247         noiseHPDVector+=towerVec; // add this tower to the vector for the RBX
00248       }
00249       if(noiseHPDVector.Mag()>0) noiseHPDVector.SetPtEtaPhi(noiseHPDVector.Pt(),0,noiseHPDVector.Phi()); // make the noise transverse
00250       else noiseHPDVector.SetPtEtaPhi(0,0,0);
00251     }
00252     //-----------FOUND a SECOND NOISY RBX-------------------
00253     if(isNoise && cntr > 0){ 
00254         CleanedMET->push_back(inCaloMet);
00255         iEvent.put(CleanedMET);         
00256         return accept2NoiseRBXEvents_; // don't try to clean these for the moment, just keep or throw away      
00257     }
00258     //----------LEADING RBX is NOT NOISY--------------------
00259     if(!isNoise && cntr == 0){ 
00260       CleanedMET->push_back(inCaloMet);
00261       iEvent.put(CleanedMET);           
00262       return true; // don't reject the event if the leading RBX isn't noise
00263     }
00264     //-----------SUBLEADING RBX is NOT NOISY: STORE INFO----
00265     if(!isNoise && nNoise>0){ //second RBX isn't noisy (and first one was), so clean 
00266       edm::RefVector<CaloTowerCollection> noiseTowers = it->rbxTowers();
00267       edm::RefVector<CaloTowerCollection>::const_iterator noiseTowersIt;
00268         for( noiseTowersIt = noiseTowers.begin(); noiseTowersIt != noiseTowers.end(); noiseTowersIt++){
00269           TVector3 towerVec;
00270           towerVec.SetPtEtaPhi((*noiseTowersIt)->pt(),(*noiseTowersIt)->eta(),(*noiseTowersIt)->phi());
00271           secondHPDVector+=towerVec; // add this tower to the vector for the RBX
00272         }
00273         if(secondHPDVector.Mag()>0) secondHPDVector.SetPtEtaPhi(secondHPDVector.Pt(),0,secondHPDVector.Phi()); // make the second transverse
00274         else secondHPDVector.SetPtEtaPhi(0,0,0);
00275         break;
00276     }
00277   } // end RBX loop
00278 
00279   if(noiseHPDVector.Mag()==0){
00280     CleanedMET->push_back(inCaloMet);
00281     iEvent.put(CleanedMET);     
00282     return true; // don't reject the event if the leading RBX isn't noise
00283   }
00284 
00285   //********************************************************************************
00286   //The Event gets here only if it had exactly 1 noisy RBX in the lead position
00287   //********************************************************************************
00288 
00289   float METsumet = met_h->front().energy();
00290 
00291   metVec+=noiseHPDVector;
00292 
00293   float ZMETsumet = METsumet-noiseHPDVector.Mag();
00294   float ZMETpt = metVec.Pt();
00295   float ZMETphi = metVec.Phi();
00296 
00297   //put the second RBX vector in the eta phi position of the leading RBX vector
00298  
00299   float SMETsumet = 0;
00300   float SMETpt = 0;
00301   float SMETphi = 0;
00302  if(secondHPDVector.Mag()>0.){
00303     secondHPDVector.SetPtEtaPhi(secondHPDVector.Pt(),noiseHPDVector.Eta(),noiseHPDVector.Phi());
00304     metVec-= secondHPDVector;
00305     SMETsumet = METsumet-noiseHPDVector.Mag();
00306     SMETpt = metVec.Pt();
00307     SMETphi = metVec.Phi();  
00308   }
00309   //Get the maximum MET:
00310   float CorMetSumEt,CorMetPt,CorMetPhi;
00311   if(ZMETpt>SMETpt){
00312     CorMetSumEt = ZMETsumet;
00313     CorMetPt = ZMETpt;
00314     CorMetPhi = ZMETphi;
00315   }else{
00316     CorMetSumEt = SMETsumet;
00317     CorMetPt = SMETpt;
00318     CorMetPhi = SMETphi;
00319   }
00320 
00321   reco::CaloMET corMet = BuildCaloMet(CorMetSumEt,CorMetPt,CorMetPhi);
00322   CleanedMET->push_back(corMet);
00323   iEvent.put(CleanedMET);
00324 
00325   return (corMet.pt() > CaloMetCut_);
00326 }
00327 
00328 reco::CaloMET HLTHcalMETNoiseCleaner::BuildCaloMet(float sumet, float pt, float phi){
00329   // Instantiate the container to hold the calorimeter specific information
00330   
00331   typedef math::XYZPoint Point;
00332   typedef math::XYZTLorentzVector LorentzVector;
00333   
00334   SpecificCaloMETData specific;
00335   // Initialise the container 
00336   specific.MaxEtInEmTowers = 0.0;         // Maximum energy in EM towers
00337   specific.MaxEtInHadTowers = 0.0;        // Maximum energy in HCAL towers
00338   specific.HadEtInHO = 0.0;          // Hadronic energy fraction in HO
00339   specific.HadEtInHB = 0.0;          // Hadronic energy in HB
00340   specific.HadEtInHF = 0.0;          // Hadronic energy in HF
00341   specific.HadEtInHE = 0.0;          // Hadronic energy in HE
00342   specific.EmEtInEB = 0.0;           // Em energy in EB
00343   specific.EmEtInEE = 0.0;           // Em energy in EE
00344   specific.EmEtInHF = 0.0;           // Em energy in HF
00345   specific.EtFractionHadronic = 0.0; // Hadronic energy fraction
00346   specific.EtFractionEm = 0.0;       // Em energy fraction
00347   specific.CaloSETInpHF = 0.0;        // CaloSET in HF+ 
00348   specific.CaloSETInmHF = 0.0;        // CaloSET in HF- 
00349   specific.CaloMETInpHF = 0.0;        // CaloMET in HF+ 
00350   specific.CaloMETInmHF = 0.0;        // CaloMET in HF- 
00351   specific.CaloMETPhiInpHF = 0.0;     // CaloMET-phi in HF+ 
00352   specific.CaloMETPhiInmHF = 0.0;     // CaloMET-phi in HF- 
00353   specific.METSignificance = 0.0;
00354 
00355   TLorentzVector p4TL;
00356   p4TL.SetPtEtaPhiM(pt,0.,phi,0.);
00357   const LorentzVector p4(p4TL.X(),p4TL.Y(),0,p4TL.T());
00358   const Point vtx( 0.0, 0.0, 0.0 );
00359   reco::CaloMET specificmet( specific, sumet, p4, vtx );
00360   return specificmet;
00361   }
00362
/data/refman/pasoursint/CMSSW_5_3_10_patch1/src/HLTrigger/JetMET/src/HLTHcalMETNoiseCleaner.cc
