EcalRecAlgos/src/EcalCleaningAlgo.cc

0001 /* Implementation of class EcalCleaningAlgo
0002    \author Stefano Argiro
0003    \date 20 Dec 2010
0004 */
0005
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0008 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0009 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0010
0011 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalCleaningAlgo.h"
0012
0013 EcalCleaningAlgo::EcalCleaningAlgo(const edm::ParameterSet& p) {
0014   cThreshold_barrel_ = p.getParameter<double>("cThreshold_barrel");
0015   ;
0016   cThreshold_endcap_ = p.getParameter<double>("cThreshold_endcap");
0017   ;
0018   e4e1_a_barrel_ = p.getParameter<double>("e4e1_a_barrel");
0019   e4e1_b_barrel_ = p.getParameter<double>("e4e1_b_barrel");
0020   e4e1_a_endcap_ = p.getParameter<double>("e4e1_a_endcap");
0021   e4e1_b_endcap_ = p.getParameter<double>("e4e1_b_endcap");
0022   e4e1Treshold_barrel_ = p.getParameter<double>("e4e1Threshold_barrel");
0023   e4e1Treshold_endcap_ = p.getParameter<double>("e4e1Threshold_endcap");
0024
0025   cThreshold_double_ = p.getParameter<double>("cThreshold_double");
0026
0027   ignoreOutOfTimeThresh_ = p.getParameter<double>("ignoreOutOfTimeThresh");
0028   tightenCrack_e1_single_ = p.getParameter<double>("tightenCrack_e1_single");
0029   tightenCrack_e4e1_single_ = p.getParameter<double>("tightenCrack_e4e1_single");
0030   tightenCrack_e1_double_ = p.getParameter<double>("tightenCrack_e1_double");
0031   tightenCrack_e6e2_double_ = p.getParameter<double>("tightenCrack_e6e2_double");
0032   e6e2thresh_ = p.getParameter<double>("e6e2thresh");
0033 }
0034
0035 /**
0036    Flag spikey channels
0037
0038    Mark single spikes. Spike definition:
0039
0040       Barrel: e> cThreshold_barrel_  &&
0041               e4e1 > e4e1_a_barrel_ * log10(e) + e4e1_b_barrel_
0042
0043       Near cracks: energy threshold is multiplied by tightenCrack_e1_single
0044                    e4e1 threshold is divided by tightenCrack_e4e1_single
0045
0046       Endcap : e> cThreshold_endcap_ &&
0047                e4e1>   e4e1_a_endcap_ * log10(e) + e4e1_b_endcap_
0048
0049    Mark double spikes    (barrel only)
0050       e> cThreshold_double_ &&
0051       e6e2 >  e6e2thresh_;
0052
0053    Near cracks:
0054           energy threshold multiplied by   tightenCrack_e1_double
0055           e6e2 threshold divided by tightenCrack_e6e2_double
0056
0057
0058    Out of time hits above e4e1_IgnoreOutOfTimeThresh_  are
0059    ignored in topological quantities
0060  */
0061
0062 EcalRecHit::Flags EcalCleaningAlgo::checkTopology(const DetId& id, const EcalRecHitCollection& rhs) {
0063   float a = 0, b = 0, e4e1thresh = 0, ethresh = 0;
0064
0065   if (id.subdetId() == EcalBarrel) {
0066     a = e4e1_a_barrel_;
0067     b = e4e1_b_barrel_;
0068     ethresh = cThreshold_barrel_;
0069
0070   } else if (id.subdetId() == EcalEndcap) {
0071     a = e4e1_a_endcap_;
0072     b = e4e1_b_endcap_;
0073     ethresh = cThreshold_endcap_;
0074   }
0075
0076   // for energies below threshold, we don't apply e4e1 cut
0077   float energy = recHitE(id, rhs, false);
0078
0079   if (energy < ethresh)
0080     return EcalRecHit::kGood;
0081   if (isNearCrack(id) && energy < ethresh * tightenCrack_e1_single_)
0082     return EcalRecHit::kGood;
0083
0084   float e4e1value = e4e1(id, rhs);
0085   e4e1thresh = a * log10(energy) + b;
0086
0087   // near cracks the cut is tighter by a factor
0088   if (isNearCrack(id)) {
0089     e4e1thresh /= tightenCrack_e4e1_single_;
0090   }
0091
0092   // identify spike
0093   if (e4e1value < e4e1thresh)
0094     return EcalRecHit::kWeird;
0095
0096   // now for double spikes
0097
0098   // no checking for double spikes in EE
0099   if (id.subdetId() == EcalEndcap)
0100     return EcalRecHit::kGood;
0101
0102   float e6e2value = e6e2(id, rhs);
0103   float e6e2thresh = e6e2thresh_;
0104   if (isNearCrack(id) && energy < cThreshold_double_ * tightenCrack_e1_double_)
0105     return EcalRecHit::kGood;
0106
0107   if (energy < cThreshold_double_)
0108     return EcalRecHit::kGood;
0109
0110   // near cracks the cut is tighter by a factor
0111   if (id.subdetId() == EcalBarrel && isNearCrack(id))
0112     e6e2thresh /= tightenCrack_e6e2_double_;
0113
0114   // identify double spike
0115   if (e6e2value < e6e2thresh)
0116     return EcalRecHit::kDiWeird;
0117
0118   return EcalRecHit::kGood;
0119 }
0120
0121 float EcalCleaningAlgo::e4e1(const DetId& id, const EcalRecHitCollection& rhs) {
0122   float s4 = 0;
0123   float e1 = recHitE(id, rhs, false);
0124
0125   if (e1 == 0)
0126     return 0;
0127   const std::vector<DetId>& neighs = neighbours(id);
0128   for (size_t i = 0; i < neighs.size(); ++i)
0129     // avoid hits out of time when making s4
0130     s4 += recHitE(neighs[i], rhs, true);
0131
0132   return s4 / e1;
0133 }
0134
0135 float EcalCleaningAlgo::e6e2(const DetId& id, const EcalRecHitCollection& rhs) {
0136   float s4_1 = 0;
0137   float s4_2 = 0;
0138   float e1 = recHitE(id, rhs, false);
0139
0140   float maxene = 0;
0141   DetId maxid;
0142
0143   if (e1 == 0)
0144     return 0;
0145
0146   const std::vector<DetId>& neighs = neighbours(id);
0147
0148   // find the most energetic neighbour ignoring time info
0149   for (size_t i = 0; i < neighs.size(); ++i) {
0150     float ene = recHitE(neighs[i], rhs, false);
0151     if (ene > maxene) {
0152       maxene = ene;
0153       maxid = neighs[i];
0154     }
0155   }
0156
0157   float e2 = maxene;
0158
0159   s4_1 = e4e1(id, rhs) * e1;
0160   s4_2 = e4e1(maxid, rhs) * e2;
0161
0162   return (s4_1 + s4_2) / (e1 + e2) - 1.;
0163 }
0164
0165 float EcalCleaningAlgo::recHitE(const DetId id, const EcalRecHitCollection& recHits, bool useTimingInfo) {
0166   if (id.rawId() == 0)
0167     return 0;
0168
0169   float threshold = e4e1Treshold_barrel_;
0170   if (id.subdetId() == EcalEndcap)
0171     threshold = e4e1Treshold_endcap_;
0172
0173   EcalRecHitCollection::const_iterator it = recHits.find(id);
0174   if (it != recHits.end()) {
0175     float ene = (*it).energy();
0176
0177     // ignore out of time in EB when making e4e1 if so configured
0178     if (useTimingInfo) {
0179       if (id.subdetId() == EcalBarrel && it->checkFlag(EcalRecHit::kOutOfTime) && ene > ignoreOutOfTimeThresh_)
0180         return 0;
0181     }
0182
0183     // ignore hits below threshold
0184     if (ene < threshold)
0185       return 0;
0186
0187     // else return the energy of this hit
0188     return ene;
0189   }
0190   return 0;
0191 }
0192
0193 /// four neighbours in the swiss cross around id
0194 const std::vector<DetId> EcalCleaningAlgo::neighbours(const DetId& id) {
0195   std::vector<DetId> ret;
0196
0197   if (id.subdetId() == EcalBarrel) {
0198     ret.push_back(EBDetId::offsetBy(id, 1, 0));
0199     ret.push_back(EBDetId::offsetBy(id, -1, 0));
0200     ret.push_back(EBDetId::offsetBy(id, 0, 1));
0201     ret.push_back(EBDetId::offsetBy(id, 0, -1));
0202   }
0203   // nobody understands what polymorphism is for, sgrunt !
0204   else if (id.subdetId() == EcalEndcap) {
0205     ret.push_back(EEDetId::offsetBy(id, 1, 0));
0206     ret.push_back(EEDetId::offsetBy(id, -1, 0));
0207     ret.push_back(EEDetId::offsetBy(id, 0, 1));
0208     ret.push_back(EEDetId::offsetBy(id, 0, -1));
0209   }
0210
0211   return ret;
0212 }
0213
0214 bool EcalCleaningAlgo::isNearCrack(const DetId& id) {
0215   if (id.subdetId() == EcalEndcap) {
0216     return EEDetId::isNextToRingBoundary(id);
0217   } else {
0218     return EBDetId::isNextToBoundary(id);
0219   }
0220 }
0221
0222 void EcalCleaningAlgo::setFlags(EcalRecHitCollection& rhs) {
0223   EcalRecHitCollection::iterator rh;
0224   //changing the collection on place
0225   for (rh = rhs.begin(); rh != rhs.end(); ++rh) {
0226     EcalRecHit::Flags state = checkTopology(rh->id(), rhs);
0227     if (state != EcalRecHit::kGood) {
0228       rh->unsetFlag(EcalRecHit::kGood);
0229       rh->setFlag(state);
0230     }
0231   }
0232 }