Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:41 

0001 #include "RecoMET/METAlgorithms/interface/HcalNoiseAlgo.h"
0002 
0003 CommonHcalNoiseRBXData::CommonHcalNoiseRBXData(const reco::HcalNoiseRBX& rbx,
0004                                                double minRecHitE,
0005                                                double minLowHitE,
0006                                                double minHighHitE,
0007                                                double TS4TS5EnergyThreshold,
0008                                                std::vector<std::pair<double, double> > const& TS4TS5UpperCut,
0009                                                std::vector<std::pair<double, double> > const& TS4TS5LowerCut,
0010                                                double minRBXRechitR45E)
0011     : r45Count_(0), r45Fraction_(0), r45EnergyFraction_(0) {
0012   // energy
0013   energy_ = rbx.recHitEnergy(minRecHitE);
0014 
0015   // ratio
0016   e2ts_ = rbx.allChargeHighest2TS();
0017   e10ts_ = rbx.allChargeTotal();
0018 
0019   // TS4TS5
0020   TS4TS5Decision_ = true;
0021   if (energy_ > TS4TS5EnergyThreshold)  // check filter
0022   {
0023     std::vector<float> AllCharge = rbx.allCharge();
0024     double BaseCharge = AllCharge[4] + AllCharge[5];
0025     if (BaseCharge < 1)
0026       BaseCharge = 1;
0027     double TS4TS5 = (AllCharge[4] - AllCharge[5]) / BaseCharge;
0028 
0029     if (CheckPassFilter(BaseCharge, TS4TS5, TS4TS5UpperCut, 1) == false)
0030       TS4TS5Decision_ = false;
0031     if (CheckPassFilter(BaseCharge, TS4TS5, TS4TS5LowerCut, -1) == false)
0032       TS4TS5Decision_ = false;
0033   } else
0034     TS4TS5Decision_ = true;
0035 
0036   // Rechit-wide R45
0037   int rbxHitCount = rbx.numRecHits(minRBXRechitR45E);
0038   r45Count_ = 0;
0039   r45Fraction_ = 0;
0040   r45EnergyFraction_ = 0;
0041   if (rbxHitCount > 0) {
0042     r45Count_ = rbx.numRecHitsFailR45(minRBXRechitR45E);
0043     r45Fraction_ = (double)(r45Count_) / (double)(rbxHitCount);
0044     r45EnergyFraction_ = rbx.recHitEnergyFailR45(minRBXRechitR45E) / rbx.recHitEnergy(minRBXRechitR45E);
0045   }
0046 
0047   // # of hits
0048   numHPDHits_ = 0;
0049   for (std::vector<reco::HcalNoiseHPD>::const_iterator it1 = rbx.HPDsBegin(); it1 != rbx.HPDsEnd(); ++it1) {
0050     int nhpdhits = it1->numRecHits(minRecHitE);
0051     if (numHPDHits_ < nhpdhits)
0052       numHPDHits_ = nhpdhits;
0053   }
0054   numRBXHits_ = rbx.numRecHits(minRecHitE);
0055   numHPDNoOtherHits_ = (numHPDHits_ == numRBXHits_) ? numHPDHits_ : 0;
0056 
0057   // # of ADC zeros
0058   numZeros_ = rbx.totalZeros();
0059 
0060   // timing
0061   minLowEHitTime_ = minHighEHitTime_ = 99999.;
0062   maxLowEHitTime_ = maxHighEHitTime_ = -99999.;
0063   lowEHitTimeSqrd_ = highEHitTimeSqrd_ = 0;
0064   numLowEHits_ = numHighEHits_ = 0;
0065   for (std::vector<reco::HcalNoiseHPD>::const_iterator it1 = rbx.HPDsBegin(); it1 != rbx.HPDsEnd(); ++it1) {
0066     edm::RefVector<HBHERecHitCollection> rechits = it1->recHits();
0067     for (edm::RefVector<HBHERecHitCollection>::const_iterator it2 = rechits.begin(); it2 != rechits.end(); ++it2) {
0068       float energy = (*it2)->energy();
0069       float time = (*it2)->time();
0070       if (energy >= minLowHitE) {
0071         if (minLowEHitTime_ > time)
0072           minLowEHitTime_ = time;
0073         if (maxLowEHitTime_ < time)
0074           maxLowEHitTime_ = time;
0075         lowEHitTimeSqrd_ += time * time;
0076         ++numLowEHits_;
0077       }
0078       if (energy >= minHighHitE) {
0079         if (minHighEHitTime_ > time)
0080           minHighEHitTime_ = time;
0081         if (maxHighEHitTime_ < time)
0082           maxHighEHitTime_ = time;
0083         highEHitTimeSqrd_ += time * time;
0084         ++numHighEHits_;
0085       }
0086     }
0087   }
0088 
0089   // emf (get the one with minimum value)
0090   HPDEMF_ = 999.;
0091   for (std::vector<reco::HcalNoiseHPD>::const_iterator it1 = rbx.HPDsBegin(); it1 != rbx.HPDsEnd(); ++it1) {
0092     double eme = it1->caloTowerEmE();
0093     double hade = it1->recHitEnergy(minRecHitE);
0094     double emf = (eme + hade) == 0 ? 999 : eme / (eme + hade);
0095     if (HPDEMF_ > emf)
0096       HPDEMF_ = emf;
0097   }
0098   double eme = rbx.caloTowerEmE();
0099   RBXEMF_ = (eme + energy_) == 0 ? 999 : eme / (eme + energy_);
0100 
0101   // calotowers
0102   rbxtowers_.clear();
0103   JoinCaloTowerRefVectorsWithoutDuplicates join;
0104   for (std::vector<reco::HcalNoiseHPD>::const_iterator it1 = rbx.HPDsBegin(); it1 != rbx.HPDsEnd(); ++it1) {
0105     join(rbxtowers_, it1->caloTowers());
0106   }
0107 
0108   return;
0109 }
0110 
0111 HcalNoiseAlgo::HcalNoiseAlgo(const edm::ParameterSet& iConfig) {
0112   pMinERatio_ = iConfig.getParameter<double>("pMinERatio");
0113   pMinEZeros_ = iConfig.getParameter<double>("pMinEZeros");
0114   pMinEEMF_ = iConfig.getParameter<double>("pMinEEMF");
0115 
0116   minERatio_ = iConfig.getParameter<double>("minERatio");
0117   minEZeros_ = iConfig.getParameter<double>("minEZeros");
0118   minEEMF_ = iConfig.getParameter<double>("minEEMF");
0119 
0120   pMinE_ = iConfig.getParameter<double>("pMinE");
0121   pMinRatio_ = iConfig.getParameter<double>("pMinRatio");
0122   pMaxRatio_ = iConfig.getParameter<double>("pMaxRatio");
0123   pMinHPDHits_ = iConfig.getParameter<int>("pMinHPDHits");
0124   pMinRBXHits_ = iConfig.getParameter<int>("pMinRBXHits");
0125   pMinHPDNoOtherHits_ = iConfig.getParameter<int>("pMinHPDNoOtherHits");
0126   pMinZeros_ = iConfig.getParameter<int>("pMinZeros");
0127   pMinLowEHitTime_ = iConfig.getParameter<double>("pMinLowEHitTime");
0128   pMaxLowEHitTime_ = iConfig.getParameter<double>("pMaxLowEHitTime");
0129   pMinHighEHitTime_ = iConfig.getParameter<double>("pMinHighEHitTime");
0130   pMaxHighEHitTime_ = iConfig.getParameter<double>("pMaxHighEHitTime");
0131   pMaxHPDEMF_ = iConfig.getParameter<double>("pMaxHPDEMF");
0132   pMaxRBXEMF_ = iConfig.getParameter<double>("pMaxRBXEMF");
0133 
0134   if (iConfig.existsAs<int>("pMinRBXRechitR45Count"))
0135     pMinRBXRechitR45Count_ = iConfig.getParameter<int>("pMinRBXRechitR45Count");
0136   else
0137     pMinRBXRechitR45Count_ = 0;
0138   if (iConfig.existsAs<double>("pMinRBXRechitR45Fraction"))
0139     pMinRBXRechitR45Fraction_ = iConfig.getParameter<double>("pMinRBXRechitR45Fraction");
0140   else
0141     pMinRBXRechitR45Fraction_ = 0;
0142   if (iConfig.existsAs<double>("pMinRBXRechitR45EnergyFraction"))
0143     pMinRBXRechitR45EnergyFraction_ = iConfig.getParameter<double>("pMinRBXRechitR45EnergyFraction");
0144   else
0145     pMinRBXRechitR45EnergyFraction_ = 0;
0146 
0147   lMinRatio_ = iConfig.getParameter<double>("lMinRatio");
0148   lMaxRatio_ = iConfig.getParameter<double>("lMaxRatio");
0149   lMinHPDHits_ = iConfig.getParameter<int>("lMinHPDHits");
0150   lMinRBXHits_ = iConfig.getParameter<int>("lMinRBXHits");
0151   lMinHPDNoOtherHits_ = iConfig.getParameter<int>("lMinHPDNoOtherHits");
0152   lMinZeros_ = iConfig.getParameter<int>("lMinZeros");
0153   lMinLowEHitTime_ = iConfig.getParameter<double>("lMinLowEHitTime");
0154   lMaxLowEHitTime_ = iConfig.getParameter<double>("lMaxLowEHitTime");
0155   lMinHighEHitTime_ = iConfig.getParameter<double>("lMinHighEHitTime");
0156   lMaxHighEHitTime_ = iConfig.getParameter<double>("lMaxHighEHitTime");
0157 
0158   if (iConfig.existsAs<std::vector<double> >("lRBXRecHitR45Cuts"))
0159     lMinRBXRechitR45Cuts_ = iConfig.getParameter<std::vector<double> >("lRBXRecHitR45Cuts");
0160   else {
0161     double defaultCut[4] = {-999, -999, -999, -999};
0162     lMinRBXRechitR45Cuts_ = std::vector<double>(defaultCut, defaultCut + 4);
0163   }
0164 
0165   tMinRatio_ = iConfig.getParameter<double>("tMinRatio");
0166   tMaxRatio_ = iConfig.getParameter<double>("tMaxRatio");
0167   tMinHPDHits_ = iConfig.getParameter<int>("tMinHPDHits");
0168   tMinRBXHits_ = iConfig.getParameter<int>("tMinRBXHits");
0169   tMinHPDNoOtherHits_ = iConfig.getParameter<int>("tMinHPDNoOtherHits");
0170   tMinZeros_ = iConfig.getParameter<int>("tMinZeros");
0171   tMinLowEHitTime_ = iConfig.getParameter<double>("tMinLowEHitTime");
0172   tMaxLowEHitTime_ = iConfig.getParameter<double>("tMaxLowEHitTime");
0173   tMinHighEHitTime_ = iConfig.getParameter<double>("tMinHighEHitTime");
0174   tMaxHighEHitTime_ = iConfig.getParameter<double>("tMaxHighEHitTime");
0175 
0176   if (iConfig.existsAs<std::vector<double> >("tRBXRecHitR45Cuts"))
0177     tMinRBXRechitR45Cuts_ = iConfig.getParameter<std::vector<double> >("tRBXRecHitR45Cuts");
0178   else {
0179     double defaultCut[4] = {-999, -999, -999, -999};
0180     tMinRBXRechitR45Cuts_ = std::vector<double>(defaultCut, defaultCut + 4);
0181   }
0182 
0183   hlMaxHPDEMF_ = iConfig.getParameter<double>("hlMaxHPDEMF");
0184   hlMaxRBXEMF_ = iConfig.getParameter<double>("hlMaxRBXEMF");
0185 }
0186 
0187 bool HcalNoiseAlgo::isProblematic(const CommonHcalNoiseRBXData& data) const {
0188   if (data.energy() > pMinE_)
0189     return true;
0190   if (data.validRatio() && data.energy() > pMinERatio_ && data.ratio() < pMinRatio_)
0191     return true;
0192   if (data.validRatio() && data.energy() > pMinERatio_ && data.ratio() > pMaxRatio_)
0193     return true;
0194   if (data.numHPDHits() >= pMinHPDHits_)
0195     return true;
0196   if (data.numRBXHits() >= pMinRBXHits_)
0197     return true;
0198   if (data.numHPDNoOtherHits() >= pMinHPDNoOtherHits_)
0199     return true;
0200   if (data.numZeros() >= pMinZeros_ && data.energy() > pMinEZeros_)
0201     return true;
0202   if (data.minLowEHitTime() < pMinLowEHitTime_)
0203     return true;
0204   if (data.maxLowEHitTime() > pMaxLowEHitTime_)
0205     return true;
0206   if (data.minHighEHitTime() < pMinHighEHitTime_)
0207     return true;
0208   if (data.maxHighEHitTime() > pMaxHighEHitTime_)
0209     return true;
0210   if (data.HPDEMF() < pMaxHPDEMF_ && data.energy() > pMinEEMF_)
0211     return true;
0212   if (data.RBXEMF() < pMaxRBXEMF_ && data.energy() > pMinEEMF_)
0213     return true;
0214   if (data.r45Count() >= pMinRBXRechitR45Count_)
0215     return true;
0216   if (data.r45Fraction() >= pMinRBXRechitR45Fraction_)
0217     return true;
0218   if (data.r45EnergyFraction() >= pMinRBXRechitR45EnergyFraction_)
0219     return true;
0220 
0221   return false;
0222 }
0223 
0224 bool HcalNoiseAlgo::passLooseNoiseFilter(const CommonHcalNoiseRBXData& data) const {
0225   return (passLooseRatio(data) && passLooseHits(data) && passLooseZeros(data) && passLooseTiming(data) &&
0226           passLooseRBXRechitR45(data));
0227 }
0228 
0229 bool HcalNoiseAlgo::passTightNoiseFilter(const CommonHcalNoiseRBXData& data) const {
0230   return (passTightRatio(data) && passTightHits(data) && passTightZeros(data) && passTightTiming(data) &&
0231           passTightRBXRechitR45(data));
0232 }
0233 
0234 bool HcalNoiseAlgo::passHighLevelNoiseFilter(const CommonHcalNoiseRBXData& data) const {
0235   if (passEMFThreshold(data)) {
0236     if (data.HPDEMF() < hlMaxHPDEMF_)
0237       return false;
0238     if (data.RBXEMF() < hlMaxRBXEMF_)
0239       return false;
0240   }
0241   return true;
0242 }
0243 
0244 bool HcalNoiseAlgo::passLooseRatio(const CommonHcalNoiseRBXData& data) const {
0245   if (passRatioThreshold(data)) {
0246     if (data.validRatio() && data.ratio() < lMinRatio_)
0247       return false;
0248     if (data.validRatio() && data.ratio() > lMaxRatio_)
0249       return false;
0250   }
0251   return true;
0252 }
0253 
0254 bool HcalNoiseAlgo::passLooseHits(const CommonHcalNoiseRBXData& data) const {
0255   if (data.numHPDHits() >= lMinHPDHits_)
0256     return false;
0257   if (data.numRBXHits() >= lMinRBXHits_)
0258     return false;
0259   if (data.numHPDNoOtherHits() >= lMinHPDNoOtherHits_)
0260     return false;
0261   return true;
0262 }
0263 
0264 bool HcalNoiseAlgo::passLooseZeros(const CommonHcalNoiseRBXData& data) const {
0265   if (passZerosThreshold(data)) {
0266     if (data.numZeros() >= lMinZeros_)
0267       return false;
0268   }
0269   return true;
0270 }
0271 
0272 bool HcalNoiseAlgo::passLooseTiming(const CommonHcalNoiseRBXData& data) const {
0273   if (data.minLowEHitTime() < lMinLowEHitTime_)
0274     return false;
0275   if (data.maxLowEHitTime() > lMaxLowEHitTime_)
0276     return false;
0277   if (data.minHighEHitTime() < lMinHighEHitTime_)
0278     return false;
0279   if (data.maxHighEHitTime() > lMaxHighEHitTime_)
0280     return false;
0281   return true;
0282 }
0283 
0284 bool HcalNoiseAlgo::passLooseRBXRechitR45(const CommonHcalNoiseRBXData& data) const {
0285   int Count = data.r45Count();
0286   double Fraction = data.r45Fraction();
0287   double EnergyFraction = data.r45EnergyFraction();
0288 
0289   for (int i = 0; i + 3 < (int)lMinRBXRechitR45Cuts_.size(); i = i + 4) {
0290     double Value = Count * lMinRBXRechitR45Cuts_[i] + Fraction * lMinRBXRechitR45Cuts_[i + 1] +
0291                    EnergyFraction * lMinRBXRechitR45Cuts_[i + 2] + lMinRBXRechitR45Cuts_[i + 3];
0292     if (Value > 0)
0293       return false;
0294   }
0295 
0296   return true;
0297 }
0298 
0299 bool HcalNoiseAlgo::passTightRatio(const CommonHcalNoiseRBXData& data) const {
0300   if (passRatioThreshold(data)) {
0301     if (data.validRatio() && data.ratio() < tMinRatio_)
0302       return false;
0303     if (data.validRatio() && data.ratio() > tMaxRatio_)
0304       return false;
0305   }
0306   return true;
0307 }
0308 
0309 bool HcalNoiseAlgo::passTightHits(const CommonHcalNoiseRBXData& data) const {
0310   if (data.numHPDHits() >= tMinHPDHits_)
0311     return false;
0312   if (data.numRBXHits() >= tMinRBXHits_)
0313     return false;
0314   if (data.numHPDNoOtherHits() >= tMinHPDNoOtherHits_)
0315     return false;
0316   return true;
0317 }
0318 
0319 bool HcalNoiseAlgo::passTightZeros(const CommonHcalNoiseRBXData& data) const {
0320   if (passZerosThreshold(data)) {
0321     if (data.numZeros() >= tMinZeros_)
0322       return false;
0323   }
0324   return true;
0325 }
0326 
0327 bool HcalNoiseAlgo::passTightTiming(const CommonHcalNoiseRBXData& data) const {
0328   if (data.minLowEHitTime() < tMinLowEHitTime_)
0329     return false;
0330   if (data.maxLowEHitTime() > tMaxLowEHitTime_)
0331     return false;
0332   if (data.minHighEHitTime() < tMinHighEHitTime_)
0333     return false;
0334   if (data.maxHighEHitTime() > tMaxHighEHitTime_)
0335     return false;
0336   return true;
0337 }
0338 
0339 bool HcalNoiseAlgo::passTightRBXRechitR45(const CommonHcalNoiseRBXData& data) const {
0340   int Count = data.r45Count();
0341   double Fraction = data.r45Fraction();
0342   double EnergyFraction = data.r45EnergyFraction();
0343 
0344   for (int i = 0; i + 3 < (int)tMinRBXRechitR45Cuts_.size(); i = i + 4) {
0345     double Value = Count * tMinRBXRechitR45Cuts_[i] + Fraction * tMinRBXRechitR45Cuts_[i + 1] +
0346                    EnergyFraction * tMinRBXRechitR45Cuts_[i + 2] + tMinRBXRechitR45Cuts_[i + 3];
0347     if (Value > 0)
0348       return false;
0349   }
0350 
0351   return true;
0352 }
0353 
0354 bool HcalNoiseAlgo::passRatioThreshold(const CommonHcalNoiseRBXData& data) const {
0355   return (data.energy() > minERatio_);
0356 }
0357 
0358 bool HcalNoiseAlgo::passZerosThreshold(const CommonHcalNoiseRBXData& data) const {
0359   return (data.energy() > minEZeros_);
0360 }
0361 
0362 bool HcalNoiseAlgo::passEMFThreshold(const CommonHcalNoiseRBXData& data) const { return (data.energy() > minEEMF_); }
0363 
0364 void JoinCaloTowerRefVectorsWithoutDuplicates::operator()(edm::RefVector<CaloTowerCollection>& v1,
0365                                                           const edm::RefVector<CaloTowerCollection>& v2) const {
0366   // combines them first into a set to get rid of duplicates and then puts them into the first vector
0367 
0368   // sorts them first to get rid of duplicates, then puts them into another RefVector
0369   twrrefset_t twrrefset;
0370   for (edm::RefVector<CaloTowerCollection>::const_iterator it = v1.begin(); it != v1.end(); ++it)
0371     twrrefset.insert(*it);
0372   for (edm::RefVector<CaloTowerCollection>::const_iterator it = v2.begin(); it != v2.end(); ++it)
0373     twrrefset.insert(*it);
0374 
0375   // clear the original refvector and put them back in
0376   v1.clear();
0377   for (twrrefset_t::const_iterator it = twrrefset.begin(); it != twrrefset.end(); ++it) {
0378     v1.push_back(*it);
0379   }
0380   return;
0381 }
0382 
0383 bool CommonHcalNoiseRBXData::CheckPassFilter(double Charge,
0384                                              double Discriminant,
0385                                              std::vector<std::pair<double, double> > const& Cuts,
0386                                              int Side) {
0387   //
0388   // Checks whether Discriminant value passes Cuts for the specified Charge.  True if pulse is good.
0389   //
0390   // The "Cuts" pairs are assumed to be sorted in terms of size from small to large,
0391   //    where each "pair" = (Charge, Discriminant)
0392   // "Side" is either positive or negative, which determines whether to discard the pulse if discriminant
0393   //    is greater or smaller than the cut value
0394   //
0395 
0396   if (Cuts.empty())  // safety check that there are some cuts defined
0397     return true;
0398 
0399   if (Charge <= Cuts[0].first)  // too small to cut on
0400     return true;
0401 
0402   int IndexLargerThanCharge = -1;  // find the range it is falling in
0403   for (int i = 1; i < (int)Cuts.size(); i++) {
0404     if (Cuts[i].first > Charge) {
0405       IndexLargerThanCharge = i;
0406       break;
0407     }
0408   }
0409 
0410   double limit = 1000000;
0411 
0412   if (IndexLargerThanCharge == -1)  // if charge is greater than the last entry, assume flat line
0413     limit = Cuts[Cuts.size() - 1].second;
0414   else  // otherwise, do a linear interpolation to find the cut position
0415   {
0416     double C1 = Cuts[IndexLargerThanCharge].first;
0417     double C2 = Cuts[IndexLargerThanCharge - 1].first;
0418     double L1 = Cuts[IndexLargerThanCharge].second;
0419     double L2 = Cuts[IndexLargerThanCharge - 1].second;
0420 
0421     limit = (Charge - C1) / (C2 - C1) * (L2 - L1) + L1;
0422   }
0423 
0424   if (Side > 0 && Discriminant > limit)
0425     return false;
0426   if (Side < 0 && Discriminant < limit)
0427     return false;
0428 
0429   return true;
0430 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!