Reference Manual

Go to the documentation of this file.
00001 #include "RecoMET/METAlgorithms/interface/HcalNoiseAlgo.h"
00002 
00003 CommonHcalNoiseRBXData::CommonHcalNoiseRBXData(const reco::HcalNoiseRBX& rbx, double minRecHitE,
00004    double minLowHitE, double minHighHitE, double TS4TS5EnergyThreshold,
00005    std::vector<std::pair<double, double> > &TS4TS5UpperCut,
00006    std::vector<std::pair<double, double> > &TS4TS5LowerCut)
00007 {
00008   // energy
00009   energy_ = rbx.recHitEnergy(minRecHitE); 
00010 
00011   // ratio
00012   e2ts_ = rbx.allChargeHighest2TS();
00013   e10ts_ = rbx.allChargeTotal();
00014 
00015   // TS4TS5
00016   TS4TS5Decision_ = true;
00017   if(energy_ > TS4TS5EnergyThreshold)   // check filter
00018   {
00019      std::vector<float> AllCharge = rbx.allCharge();
00020      double BaseCharge = AllCharge[4] + AllCharge[5];
00021      if(BaseCharge < 1)
00022         BaseCharge = 1;
00023      double TS4TS5 = (AllCharge[4] - AllCharge[5]) / BaseCharge;
00024 
00025      if(CheckPassFilter(BaseCharge, TS4TS5, TS4TS5UpperCut, 1) == false)
00026         TS4TS5Decision_ = false;
00027      if(CheckPassFilter(BaseCharge, TS4TS5, TS4TS5LowerCut, -1) == false)
00028         TS4TS5Decision_ = false;
00029   }
00030   else
00031      TS4TS5Decision_ = true;
00032 
00033   // # of hits
00034   numHPDHits_ = 0;
00035   for(std::vector<reco::HcalNoiseHPD>::const_iterator it1=rbx.HPDsBegin(); it1!=rbx.HPDsEnd(); ++it1) {
00036     int nhpdhits=it1->numRecHits(minRecHitE);
00037     if(numHPDHits_ < nhpdhits) numHPDHits_ = nhpdhits;
00038   }
00039   numRBXHits_ = rbx.numRecHits(minRecHitE);
00040   numHPDNoOtherHits_ = (numHPDHits_ == numRBXHits_) ? numHPDHits_ : 0;    
00041 
00042   // # of ADC zeros
00043   numZeros_ = rbx.totalZeros();
00044 
00045   // timing
00046   minLowEHitTime_ = minHighEHitTime_ = 99999.;
00047   maxLowEHitTime_ = maxHighEHitTime_ = -99999.;
00048   lowEHitTimeSqrd_ = highEHitTimeSqrd_ = 0;
00049   numLowEHits_ = numHighEHits_ = 0;
00050   for(std::vector<reco::HcalNoiseHPD>::const_iterator it1=rbx.HPDsBegin(); it1!=rbx.HPDsEnd(); ++it1) {
00051     edm::RefVector<HBHERecHitCollection> rechits=it1->recHits();
00052     for(edm::RefVector<HBHERecHitCollection>::const_iterator it2=rechits.begin(); it2!=rechits.end(); ++it2) {
00053       float energy=(*it2)->energy();
00054       float time=(*it2)->time();
00055       if(energy>=minLowHitE) {
00056         if(minLowEHitTime_ > time) minLowEHitTime_ = time;
00057         if(maxLowEHitTime_ < time) maxLowEHitTime_ = time;
00058         lowEHitTimeSqrd_ += time*time;
00059         ++numLowEHits_;
00060       }
00061       if(energy>=minHighHitE) {
00062         if(minHighEHitTime_ > time) minHighEHitTime_ = time;
00063         if(maxHighEHitTime_ < time) maxHighEHitTime_ = time;
00064         highEHitTimeSqrd_ += time*time;
00065         ++numHighEHits_;
00066       }
00067     }
00068   }
00069 
00070   // emf
00071   HPDEMF_ = 999.;
00072   for(std::vector<reco::HcalNoiseHPD>::const_iterator it1=rbx.HPDsBegin(); it1!=rbx.HPDsEnd(); ++it1) {
00073     double eme=it1->caloTowerEmE();
00074     double hade=it1->recHitEnergy(minRecHitE);
00075     double emf=(eme+hade)==0 ? 999 : eme/(eme+hade);
00076     if(HPDEMF_ > emf) emf = HPDEMF_;
00077   }
00078   double eme=rbx.caloTowerEmE();
00079   RBXEMF_ = (eme+energy_)==0 ? 999 : eme/(eme+energy_);
00080 
00081   // calotowers
00082   rbxtowers_.clear();
00083   JoinCaloTowerRefVectorsWithoutDuplicates join;
00084   for(std::vector<reco::HcalNoiseHPD>::const_iterator it1=rbx.HPDsBegin(); it1!=rbx.HPDsEnd(); ++it1) {
00085     join(rbxtowers_, it1->caloTowers());
00086   }
00087 
00088   return;
00089 }
00090 
00091 HcalNoiseAlgo::HcalNoiseAlgo(const edm::ParameterSet& iConfig)
00092 {
00093   pMinERatio_ = iConfig.getParameter<double>("pMinERatio");
00094   pMinEZeros_ = iConfig.getParameter<double>("pMinEZeros");
00095   pMinEEMF_ = iConfig.getParameter<double>("pMinEEMF");
00096 
00097   minERatio_ = iConfig.getParameter<double>("minERatio");
00098   minEZeros_ = iConfig.getParameter<double>("minEZeros");
00099   minEEMF_ = iConfig.getParameter<double>("minEEMF");  
00100 
00101   pMinE_ = iConfig.getParameter<double>("pMinE");
00102   pMinRatio_ = iConfig.getParameter<double>("pMinRatio");
00103   pMaxRatio_ = iConfig.getParameter<double>("pMaxRatio");
00104   pMinHPDHits_ = iConfig.getParameter<int>("pMinHPDHits");
00105   pMinRBXHits_ = iConfig.getParameter<int>("pMinRBXHits");
00106   pMinHPDNoOtherHits_ = iConfig.getParameter<int>("pMinHPDNoOtherHits");
00107   pMinZeros_ = iConfig.getParameter<int>("pMinZeros");
00108   pMinLowEHitTime_ = iConfig.getParameter<double>("pMinLowEHitTime");
00109   pMaxLowEHitTime_ = iConfig.getParameter<double>("pMaxLowEHitTime");
00110   pMinHighEHitTime_ = iConfig.getParameter<double>("pMinHighEHitTime");
00111   pMaxHighEHitTime_ = iConfig.getParameter<double>("pMaxHighEHitTime");
00112   pMaxHPDEMF_ = iConfig.getParameter<double>("pMaxHPDEMF");
00113   pMaxRBXEMF_ = iConfig.getParameter<double>("pMaxRBXEMF");
00114 
00115   lMinRatio_ = iConfig.getParameter<double>("lMinRatio");
00116   lMaxRatio_ = iConfig.getParameter<double>("lMaxRatio");
00117   lMinHPDHits_ = iConfig.getParameter<int>("lMinHPDHits");
00118   lMinRBXHits_ = iConfig.getParameter<int>("lMinRBXHits");
00119   lMinHPDNoOtherHits_ = iConfig.getParameter<int>("lMinHPDNoOtherHits");
00120   lMinZeros_ = iConfig.getParameter<int>("lMinZeros");
00121   lMinLowEHitTime_ = iConfig.getParameter<double>("lMinLowEHitTime");
00122   lMaxLowEHitTime_ = iConfig.getParameter<double>("lMaxLowEHitTime");
00123   lMinHighEHitTime_ = iConfig.getParameter<double>("lMinHighEHitTime");
00124   lMaxHighEHitTime_ = iConfig.getParameter<double>("lMaxHighEHitTime");
00125 
00126   tMinRatio_ = iConfig.getParameter<double>("tMinRatio");
00127   tMaxRatio_ = iConfig.getParameter<double>("tMaxRatio");
00128   tMinHPDHits_ = iConfig.getParameter<int>("tMinHPDHits");
00129   tMinRBXHits_ = iConfig.getParameter<int>("tMinRBXHits");
00130   tMinHPDNoOtherHits_ = iConfig.getParameter<int>("tMinHPDNoOtherHits");
00131   tMinZeros_ = iConfig.getParameter<int>("tMinZeros");
00132   tMinLowEHitTime_ = iConfig.getParameter<double>("tMinLowEHitTime");
00133   tMaxLowEHitTime_ = iConfig.getParameter<double>("tMaxLowEHitTime");
00134   tMinHighEHitTime_ = iConfig.getParameter<double>("tMinHighEHitTime");
00135   tMaxHighEHitTime_ = iConfig.getParameter<double>("tMaxHighEHitTime");
00136 
00137   hlMaxHPDEMF_ = iConfig.getParameter<double>("hlMaxHPDEMF");
00138   hlMaxRBXEMF_ = iConfig.getParameter<double>("hlMaxRBXEMF");
00139 }
00140 
00141 bool HcalNoiseAlgo::isProblematic(const CommonHcalNoiseRBXData& data) const
00142 {
00143   if(data.energy()>pMinE_) return true;
00144   if(data.validRatio() && data.energy()>pMinERatio_ && data.ratio()<pMinRatio_) return true;
00145   if(data.validRatio() && data.energy()>pMinERatio_ && data.ratio()>pMaxRatio_) return true;
00146   if(data.numHPDHits()>=pMinHPDHits_) return true;
00147   if(data.numRBXHits()>=pMinRBXHits_) return true;
00148   if(data.numHPDNoOtherHits()>=pMinHPDNoOtherHits_) return true;
00149   if(data.numZeros()>=pMinZeros_ && data.energy()>pMinEZeros_) return true;
00150   if(data.minLowEHitTime()<pMinLowEHitTime_) return true;
00151   if(data.maxLowEHitTime()>pMaxLowEHitTime_) return true;
00152   if(data.minHighEHitTime()<pMinHighEHitTime_) return true;
00153   if(data.maxHighEHitTime()>pMaxHighEHitTime_) return true;
00154   if(data.HPDEMF()<pMaxHPDEMF_ && data.energy()>pMinEEMF_) return true;
00155   if(data.RBXEMF()<pMaxRBXEMF_ && data.energy()>pMinEEMF_) return true;  return false;
00156 }
00157 
00158 
00159 bool HcalNoiseAlgo::passLooseNoiseFilter(const CommonHcalNoiseRBXData& data) const
00160 {
00161   return (passLooseRatio(data) && passLooseHits(data) && passLooseZeros(data) && passLooseTiming(data));
00162 }
00163 
00164 bool HcalNoiseAlgo::passTightNoiseFilter(const CommonHcalNoiseRBXData& data) const
00165 {
00166   return (passTightRatio(data) && passTightHits(data) && passTightZeros(data) && passTightTiming(data));
00167 }
00168 
00169 bool HcalNoiseAlgo::passHighLevelNoiseFilter(const CommonHcalNoiseRBXData& data) const
00170 {
00171   if(passEMFThreshold(data)) {
00172     if(data.HPDEMF()<hlMaxHPDEMF_) return false;
00173     if(data.RBXEMF()<hlMaxRBXEMF_) return false;
00174   }
00175   return true;
00176 }
00177 
00178 bool HcalNoiseAlgo::passLooseRatio(const CommonHcalNoiseRBXData& data) const
00179 {
00180   if(passRatioThreshold(data)) {
00181     if(data.validRatio() && data.ratio()<lMinRatio_) return false;
00182     if(data.validRatio() && data.ratio()>lMaxRatio_) return false;
00183   }
00184   return true;
00185 }
00186 
00187 bool HcalNoiseAlgo::passLooseHits(const CommonHcalNoiseRBXData& data) const
00188 {
00189   if(data.numHPDHits()>=lMinHPDHits_) return false;
00190   if(data.numRBXHits()>=lMinRBXHits_) return false;
00191   if(data.numHPDNoOtherHits()>=lMinHPDNoOtherHits_) return false;
00192   return true;
00193 }
00194 
00195 bool HcalNoiseAlgo::passLooseZeros(const CommonHcalNoiseRBXData& data) const
00196 {
00197   if(passZerosThreshold(data)) {
00198     if(data.numZeros()>=lMinZeros_) return false;
00199   }
00200   return true;
00201 }
00202 
00203 bool HcalNoiseAlgo::passLooseTiming(const CommonHcalNoiseRBXData& data) const
00204 {
00205   if(data.minLowEHitTime()<lMinLowEHitTime_) return false;
00206   if(data.maxLowEHitTime()>lMaxLowEHitTime_) return false;
00207   if(data.minHighEHitTime()<lMinHighEHitTime_) return false;
00208   if(data.maxHighEHitTime()>lMaxHighEHitTime_) return false;
00209   return true;
00210 }
00211 
00212 bool HcalNoiseAlgo::passTightRatio(const CommonHcalNoiseRBXData& data) const
00213 {
00214   if(passRatioThreshold(data)) {
00215     if(data.validRatio() && data.ratio()<tMinRatio_) return false;
00216     if(data.validRatio() && data.ratio()>tMaxRatio_) return false;
00217   }
00218   return true;
00219 }
00220 
00221 bool HcalNoiseAlgo::passTightHits(const CommonHcalNoiseRBXData& data) const
00222 {
00223   if(data.numHPDHits()>=tMinHPDHits_) return false;
00224   if(data.numRBXHits()>=tMinRBXHits_) return false;
00225   if(data.numHPDNoOtherHits()>=tMinHPDNoOtherHits_) return false;
00226   return true;
00227 }
00228 
00229 bool HcalNoiseAlgo::passTightZeros(const CommonHcalNoiseRBXData& data) const
00230 {
00231   if(passZerosThreshold(data)) {
00232     if(data.numZeros()>=tMinZeros_) return false;
00233   }
00234   return true;
00235 }
00236 
00237 bool HcalNoiseAlgo::passTightTiming(const CommonHcalNoiseRBXData& data) const
00238 {
00239   if(data.minLowEHitTime()<tMinLowEHitTime_) return false;
00240   if(data.maxLowEHitTime()>tMaxLowEHitTime_) return false;
00241   if(data.minHighEHitTime()<tMinHighEHitTime_) return false;
00242   if(data.maxHighEHitTime()>tMaxHighEHitTime_) return false;
00243   return true;
00244 }
00245 
00246 bool HcalNoiseAlgo::passRatioThreshold(const CommonHcalNoiseRBXData& data) const
00247 {
00248   return (data.energy()>minERatio_);
00249 }
00250 
00251 bool HcalNoiseAlgo::passZerosThreshold(const CommonHcalNoiseRBXData& data) const
00252 {
00253   return (data.energy()>minEZeros_);
00254 }
00255 
00256 bool HcalNoiseAlgo::passEMFThreshold(const CommonHcalNoiseRBXData& data) const
00257 {
00258   return (data.energy()>minEEMF_);
00259 }
00260 
00261 void JoinCaloTowerRefVectorsWithoutDuplicates::operator()(edm::RefVector<CaloTowerCollection>& v1, const edm::RefVector<CaloTowerCollection>& v2) const
00262 {
00263   // combines them first into a set to get rid of duplicates and then puts them into the first vector
00264 
00265   // sorts them first to get rid of duplicates, then puts them into another RefVector
00266   twrrefset_t twrrefset;
00267   for(edm::RefVector<CaloTowerCollection>::const_iterator it=v1.begin(); it!=v1.end(); ++it)
00268     twrrefset.insert(*it);
00269   for(edm::RefVector<CaloTowerCollection>::const_iterator it=v2.begin(); it!=v2.end(); ++it)
00270     twrrefset.insert(*it);
00271 
00272   // clear the original refvector and put them back in
00273   v1.clear();
00274   for(twrrefset_t::const_iterator it=twrrefset.begin(); it!=twrrefset.end(); ++it) {
00275     v1.push_back(*it);
00276   }
00277   return;
00278 }
00279 
00280 bool CommonHcalNoiseRBXData::CheckPassFilter(double Charge, double Discriminant,
00281    std::vector<std::pair<double, double> > &Cuts, int Side)
00282 {
00283    //
00284    // Checks whether Discriminant value passes Cuts for the specified Charge.  True if pulse is good.
00285    //
00286    // The "Cuts" pairs are assumed to be sorted in terms of size from small to large,
00287    //    where each "pair" = (Charge, Discriminant)
00288    // "Side" is either positive or negative, which determines whether to discard the pulse if discriminant
00289    //    is greater or smaller than the cut value
00290    //
00291 
00292    if(Cuts.size() == 0)   // safety check that there are some cuts defined
00293       return true;
00294 
00295    if(Charge <= Cuts[0].first)   // too small to cut on
00296       return true;
00297 
00298    int IndexLargerThanCharge = -1;   // find the range it is falling in
00299    for(int i = 1; i < (int)Cuts.size(); i++)
00300    {
00301       if(Cuts[i].first > Charge)
00302       {
00303          IndexLargerThanCharge = i;
00304          break;
00305       }
00306    }
00307 
00308    double limit = 1000000;
00309 
00310    if(IndexLargerThanCharge == -1)   // if charge is greater than the last entry, assume flat line
00311       limit = Cuts[Cuts.size()-1].second;
00312    else   // otherwise, do a linear interpolation to find the cut position
00313    {
00314       double C1 = Cuts[IndexLargerThanCharge].first;
00315       double C2 = Cuts[IndexLargerThanCharge-1].first;
00316       double L1 = Cuts[IndexLargerThanCharge].second;
00317       double L2 = Cuts[IndexLargerThanCharge-1].second;
00318 
00319       limit = (Charge - C1) / (C2 - C1) * (L2 - L1) + L1;
00320    }
00321 
00322    if(Side > 0 && Discriminant > limit)
00323       return false;
00324    if(Side < 0 && Discriminant < limit)
00325       return false;
00326 
00327    return true;
00328 }
00329
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoMET/METAlgorithms/src/HcalNoiseAlgo.cc
