HcalRecAlgos/test/HcalRecHitReflagger.cc

0001 // -*- C++ -*-

0002 //

0003 // Package:    HcalRecHitReflagger

0004 // Class:      HcalRecHitReflagger

0005 //

0006 /**\class HcalRecHitReflagger HcalRecHitReflagger.cc

0007

0008  Description: [Creates a new collection of rechits from an existing collection, and changes rechit flags according to new set of user-supplied conditions.  ]

0009

0010 */
0011 //

0012 // Original Author:  Dinko Ferencek,8 R-004,+41227676479,  Jeff Temple, 6-1-027

0013 //         Created:  Thu Mar 11 13:42:11 CET 2010

0014 //

0015 //

0016
0017 // system include files

0018 #include <iostream>
0019 #include <memory>
0020
0021 // user include files

0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/one/EDProducer.h"
0024
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030
0031 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0032 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0033 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0034 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0035
0036 #include "CondFormats/HcalObjects/interface/HcalChannelStatus.h"
0037 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0038 #include "CondFormats/HcalObjects/interface/HcalCondObjectContainer.h"
0039 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
0040
0041 using namespace std;
0042 using namespace edm;
0043
0044 //

0045 // class declaration

0046 //

0047
0048 class HcalRecHitReflagger : public edm::one::EDProducer<> {
0049 public:
0050   explicit HcalRecHitReflagger(const edm::ParameterSet&);
0051   ~HcalRecHitReflagger();
0052
0053 private:
0054   void beginJob() override;
0055   void produce(edm::Event&, const edm::EventSetup&) override;
0056   void endJob() override;
0057
0058   // Threshold function gets values from polynomial-parameterized functions

0059   double GetThreshold(const int base, const std::vector<double>& params);
0060   double GetThreshold(const double base, const std::vector<double>& params);
0061
0062   // Perform a check of a rechit's S9/S1 value compared to a given threshold

0063   bool CheckS9S1(const HFRecHit& hf);
0064
0065   // Perform a check of a rechit's |L-S|/|L+S| ratio, compared to a paremeterized threshold

0066   bool CheckPET(const HFRecHit& hf);
0067
0068   // Get S9S1, PET values

0069   double GetS9S1value(const HFRecHit& hf);
0070   double GetPETvalue(const HFRecHit& hf);
0071   double GetSlope(const int ieta, const std::vector<double>& params);
0072
0073   // ----------member data ---------------------------

0074   const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tokTopo_;
0075   const edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> tokChan_;
0076   const HcalTopology* topo;
0077   edm::InputTag hfInputLabel_;
0078   edm::EDGetTokenT<HFRecHitCollection> tok_hf_;
0079   int hfFlagBit_;
0080
0081   // Select the test you wish to run

0082   bool hfBitAlwaysOn_;
0083   bool hfBitAlwaysOff_;
0084   bool hf_Algo3test_;
0085   bool hf_Algo2test_;
0086   bool hf_Algo3_AND_Algo2_;
0087   bool hf_Algo3_OR_Algo2_;
0088
0089   std::vector<double> S9S1_EnergyThreshLong_;
0090   std::vector<double> S9S1_ETThreshLong_;
0091   std::vector<double> S9S1_EnergyThreshShort_;
0092   std::vector<double> S9S1_ETThreshShort_;
0093   std::vector<double> S9S1_optimumslope_;
0094   std::vector<double> PET_EnergyThreshLong_;
0095   std::vector<double> PET_ETThreshLong_;
0096   std::vector<double> PET_EnergyThreshShort_;
0097   std::vector<double> PET_ETThreshShort_;
0098   std::vector<double> PET_ratiocut_;
0099   int debug_;
0100
0101   // Store channels marked as bad in the channel status map

0102   std::map<HcalDetId, unsigned int> badstatusmap;
0103
0104   edm::Handle<HFRecHitCollection> hfRecHits;
0105 };
0106
0107 //

0108 // constants, enums and typedefs

0109 //

0110
0111 //

0112 // static data member definitions

0113 //

0114
0115 //

0116 // constructors and destructor

0117 //

0118 HcalRecHitReflagger::HcalRecHitReflagger(const edm::ParameterSet& ps)
0119     : tokTopo_(esConsumes<HcalTopology, HcalRecNumberingRecord>()),
0120       tokChan_(esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"))) {
0121   //register your products

0122   produces<HFRecHitCollection>();
0123
0124   hfInputLabel_ = ps.getUntrackedParameter<InputTag>("hfInputLabel", edm::InputTag("hfreco"));
0125   tok_hf_ = consumes<HFRecHitCollection>(hfInputLabel_);
0126   hfFlagBit_ = ps.getUntrackedParameter<int>("hfFlagBit", HcalCaloFlagLabels::UserDefinedBit0);
0127   debug_ = ps.getUntrackedParameter<int>("debug", 0);
0128
0129   // sanity check bits -- turn bit on or off always

0130   hfBitAlwaysOn_ = ps.getUntrackedParameter<bool>("hfBitAlwaysOn", false);
0131   hfBitAlwaysOff_ = ps.getUntrackedParameter<bool>("hfBitAlwaysOff", false);
0132
0133   // Enable/disable various tests

0134   hf_Algo3test_ = ps.getUntrackedParameter<bool>("hf_Algo3test", false);
0135   hf_Algo2test_ = ps.getUntrackedParameter<bool>("hf_Algo2test", false);
0136   hf_Algo3_AND_Algo2_ = ps.getUntrackedParameter<bool>("hf_Algo3_AND_Algo2", false);
0137   hf_Algo3_OR_Algo2_ = ps.getUntrackedParameter<bool>("hf_Algo3_OR_Algo2", false);
0138
0139   // Get thresholds for each test

0140   const edm::ParameterSet& S9S1 = ps.getParameter<edm::ParameterSet>("hf_S9S1_params");
0141   S9S1_EnergyThreshLong_ = S9S1.getUntrackedParameter<std::vector<double> >("S9S1_EnergyThreshLong");
0142   S9S1_ETThreshLong_ = S9S1.getUntrackedParameter<std::vector<double> >("S9S1_ETThreshLong");
0143   S9S1_EnergyThreshShort_ = S9S1.getUntrackedParameter<std::vector<double> >("S9S1_EnergyThreshShort");
0144   S9S1_ETThreshShort_ = S9S1.getUntrackedParameter<std::vector<double> >("S9S1_ETThreshShort");
0145   S9S1_optimumslope_ = S9S1.getParameter<std::vector<double> >("S9S1_optimumslope");
0146
0147   const edm::ParameterSet& PET = ps.getParameter<edm::ParameterSet>("hf_PET_params");
0148   PET_EnergyThreshLong_ = PET.getUntrackedParameter<std::vector<double> >("PET_EnergyThreshLong");
0149   PET_ETThreshLong_ = PET.getUntrackedParameter<std::vector<double> >("PET_ETThreshLong");
0150   PET_EnergyThreshShort_ = PET.getUntrackedParameter<std::vector<double> >("PET_EnergyThreshShort");
0151   PET_ETThreshShort_ = PET.getUntrackedParameter<std::vector<double> >("PET_ETThreshShort");
0152   PET_ratiocut_ = PET.getParameter<std::vector<double> >("PET_ratiocut");
0153 }
0154
0155 HcalRecHitReflagger::~HcalRecHitReflagger() {
0156   // do anything here that needs to be done at desctruction time

0157   // (e.g. close files, deallocate resources etc.)

0158 }
0159
0160 //

0161 // member functions

0162 //

0163 // ------------ method called to produce the data  ------------

0164 void HcalRecHitReflagger::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0165   // read HF RecHits

0166   if (!iEvent.getByToken(tok_hf_, hfRecHits)) {
0167     if (debug_ > 0)
0168       std::cout << "Unable to find HFRecHitCollection with label '" << hfInputLabel_ << "'" << std::endl;
0169     return;
0170   }
0171
0172   //get the hcal topology

0173   topo = &iSetup.getData(tokTopo_);
0174
0175   //get channel quality conditions

0176   const HcalChannelQuality& chanquality_ = iSetup.getData(tokChan_);
0177
0178   std::vector<DetId> mydetids = chanquality_.getAllChannels();
0179   for (std::vector<DetId>::const_iterator i = mydetids.begin(); i != mydetids.end(); ++i) {
0180     if (i->det() != DetId::Hcal)
0181       continue;  // not an hcal cell

0182     HcalDetId id = HcalDetId(*i);
0183     int status = (chanquality_.getValues(id))->getValue();
0184     if ((status & (1 << HcalChannelStatus::HcalCellDead)) == 0)
0185       continue;
0186     badstatusmap[id] = status;
0187   }
0188
0189   // prepare the output HF RecHit collection

0190   auto pOut = std::make_unique<HFRecHitCollection>();
0191
0192   // loop over rechits, and set the new bit you wish to use

0193   for (HFRecHitCollection::const_iterator recHit = hfRecHits->begin(); recHit != hfRecHits->end(); ++recHit) {
0194     HFRecHit newhit = (HFRecHit)(*recHit);
0195     // Set bit to be on for all hits

0196     if (hfBitAlwaysOn_)
0197       newhit.setFlagField(1, hfFlagBit_);
0198
0199     // Set bit to be off for all hits

0200     else if (hfBitAlwaysOff_)
0201       newhit.setFlagField(0, hfFlagBit_);
0202
0203     else {
0204       // Check booleans -- true booleans mean that the channel seems noisy

0205       bool Algo2bool = false;  // Algorithm 2 is PET testing for both long & short

0206       bool Algo3bool = false;  // Algorithm 3 is PET testing for short, S9S1 for long

0207
0208       HcalDetId id(newhit.detid().rawId());
0209       int depth = newhit.id().depth();
0210       int ieta = newhit.id().ieta();
0211
0212       // Code below performs Algo2 or Algo3 tests;  modify as you see fit

0213
0214       if (depth == 2)  // short fibers -- use PET test everywhere

0215       {
0216         Algo2bool = CheckPET(newhit);
0217         Algo3bool = Algo2bool;
0218       } else if (depth == 1)  // long fibers

0219       {
0220         Algo2bool = CheckPET(newhit);
0221         abs(ieta) == 29 ? Algo3bool = Algo2bool : Algo3bool = CheckS9S1(newhit);
0222       }
0223
0224       int flagvalue = -1;  // new value for flag  (-1 means flag not changed)

0225
0226       // run through test combinations

0227       if (hf_Algo3_AND_Algo2_ == true)
0228         (Algo2bool && Algo3bool) == true ? flagvalue = 1 : flagvalue = 0;
0229       else if (hf_Algo3_OR_Algo2_ == true)
0230         (Algo2bool || Algo3bool) == true ? flagvalue = 1 : flagvalue = 0;
0231       else if (hf_Algo3test_ == true)
0232         Algo3bool == true ? flagvalue = 1 : flagvalue = 0;
0233       else if (hf_Algo2test_ == true)
0234         Algo2bool == true ? flagvalue = 1 : flagvalue = 0;
0235
0236       // Set flag bit based on test; if no tests set, don't change flag

0237       if (flagvalue != -1)
0238         newhit.setFlagField(flagvalue, hfFlagBit_);
0239     }  // hfBitAlwaysOn/Off bits not set; run other tests

0240
0241     // Add rechit to collection

0242     pOut->push_back(newhit);
0243   }
0244
0245   // put the re-flagged HF RecHit collection into the Event

0246   iEvent.put(std::move(pOut));
0247 }
0248
0249 // ------------ method called once each job just before starting event loop  ------------

0250 void HcalRecHitReflagger::beginJob() {}
0251
0252 // ------------ method called once each job just after ending the event loop  ------------

0253 void HcalRecHitReflagger::endJob() {}
0254
0255 bool HcalRecHitReflagger::CheckPET(const HFRecHit& hf) {
0256   /*

0257     Runs 'PET' test:  Computes ratio |L-S|/(L+S) for channel passing ieta-dependent energy and ET cuts.

0258     Channel is marked as noise if ratio exceeds ratio R(E)

0259   */
0260
0261   HcalDetId id(hf.detid().rawId());
0262   int iphi = id.iphi();
0263   int ieta = id.ieta();
0264   double energy = hf.energy();
0265   int depth = id.depth();
0266   std::pair<double, double> etas = topo->etaRange(HcalForward, abs(ieta));
0267   double fEta = 0.5 * (etas.first + etas.second);  // calculate eta as average of eta values at ieta boundaries

0268   double ET = energy / cosh(fEta);
0269   double threshold = 0;
0270
0271   if (debug_ > 0)
0272     std::cout << "<RechitReflagger::CheckPET>  energy = " << energy << "  ieta = " << ieta << std::endl;
0273
0274   if (depth == 1)  //long fibers

0275   {
0276     // Check that energy, ET above threshold

0277     threshold = GetThreshold(ieta, PET_EnergyThreshLong_);
0278     if (energy < threshold)
0279       return false;
0280     threshold = GetThreshold(ieta, PET_ETThreshLong_);
0281     if (ET < threshold)
0282       return false;
0283
0284     // Get threshold cut to use at this energy, and compute the |L-S|/(L+S) ratio

0285
0286     double PETcut = GetThreshold(energy, PET_ratiocut_);
0287     double PETvalue = GetPETvalue(hf);
0288     if (debug_ > 1)
0289       std::cout << "  <HcalRecHitReflagger::CheckPET>  ieta = " << ieta
0290                 << "  energy threshold = " << GetThreshold(ieta, PET_EnergyThreshLong_)
0291                 << "  ET threshold = " << GetThreshold(ieta, PET_ETThreshLong_) << "   ENERGY = " << energy
0292                 << "   PET threshold = " << PETcut << std::endl;
0293
0294     if (PETvalue > PETcut) {
0295       if (debug_ > 2)
0296         std::cout << "***** FOUND CHANNEL FAILING PET CUT!  CHANNEL = " << ieta << ",  " << iphi << ",  " << depth
0297                   << std::endl;
0298       return true;  // channel looks noisy

0299     } else
0300       return false;  // channel doesn't look noisy

0301   }  // if (depth==1)

0302
0303   else if (depth == 2)  //short fibers

0304   {
0305     threshold = GetThreshold(ieta, PET_EnergyThreshShort_);
0306     if (energy < threshold)
0307       return false;
0308     threshold = GetThreshold(ieta, PET_ETThreshShort_);
0309     if (ET < threshold)
0310       return false;
0311
0312     double PETcut = GetThreshold(energy, PET_ratiocut_);
0313     double PETvalue = GetPETvalue(hf);
0314     if (debug_ > 1)
0315       std::cout << "  <HcalRecHitReflagger::CheckPET>  ieta = " << ieta
0316                 << "  energy threshold = " << GetThreshold(ieta, PET_EnergyThreshShort_)
0317                 << "  ET threshold = " << GetThreshold(ieta, PET_ETThreshShort_) << "   ENERGY = " << energy
0318                 << "   PET threshold = " << PETcut << endl;
0319
0320     if (PETvalue > PETcut) {
0321       if (debug_ > 2)
0322         std::cout << "***** FOUND CHANNEL FAILING PET CUT!  CHANNEL = " << ieta << ",  " << iphi << ",  " << depth
0323                   << std::endl;
0324
0325       return true;  // looks noisy

0326     } else
0327       return false;  //doesn't look noisy

0328   }  // else if (depth==2)

0329   return false;  // should never reach this part

0330 }  // CheckPET(HFRecHit& hf)

0331
0332 bool HcalRecHitReflagger::CheckS9S1(const HFRecHit& hf) {
0333   // Apply S9/S1 test  -- only sensible for long fibers at the moment!

0334
0335   HcalDetId id(hf.detid().rawId());
0336   int iphi = id.iphi();
0337   int ieta = id.ieta();
0338   double energy = hf.energy();
0339   int depth = id.depth();
0340   std::pair<double, double> etas = topo->etaRange(HcalForward, ieta);
0341   double fEta = 0.5 * (etas.first + etas.second);  // calculate eta as average of eta values at ieta boundaries

0342   double ET = energy / cosh(fEta);
0343   double threshold = 0;
0344
0345   //double S9S1slope=GetSlope(ieta,S9S1_optimumslope_ ); // get optimized slope for this ieta

0346   double S9S1slope = S9S1_optimumslope_[abs(ieta) - 29];
0347   double S9S1value = GetS9S1value(hf);  // get the ratio S9/S1

0348
0349   if (debug_ > 0)
0350     std::cout << "<RechitReflagger::CheckS9S1>  energy = " << energy << "  ieta = " << ieta
0351               << "  S9S1slope = " << S9S1slope << "  S9S1value = " << S9S1value << std::endl;
0352
0353   if (depth == 1)  // long fiber

0354   {
0355     // Step 1:  Check that energy above threshold

0356     threshold = GetThreshold(ieta, S9S1_EnergyThreshLong_);
0357     if (energy < threshold)
0358       return false;
0359     threshold = GetThreshold(ieta, S9S1_ETThreshLong_);
0360     if (ET < threshold)
0361       return false;
0362
0363     // S9S1 cut has form [0]+[1]*log[E],

0364     double S9S1cut = -1. * S9S1slope * log(GetThreshold(ieta, PET_EnergyThreshLong_)) + S9S1slope * (log(energy));
0365
0366     if (S9S1value < S9S1cut) {
0367       if (debug_ > 0)
0368         std::cout << "***** FOUND CHANNEL FAILING S9S1 CUT!  CHANNEL = " << ieta << ",  " << iphi << ",  " << depth
0369                   << std::endl;
0370       return true;
0371     } else
0372       return false;
0373   }
0374
0375   // ****************WARNING!!!!!****************  SHORT FIBERS DON'T HAVE TUNED THRESHOLDS!  USE PET CUT FOR SHORT FIBERS! ******* //

0376   else if (depth == 2)  //short fiber

0377   {
0378     // Step 1:  Check that energy above threshold

0379     threshold = GetThreshold(ieta, S9S1_EnergyThreshShort_);
0380     if (energy < threshold)
0381       return false;
0382     threshold = GetThreshold(ieta, S9S1_ETThreshShort_);
0383     if (ET < threshold)
0384       return false;
0385
0386     // S9S1 cut has form [0]+[1]*log[E],

0387     double S9S1cut = -1. * S9S1slope * log(GetThreshold(ieta, PET_EnergyThreshShort_)) + S9S1slope * (log(energy));
0388
0389     if (S9S1value < S9S1cut) {
0390       if (debug_ > 0)
0391         std::cout << "***** FOUND CHANNEL FAILING S9S1 CUT!  CHANNEL = " << ieta << ",  " << iphi << ",  " << depth
0392                   << std::endl;
0393       return true;
0394     } else
0395       return false;
0396   }
0397   return false;
0398 }
0399
0400 double HcalRecHitReflagger::GetS9S1value(const HFRecHit& hf) {
0401   // sum all energies around cell hf (including hf itself)

0402   // Then subtract hf energy, and divide by this energy to form ratio S9/S1

0403
0404   HcalDetId id(hf.detid().rawId());
0405   int iphi = id.iphi();
0406   int ieta = id.ieta();
0407   double energy = hf.energy();
0408   int depth = id.depth();
0409   if (debug_ > 0)
0410     std::cout << "\t<HcalRecHitReflagger::GetS9S1value>  Channel = (" << iphi << ",  " << ieta << ",  " << depth
0411               << ")  :  Energy = " << energy << std::endl;
0412
0413   double S9S1 = 0;  // starting value

0414
0415   int testphi = -99;
0416
0417   // Part A:  Check fixed iphi, and vary ieta

0418   for (int d = 1; d <= 2; ++d)  // depth loop

0419   {
0420     for (int i = ieta - 1; i <= ieta + 1; ++i)  // ieta loop

0421     {
0422       testphi = iphi;
0423       // Special case when ieta=39, since ieta=40 only has phi values at 3,7,11,...

0424       // phi=3 covers 3,4,5,6

0425       if (abs(ieta) == 39 && abs(i) > 39 && testphi % 4 == 1)
0426         testphi -= 2;
0427       while (testphi < 0)
0428         testphi += 72;
0429       if (i == ieta && d == depth)
0430         continue;  // don't add the cell itself

0431       // Look to see if neighbor is in rechit collection

0432       HcalDetId neighbor(HcalForward, i, testphi, d);
0433       HFRecHitCollection::const_iterator neigh = hfRecHits->find(neighbor);
0434       if (neigh != hfRecHits->end())
0435         S9S1 += neigh->energy();
0436     }
0437   }
0438
0439   // Part B: Fix ieta, and loop over iphi.  A bit more tricky, because of iphi wraparound and different segmentation at 40, 41

0440
0441   int phiseg = 2;  // 10 degree segmentation for most of HF (1 iphi unit = 5 degrees)

0442   if (abs(ieta) > 39)
0443     phiseg = 4;  // 20 degree segmentation for |ieta|>39

0444   for (int d = 1; d <= 2; ++d) {
0445     for (int i = iphi - phiseg; i <= iphi + phiseg; i += phiseg) {
0446       if (i == iphi)
0447         continue;  // don't add the cell itself, or its depthwise partner (which is already counted above)

0448       testphi = i;
0449       // Our own modular function, since default produces results -1%72 = -1

0450       while (testphi < 0)
0451         testphi += 72;
0452       while (testphi > 72)
0453         testphi -= 72;
0454       // Look to see if neighbor is in rechit collection

0455       HcalDetId neighbor(HcalForward, ieta, testphi, d);
0456       HFRecHitCollection::const_iterator neigh = hfRecHits->find(neighbor);
0457       if (neigh != hfRecHits->end())
0458         S9S1 += neigh->energy();
0459     }
0460   }
0461
0462   if (abs(ieta) == 40)  // add extra cells for 39/40 boundary due to increased phi size at ieta=40.

0463   {
0464     for (int d = 1; d <= 2; ++d)  // add cells from both depths!

0465     {
0466       HcalDetId neighbor(HcalForward, 39 * abs(ieta) / ieta, (iphi + 2) % 72, d);
0467       HFRecHitCollection::const_iterator neigh = hfRecHits->find(neighbor);
0468       if (neigh != hfRecHits->end())
0469         S9S1 += neigh->energy();
0470     }
0471   }
0472
0473   S9S1 /= energy;  // divide to form actual ratio

0474   return S9S1;
0475 }  // GetS9S1value

0476
0477 double HcalRecHitReflagger::GetPETvalue(const HFRecHit& hf) {
0478   HcalDetId id(hf.detid().rawId());
0479   int iphi = id.iphi();
0480   int ieta = id.ieta();
0481   int depth = id.depth();
0482   double energy = hf.energy();
0483   if (debug_ > 0)
0484     std::cout << "\t<HcalRecHitReflagger::GetPETvalue>  Channel = (" << iphi << ",  " << ieta << ",  " << depth
0485               << ")  :  Energy = " << energy << std::endl;
0486
0487   HcalDetId pId(HcalForward, ieta, iphi, 3 - depth);  // get partner;

0488   // Check if partner is in known dead cell list; if so, don't flag

0489   if (badstatusmap.find(pId) != badstatusmap.end())
0490     return 0;
0491
0492   double epartner = 0;  // assume no partner

0493   // TO DO:  Add protection against cells marked dead in status database

0494   HFRecHitCollection::const_iterator part = hfRecHits->find(pId);
0495   if (part != hfRecHits->end())
0496     epartner = part->energy();
0497   return 1. * (energy - epartner) / (energy + epartner);
0498 }
0499
0500 double HcalRecHitReflagger::GetThreshold(const int base, const std::vector<double>& params) {
0501   double thresh = 0;
0502   for (unsigned int i = 0; i < params.size(); ++i) {
0503     thresh += params[i] * pow(fabs(base), (int)i);
0504   }
0505   return thresh;
0506 }
0507
0508 double HcalRecHitReflagger::GetThreshold(const double base, const std::vector<double>& params) {
0509   double thresh = 0;
0510   for (unsigned int i = 0; i < params.size(); ++i) {
0511     thresh += params[i] * pow(fabs(base), (int)i);
0512   }
0513   return thresh;
0514 }
0515
0516 double HcalRecHitReflagger::GetSlope(const int ieta, const std::vector<double>& params) {
0517   double slope = 0;
0518   if (abs(ieta) == 40 && params.size() > 0)
0519     slope = params[0];
0520   else if (abs(ieta) == 41 && params.size() > 1)
0521     slope = params[1];
0522   else if (params.size() > 2) {
0523     for (unsigned int i = 2; i < params.size(); ++i)
0524       slope += params[i] * pow(static_cast<double>(abs(ieta)), static_cast<int>(i) - 2);
0525   }
0526   if (debug_ > 1)
0527     std::cout << "<HcalRecHitReflagger::GetSlope>  ieta = " << ieta << "  slope = " << slope << std::endl;
0528   return slope;
0529 }
0530
0531 //define this as a plug-in

0532 DEFINE_FWK_MODULE(HcalRecHitReflagger);