EcalRecProducers/plugins/EcalUncalibRecHitWorkerFixedAlphaBetaFit.cc

0001 /** \class EcalAnalFitUncalibRecHitProducer
0002  *   produce ECAL uncalibrated rechits from dataframes with the analytic specific fit method
0003  *   with alfa and beta fixed.
0004  *
0005  *  \author A. Ghezzi, Mar 2006
0006  *
0007  */
0008
0009 #include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
0010 #include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"
0011 #include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
0012 #include "CondFormats/EcalObjects/interface/EcalMGPAGainRatio.h"
0013 #include "CondFormats/EcalObjects/interface/EcalPedestals.h"
0014 #include "DataFormats/Common/interface/Handle.h"
0015 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0016 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0017 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0018 #include "FWCore/Framework/interface/ESHandle.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0024 #include "FWCore/Utilities/interface/ESGetToken.h"
0025 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitFixedAlphaBetaAlgo.h"
0026 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitWorkerRunOneDigiBase.h"
0027
0028 #include <cmath>
0029 #include <fstream>
0030 #include <iostream>
0031
0032 class EcalUncalibRecHitWorkerFixedAlphaBetaFit : public EcalUncalibRecHitWorkerRunOneDigiBase {
0033 public:
0034   EcalUncalibRecHitWorkerFixedAlphaBetaFit(const edm::ParameterSet& ps, edm::ConsumesCollector&);
0035   EcalUncalibRecHitWorkerFixedAlphaBetaFit(){};
0036   ~EcalUncalibRecHitWorkerFixedAlphaBetaFit() override{};
0037
0038   void set(const edm::EventSetup& es) override;
0039   bool run(const edm::Event& evt,
0040            const EcalDigiCollection::const_iterator& digi,
0041            EcalUncalibratedRecHitCollection& result) override;
0042
0043   edm::ParameterSetDescription getAlgoDescription() override;
0044
0045 private:
0046   double AmplThrEB_;
0047   double AmplThrEE_;
0048
0049   EcalUncalibRecHitFixedAlphaBetaAlgo<EBDataFrame> algoEB_;
0050   EcalUncalibRecHitFixedAlphaBetaAlgo<EEDataFrame> algoEE_;
0051
0052   double alphaEB_;
0053   double betaEB_;
0054   double alphaEE_;
0055   double betaEE_;
0056   std::vector<std::vector<std::pair<double, double> > >
0057       alphaBetaValues_;  // List of alpha and Beta values [SM#][CRY#](alpha, beta)
0058   bool useAlphaBetaArray_;
0059   std::string alphabetaFilename_;
0060
0061   bool setAlphaBeta();  // Sets the alphaBetaValues_ vectors by the values provided in alphabetaFilename_
0062
0063   edm::ESHandle<EcalGainRatios> pRatio;
0064   edm::ESHandle<EcalPedestals> pedHandle;
0065   edm::ESGetToken<EcalGainRatios, EcalGainRatiosRcd> ratiosToken_;
0066   edm::ESGetToken<EcalPedestals, EcalPedestalsRcd> pedestalsToken_;
0067 };
0068
0069 EcalUncalibRecHitWorkerFixedAlphaBetaFit::EcalUncalibRecHitWorkerFixedAlphaBetaFit(const edm::ParameterSet& ps,
0070                                                                                    edm::ConsumesCollector& c)
0071     : EcalUncalibRecHitWorkerRunOneDigiBase(ps, c),
0072       ratiosToken_(c.esConsumes<EcalGainRatios, EcalGainRatiosRcd>()),
0073       pedestalsToken_(c.esConsumes<EcalPedestals, EcalPedestalsRcd>()) {
0074   alphaEB_ = ps.getParameter<double>("alphaEB");
0075   betaEB_ = ps.getParameter<double>("betaEB");
0076   alphaEE_ = ps.getParameter<double>("alphaEE");
0077   betaEE_ = ps.getParameter<double>("betaEE");
0078
0079   alphabetaFilename_ = ps.getUntrackedParameter<std::string>("AlphaBetaFilename");
0080   useAlphaBetaArray_ = setAlphaBeta();  // set crystalwise values of alpha and beta
0081   if (!useAlphaBetaArray_) {
0082     edm::LogInfo("EcalUncalibRecHitError") << " No alfa-beta file found. Using the deafult values.";
0083   }
0084
0085   algoEB_.SetMinAmpl(ps.getParameter<double>("MinAmplBarrel"));
0086   algoEE_.SetMinAmpl(ps.getParameter<double>("MinAmplEndcap"));
0087
0088   bool dyn_pede = ps.getParameter<bool>("UseDynamicPedestal");
0089   algoEB_.SetDynamicPedestal(dyn_pede);
0090   algoEE_.SetDynamicPedestal(dyn_pede);
0091 }
0092
0093 void EcalUncalibRecHitWorkerFixedAlphaBetaFit::set(const edm::EventSetup& es) {
0094   // Gain Ratios
0095   LogDebug("EcalUncalibRecHitDebug") << "fetching gainRatios....";
0096   pRatio = es.getHandle(ratiosToken_);
0097   LogDebug("EcalUncalibRecHitDebug") << "done.";
0098
0099   // fetch the pedestals from the cond DB via EventSetup
0100   LogDebug("EcalUncalibRecHitDebug") << "fetching pedestals....";
0101   pedHandle = es.getHandle(pedestalsToken_);
0102   LogDebug("EcalUncalibRecHitDebug") << "done.";
0103 }
0104
0105 //Sets the alphaBetaValues_ vectors by the values provided in alphabetaFilename_
0106 bool EcalUncalibRecHitWorkerFixedAlphaBetaFit::setAlphaBeta() {
0107   std::ifstream file(alphabetaFilename_.c_str());
0108   if (!file.is_open())
0109     return false;
0110
0111   alphaBetaValues_.resize(36);
0112
0113   char buffer[100];
0114   int sm, cry, ret;
0115   float a, b;
0116   std::pair<double, double> p(-1, -1);
0117
0118   while (!file.getline(buffer, 100).eof()) {
0119     ret = sscanf(buffer, "%d %d %f %f", &sm, &cry, &a, &b);
0120     if ((ret != 4) || (sm <= 0) || (sm > 36) || (cry <= 0) || (cry > 1700)) {
0121       // send warning
0122       continue;
0123     }
0124
0125     if (alphaBetaValues_[sm - 1].empty()) {
0126       alphaBetaValues_[sm - 1].resize(1700, p);
0127     }
0128     alphaBetaValues_[sm - 1][cry - 1].first = a;
0129     alphaBetaValues_[sm - 1][cry - 1].second = b;
0130   }
0131
0132   file.close();
0133   return true;
0134 }
0135
0136 bool EcalUncalibRecHitWorkerFixedAlphaBetaFit::run(const edm::Event& evt,
0137                                                    const EcalDigiCollection::const_iterator& itdg,
0138                                                    EcalUncalibratedRecHitCollection& result) {
0139   const EcalGainRatioMap& gainMap = pRatio.product()->getMap();  // map of gain ratios
0140   EcalGainRatioMap::const_iterator gainIter;                     // gain iterator
0141   EcalMGPAGainRatio aGain;                                       // gain object for a single xtal
0142
0143   const EcalPedestalsMap& pedMap = pedHandle.product()->getMap();  // map of pedestals
0144   EcalPedestalsMapIterator pedIter;                                // pedestal iterator
0145   EcalPedestals::Item aped;                                        // pedestal object for a single xtal
0146
0147   DetId detid(itdg->id());
0148
0149   // find pedestals for this channel
0150   //LogDebug("EcalUncalibRecHitDebug") << "looking up pedestal for crystal: " << itdg->id();
0151   pedIter = pedMap.find(itdg->id());
0152   if (pedIter != pedMap.end()) {
0153     aped = (*pedIter);
0154   } else {
0155     edm::LogError("EcalUncalibRecHitWorkerFixedAlphaBetaFit") << "error!! could not find pedestals for channel: ";
0156     if (detid.subdetId() == EcalBarrel) {
0157       edm::LogError("EcalUncalibRecHitWorkerFixedAlphaBetaFit") << EBDetId(detid);
0158     } else {
0159       edm::LogError("EcalUncalibRecHitWorkerFixedAlphaBetaFit") << EEDetId(detid);
0160     }
0161     edm::LogError("EcalUncalibRecHitWorkerFixedAlphaBetaFit") << "\n  no uncalib rechit will be made for this digi!";
0162     return false;
0163   }
0164   double pedVec[3];
0165   pedVec[0] = aped.mean_x12;
0166   pedVec[1] = aped.mean_x6;
0167   pedVec[2] = aped.mean_x1;
0168
0169   // find gain ratios
0170   //LogDebug("EcalUncalibRecHitDebug") << "looking up gainRatios for crystal: " << EBDetId(itdg->id()) ; // FIXME!!!!!!!!
0171   gainIter = gainMap.find(itdg->id());
0172   if (gainIter != gainMap.end()) {
0173     aGain = (*gainIter);
0174   } else {
0175     edm::LogError("EcalUncalibRecHitWorkerFixedAlphaBetaFit") << "error!! could not find gain ratios for channel: ";
0176     if (detid.subdetId() == EcalBarrel) {
0177       edm::LogError("EcalUncalibRecHitWorkerFixedAlphaBetaFit") << EBDetId(detid);
0178     } else {
0179       edm::LogError("EcalUncalibRecHitWorkerFixedAlphaBetaFit") << EEDetId(detid);
0180     }
0181     edm::LogError("EcalUncalibRecHitWorkerFixedAlphaBetaFit") << "\n  no uncalib rechit will be made for this digi!";
0182     return false;
0183   }
0184   double gainRatios[3];
0185   gainRatios[0] = 1.;
0186   gainRatios[1] = aGain.gain12Over6();
0187   gainRatios[2] = aGain.gain6Over1() * aGain.gain12Over6();
0188
0189   if (detid.subdetId() == EcalBarrel) {
0190     // Define Alpha and Beta either by stored values or by default universal values
0191     EBDetId ebDetId(detid);
0192     double a, b;
0193     if (useAlphaBetaArray_) {
0194       if (!alphaBetaValues_[ebDetId.ism() - 1].empty()) {
0195         a = alphaBetaValues_[ebDetId.ism() - 1][ebDetId.ic() - 1].first;
0196         b = alphaBetaValues_[ebDetId.ism() - 1][ebDetId.ic() - 1].second;
0197         if ((a == -1) && (b == -1)) {
0198           a = alphaEB_;
0199           b = betaEB_;
0200         }
0201       } else {
0202         a = alphaEB_;
0203         b = betaEB_;
0204       }
0205     } else {
0206       a = alphaEB_;
0207       b = betaEB_;
0208     }
0209     algoEB_.SetAlphaBeta(a, b);
0210     result.push_back(algoEB_.makeRecHit(*itdg, pedVec, gainRatios, nullptr, nullptr));
0211   } else {
0212     //FIX ME load in a and b from a file
0213     algoEE_.SetAlphaBeta(alphaEE_, betaEE_);
0214     result.push_back(algoEE_.makeRecHit(*itdg, pedVec, gainRatios, nullptr, nullptr));
0215   }
0216   return true;
0217 }
0218
0219 edm::ParameterSetDescription EcalUncalibRecHitWorkerFixedAlphaBetaFit::getAlgoDescription() {
0220   edm::ParameterSetDescription psd;
0221
0222   psd.addNode(edm::ParameterDescription<double>("alphaEB", 1.138, true) and
0223               edm::ParameterDescription<double>("alphaEE", 1.89, true) and
0224               edm::ParameterDescription<std::string>("AlphaBetaFilename", "NOFILE", false) and
0225               edm::ParameterDescription<double>("betaEB", 1.655, true) and
0226               edm::ParameterDescription<double>("MinAmplEndcap", 14.0, true) and
0227               edm::ParameterDescription<double>("MinAmplBarrel", 8.0, true) and
0228               edm::ParameterDescription<double>("betaEE", 1.4, true) and
0229               edm::ParameterDescription<bool>("UseDynamicPedestal", true, true));
0230
0231   return psd;
0232 }
0233
0234 #include "FWCore/Framework/interface/MakerMacros.h"
0235 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitWorkerFactory.h"
0236 DEFINE_EDM_PLUGIN(EcalUncalibRecHitWorkerFactory,
0237                   EcalUncalibRecHitWorkerFixedAlphaBetaFit,
0238                   "EcalUncalibRecHitWorkerFixedAlphaBetaFit");
0239 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitFillDescriptionWorkerFactory.h"
0240 DEFINE_EDM_PLUGIN(EcalUncalibRecHitFillDescriptionWorkerFactory,
0241                   EcalUncalibRecHitWorkerFixedAlphaBetaFit,
0242                   "EcalUncalibRecHitWorkerFixedAlphaBetaFit");