EcalLaserCorrection/src/EcalLaserDbService.cc

0001 #include <iostream>
0002
0003 #include "FWCore/Utilities/interface/typelookup.h"
0004
0005 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
0006
0007 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEEBGeom.h"
0008 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEEEGeom.h"
0009 // #include "CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h"
0010
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 using namespace std;
0014
0015 EcalLaserDbService::EcalLaserDbService()
0016     : mAlphas_(nullptr),
0017       mAPDPNRatiosRef_(nullptr),
0018       mAPDPNRatios_(nullptr),
0019       mLinearCorrections_(nullptr),
0020       maxExtrapolationTime_(0) {}
0021
0022 const EcalLaserAlphas* EcalLaserDbService::getAlphas() const { return mAlphas_; }
0023
0024 const EcalLaserAPDPNRatiosRef* EcalLaserDbService::getAPDPNRatiosRef() const { return mAPDPNRatiosRef_; }
0025
0026 const EcalLaserAPDPNRatios* EcalLaserDbService::getAPDPNRatios() const { return mAPDPNRatios_; }
0027
0028 const EcalLinearCorrections* EcalLaserDbService::getLinearCorrections() const { return mLinearCorrections_; }
0029
0030 float EcalLaserDbService::getLaserCorrection(DetId const& xid, edm::Timestamp const& iTime) const {
0031   float correctionFactor = 1.0;
0032
0033   const EcalLaserAPDPNRatios::EcalLaserAPDPNRatiosMap& laserRatiosMap = mAPDPNRatios_->getLaserMap();
0034   const EcalLaserAPDPNRatios::EcalLaserTimeStampMap& laserTimeMap = mAPDPNRatios_->getTimeMap();
0035   const EcalLaserAPDPNRatiosRefMap& laserRefMap = mAPDPNRatiosRef_->getMap();
0036   const EcalLaserAlphaMap& laserAlphaMap = mAlphas_->getMap();
0037   const EcalLinearCorrections::EcalValueMap& linearValueMap = mLinearCorrections_->getValueMap();
0038   const EcalLinearCorrections::EcalTimeMap& linearTimeMap = mLinearCorrections_->getTimeMap();
0039
0040   EcalLaserAPDPNRatios::EcalLaserAPDPNpair apdpnpair;
0041   EcalLaserAPDPNRatios::EcalLaserTimeStamp timestamp;
0042   EcalLaserAPDPNref apdpnref;
0043   EcalLaserAlpha alpha;
0044   EcalLinearCorrections::Values linValues;
0045   EcalLinearCorrections::Times linTimes;
0046
0047   if (xid.det() == DetId::Ecal) {
0048   } else {
0049     edm::LogError("EcalLaserDbService") << " DetId is NOT in ECAL" << endl;
0050     return correctionFactor;
0051   }
0052
0053   int iLM;
0054   int xind;
0055   bool isBarrel = true;
0056   if (xid.subdetId() == EcalBarrel) {
0057     EBDetId ebid(xid.rawId());
0058     xind = ebid.hashedIndex();
0059     iLM = MEEBGeom::lmr(ebid.ieta(), ebid.iphi());
0060   } else if (xid.subdetId() == EcalEndcap) {
0061     isBarrel = false;
0062     EEDetId eeid(xid.rawId());
0063     xind = eeid.hashedIndex();
0064     // SuperCrystal coordinates
0065     MEEEGeom::SuperCrysCoord iX = (eeid.ix() - 1) / 5 + 1;
0066     MEEEGeom::SuperCrysCoord iY = (eeid.iy() - 1) / 5 + 1;
0067     iLM = MEEEGeom::lmr(iX, iY, eeid.zside());
0068   } else {
0069     edm::LogError("EcalLaserDbService") << " DetId is NOT in ECAL Barrel or Endcap" << endl;
0070     return correctionFactor;
0071   }
0072   //  std::cout << " LM num ====> " << iLM << endl;
0073
0074   // get alpha, apd/pn ref, apd/pn pairs and timestamps for interpolation
0075
0076 #ifdef VERIFY_LASER
0077   EcalLaserAPDPNRatios::EcalLaserAPDPNRatiosMap::const_iterator itratio = laserRatiosMap.find(xid);
0078   if (itratio != laserRatiosMap.end()) {
0079     apdpnpair = (*itratio);
0080   } else {
0081     edm::LogError("EcalLaserDbService") << "error with laserRatiosMap!" << endl;
0082     return correctionFactor;
0083   }
0084
0085   if (iLM - 1 < (int)laserTimeMap.size()) {
0086     timestamp = laserTimeMap[iLM - 1];
0087   } else {
0088     edm::LogError("EcalLaserDbService") << "error with laserTimeMap!" << endl;
0089     return correctionFactor;
0090   }
0091
0092   EcalLinearCorrections::EcalValueMap::const_iterator itlin = linearValueMap.find(xid);
0093   if (itlin != linearValueMap.end()) {
0094     linValues = (*itlin);
0095   } else {
0096     edm::LogError("EcalLaserDbService") << "error with linearValueMap!" << endl;
0097     return correctionFactor;
0098   }
0099
0100   if (iLM - 1 < (int)linearTimeMap.size()) {
0101     linTimes = linearTimeMap[iLM - 1];
0102   } else {
0103     edm::LogError("EcalLaserDbService") << "error with laserTimeMap!" << endl;
0104     return correctionFactor;
0105   }
0106
0107   EcalLaserAPDPNRatiosRefMap::const_iterator itref = laserRefMap.find(xid);
0108   if (itref != laserRefMap.end()) {
0109     apdpnref = (*itref);
0110   } else {
0111     edm::LogError("EcalLaserDbService") << "error with laserRefMap!" << endl;
0112     return correctionFactor;
0113   }
0114
0115   EcalLaserAlphaMap::const_iterator italpha = laserAlphaMap.find(xid);
0116   if (italpha != laserAlphaMap.end()) {
0117     alpha = (*italpha);
0118   } else {
0119     edm::LogError("EcalLaserDbService") << "error with laserAlphaMap!" << endl;
0120     return correctionFactor;
0121   }
0122
0123 #else
0124
0125   // waiting for templated lambdas
0126   auto getCond = [=](EcalFloatCondObjectContainer const& cond) -> float {
0127     return isBarrel ? cond.barrel(xind) : cond.endcap(xind);
0128   };
0129
0130   auto getPair =
0131       [=](EcalLaserAPDPNRatios::EcalLaserAPDPNRatiosMap const& cond) -> EcalLaserAPDPNRatios::EcalLaserAPDPNpair {
0132     return isBarrel ? cond.barrel(xind) : cond.endcap(xind);
0133   };
0134
0135   auto getLinear = [=](EcalLinearCorrections::EcalValueMap const& cond) -> EcalLinearCorrections::Values {
0136     return isBarrel ? cond.barrel(xind) : cond.endcap(xind);
0137   };
0138
0139   apdpnpair = getPair(laserRatiosMap);
0140   linValues = getLinear(linearValueMap);
0141   apdpnref = getCond(laserRefMap);
0142   alpha = getCond(laserAlphaMap);
0143
0144   if (iLM - 1 < (int)laserTimeMap.size()) {
0145     timestamp = laserTimeMap[iLM - 1];
0146   } else {
0147     edm::LogError("EcalLaserDbService") << "error with laserTimeMap!" << endl;
0148     return correctionFactor;
0149   }
0150
0151   if (iLM - 1 < (int)linearTimeMap.size()) {
0152     linTimes = linearTimeMap[iLM - 1];
0153   } else {
0154     edm::LogError("EcalLaserDbService") << "error with laserTimeMap!" << endl;
0155     return correctionFactor;
0156   }
0157
0158 #endif
0159
0160   // should implement some default in case of error...
0161
0162   // should do some quality checks first
0163   // ...
0164
0165   // we will need to treat time differently...
0166   // is time in DB same format as in MC?  probably not...
0167
0168   // interpolation
0169
0170   edm::TimeValue_t t = iTime.value();
0171   long long t_i = 0, t_f = 0;
0172   float p_i = 0, p_f = 0;
0173   long long lt_i = 0, lt_f = 0;
0174   float lp_i = 0, lp_f = 0;
0175
0176   if (t >= timestamp.t1.value() && t < timestamp.t2.value()) {
0177     t_i = timestamp.t1.value();
0178     t_f = timestamp.t2.value();
0179     p_i = apdpnpair.p1;
0180     p_f = apdpnpair.p2;
0181   } else if (t >= timestamp.t2.value() && t <= timestamp.t3.value()) {
0182     t_i = timestamp.t2.value();
0183     t_f = timestamp.t3.value();
0184     p_i = apdpnpair.p2;
0185     p_f = apdpnpair.p3;
0186   } else if (t < timestamp.t1.value()) {
0187     t_i = timestamp.t1.value();
0188     t_f = timestamp.t2.value();
0189     p_i = apdpnpair.p1;
0190     p_f = apdpnpair.p2;
0191
0192   } else if (t > timestamp.t3.value()) {
0193     t_i = timestamp.t2.value();
0194     t_f = timestamp.t3.value();
0195     p_i = apdpnpair.p2;
0196     p_f = apdpnpair.p3;
0197   }
0198
0199   long long t_laser = t;
0200   if (t > timestamp.t3.value() + maxExtrapolationTime_)
0201     t_laser = ((long long)timestamp.t3.value()) + maxExtrapolationTime_;
0202
0203   if (t >= linTimes.t1.value() && t < linTimes.t2.value()) {
0204     lt_i = linTimes.t1.value();
0205     lt_f = linTimes.t2.value();
0206     lp_i = linValues.p1;
0207     lp_f = linValues.p2;
0208   } else if (t >= linTimes.t2.value() && t <= linTimes.t3.value()) {
0209     lt_i = linTimes.t2.value();
0210     lt_f = linTimes.t3.value();
0211     lp_i = linValues.p2;
0212     lp_f = linValues.p3;
0213   } else if (t < linTimes.t1.value()) {
0214     lt_i = linTimes.t1.value();
0215     lt_f = linTimes.t2.value();
0216     lp_i = linValues.p1;
0217     lp_f = linValues.p2;
0218
0219   } else if (t > linTimes.t3.value()) {
0220     lt_i = linTimes.t2.value();
0221     lt_f = linTimes.t3.value();
0222     lp_i = linValues.p2;
0223     lp_f = linValues.p3;
0224   }
0225
0226   if (apdpnref != 0 && (t_i - t_f) != 0 && (lt_i - lt_f) != 0) {
0227     long long tt = t;  // never subtract two unsigned!
0228     float interpolatedLaserResponse =
0229         p_i / apdpnref + float(t_laser - t_i) * (p_f - p_i) / (apdpnref * float(t_f - t_i));
0230     float interpolatedLinearResponse =
0231         lp_i / apdpnref + float(tt - lt_i) * (lp_f - lp_i) / (apdpnref * float(lt_f - lt_i));  // FIXED BY FC
0232
0233     if (interpolatedLinearResponse > 2.f || interpolatedLinearResponse < 0.1f)
0234       interpolatedLinearResponse = 1.f;
0235     if (interpolatedLaserResponse <= 0.) {
0236       // print message only if it is the first time we see < 0
0237       // on this detid
0238       if (channelsWithInvalidCorrection_.insert(xid.rawId()).second) {
0239         edm::LogError("EcalLaserDbService") << "Interpolated Laser correction <0 for detid " << xid.rawId();
0240       }
0241
0242       return correctionFactor;
0243
0244     } else {
0245       float interpolatedTransparencyResponse = interpolatedLaserResponse / interpolatedLinearResponse;
0246
0247       correctionFactor = 1.f / (std::pow(interpolatedTransparencyResponse, alpha) * interpolatedLinearResponse);
0248     }
0249
0250   } else {
0251     edm::LogError("EcalLaserDbService") << "apdpnref (" << apdpnref << ") "
0252                                         << "or t_i-t_f (" << (t_i - t_f) << " is zero!";
0253     return correctionFactor;
0254   }
0255
0256   return correctionFactor;
0257 }
0258
0259 TYPELOOKUP_DATA_REG(EcalLaserDbService);