plugins/alpaka/EnergyComputationKernels.dev.cc

0001 #include "CondFormats/EcalObjects/interface/EcalChannelStatusCode.h"
0002 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0003 #include "DataFormats/EcalDigi/interface/EcalConstants.h"
0004 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0005 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0006
0007 #include "EnergyComputationKernels.h"
0008
0009 namespace ALPAKA_ACCELERATOR_NAMESPACE::ecal::rechit {
0010
0011   ALPAKA_FN_ACC void makeRecHit(int const inputCh,
0012                                 uint32_t const* didCh,
0013                                 float const* amplitude,
0014                                 float const* amplitudeError,
0015                                 float const* jitter,
0016                                 uint32_t const* aux,
0017                                 float const* chi2_in,
0018                                 uint32_t const* flags_in,
0019                                 uint32_t* did,
0020                                 float* energy,
0021                                 float* time,
0022                                 uint32_t* flagBits,
0023                                 uint32_t* extra,
0024                                 EcalRecHitConditionsDevice::ConstView conditionsDev,
0025                                 EcalRecHitParametersDevice::Product const* parametersDev,
0026                                 // time, used for time dependent corrections
0027                                 edm::TimeValue_t const& eventTime,
0028                                 // configuration
0029                                 bool const isPhase2,
0030                                 bool const killDeadChannels,
0031                                 bool const recoverIsolatedChannels,
0032                                 bool const recoverVFE,
0033                                 bool const recoverFE,
0034                                 float const laserMIN,
0035                                 float const laserMAX,
0036                                 uint32_t flagmask) {
0037     // simple copy of input det id to output
0038     did[inputCh] = didCh[inputCh];
0039
0040     auto const did_to_use = DetId{didCh[inputCh]};
0041
0042     auto const isBarrel = did_to_use.subdetId() == EcalBarrel;
0043     auto const hashedId = isBarrel ? ecal::reconstruction::hashedIndexEB(did_to_use.rawId())
0044                                    : conditionsDev.offsetEE() + ecal::reconstruction::hashedIndexEE(did_to_use.rawId());
0045
0046     auto const intercalib = conditionsDev.intercalibConstants()[hashedId];
0047
0048     // only two values for ADC to GeV, EB or EE
0049     auto const adc2gev_to_use = isBarrel ? conditionsDev.adcToGeVConstantEB() : conditionsDev.adcToGeVConstantEE();
0050
0051     auto const timeCalib = conditionsDev.timeCalibConstants()[hashedId];
0052     auto const timeOffset = isBarrel ? conditionsDev.timeOffsetConstantEB() : conditionsDev.timeOffsetConstantEE();
0053
0054     int iLM = 1;
0055     if (isBarrel) {
0056       iLM = ecal::reconstruction::laserMonitoringRegionEB(did_to_use.rawId());
0057     } else {
0058       iLM = ecal::reconstruction::laserMonitoringRegionEE(did_to_use.rawId());
0059     }
0060
0061     long long t_i = 0, t_f = 0;
0062     float p_i = 0, p_f = 0;
0063     auto const t1 = conditionsDev.laserAPDPNRatios_t1()[iLM - 1];
0064     auto const t2 = conditionsDev.laserAPDPNRatios_t2()[iLM - 1];
0065     auto const t3 = conditionsDev.laserAPDPNRatios_t3()[iLM - 1];
0066     auto const p1 = conditionsDev.laserAPDPNRatios_p1()[hashedId];
0067     auto const p2 = conditionsDev.laserAPDPNRatios_p2()[hashedId];
0068     auto const p3 = conditionsDev.laserAPDPNRatios_p3()[hashedId];
0069
0070     // laser
0071     if (eventTime >= t1 && eventTime < t2) {
0072       t_i = t1;
0073       t_f = t2;
0074       p_i = p1;
0075       p_f = p2;
0076     } else if (eventTime >= t2 && eventTime <= t3) {
0077       t_i = t2;
0078       t_f = t3;
0079       p_i = p2;
0080       p_f = p3;
0081     } else if (eventTime < t1) {
0082       t_i = t1;
0083       t_f = t2;
0084       p_i = p1;
0085       p_f = p2;
0086     } else if (eventTime > t3) {
0087       t_i = t2;
0088       t_f = t3;
0089       p_i = p2;
0090       p_f = p3;
0091     }
0092
0093     long long lt_i = 0, lt_f = 0;
0094     float lp_i = 0, lp_f = 0;
0095     auto const lt1 = conditionsDev.linearCorrections_t1()[iLM - 1];
0096     auto const lt2 = conditionsDev.linearCorrections_t2()[iLM - 1];
0097     auto const lt3 = conditionsDev.linearCorrections_t3()[iLM - 1];
0098     auto const lp1 = conditionsDev.linearCorrections_p1()[hashedId];
0099     auto const lp2 = conditionsDev.linearCorrections_p2()[hashedId];
0100     auto const lp3 = conditionsDev.linearCorrections_p3()[hashedId];
0101
0102     // linear corrections
0103     if (eventTime >= lt1 && eventTime < lt2) {
0104       lt_i = lt1;
0105       lt_f = lt2;
0106       lp_i = lp1;
0107       lp_f = lp2;
0108     } else if (eventTime >= lt2 && eventTime <= lt3) {
0109       lt_i = lt2;
0110       lt_f = lt3;
0111       lp_i = lp2;
0112       lp_f = lp3;
0113     } else if (eventTime < lt1) {
0114       lt_i = lt1;
0115       lt_f = lt2;
0116       lp_i = lp1;
0117       lp_f = lp2;
0118     } else if (eventTime > lt3) {
0119       lt_i = lt2;
0120       lt_f = lt3;
0121       lp_i = lp2;
0122       lp_f = lp3;
0123     }
0124
0125     // apdpnref and alpha
0126     auto const apdpnref = conditionsDev.laserAPDPNref()[hashedId];
0127     auto const alpha = conditionsDev.laserAlpha()[hashedId];
0128
0129     // now calculate transparency correction
0130     float lasercalib = 1.;
0131     if (apdpnref != 0 && (t_i - t_f) != 0 && (lt_i - lt_f) != 0) {
0132       long long tt = eventTime;  // never subtract two unsigned!
0133       auto const interpolatedLaserResponse =
0134           p_i / apdpnref + static_cast<float>(tt - t_i) * (p_f - p_i) / (apdpnref * static_cast<float>(t_f - t_i));
0135
0136       auto interpolatedLinearResponse =
0137           lp_i / apdpnref +
0138           static_cast<float>(tt - lt_i) * (lp_f - lp_i) / (apdpnref * static_cast<float>(lt_f - lt_i));  // FIXED BY FC
0139
0140       if (interpolatedLinearResponse > 2.f || interpolatedLinearResponse < 0.1f) {
0141         interpolatedLinearResponse = 1.f;
0142       }
0143       if (interpolatedLaserResponse <= 0.) {
0144         lasercalib = 1.;
0145       } else {
0146         auto const interpolatedTransparencyResponse = interpolatedLaserResponse / interpolatedLinearResponse;
0147         lasercalib = 1.f / (std::pow(interpolatedTransparencyResponse, alpha) * interpolatedLinearResponse);
0148       }
0149     }
0150
0151     //
0152     // Check for channels to be excluded from reconstruction
0153     //
0154     // Default energy not to be updated if "channelStatusToBeExcluded"
0155     energy[inputCh] = -1.;  //un-physical default
0156
0157     // default values for the flags
0158     flagBits[inputCh] = 0;
0159     extra[inputCh] = 0;
0160
0161     auto const dbChStatus = static_cast<EcalChannelStatusCode::Code>(conditionsDev.channelStatus()[hashedId] &
0162                                                                      EcalChannelStatusCode::chStatusMask);
0163     auto const& exclChStatCodes = parametersDev->channelStatusCodesToBeExcluded;
0164     if (exclChStatCodes[dbChStatus]) {
0165       // skip the channel if the channel status bit is set and should be excluded
0166       return;
0167     }
0168
0169     // Take our association map of dbChStatuses-> recHit flagbits and return the appropriate flagbit word
0170     auto const& recoFlagBits = parametersDev->recoFlagBits[dbChStatus];
0171     flagBits[inputCh] = recoFlagBits;
0172
0173     if ((flagmask & recoFlagBits) && killDeadChannels) {
0174       // skip this channel
0175       return;
0176     }
0177
0178     //
0179     // multiply the adc counts with factors to get the energy in GeV
0180     //
0181     energy[inputCh] = amplitude[inputCh] * adc2gev_to_use * intercalib * lasercalib;
0182
0183     // time
0184     auto const sampPeriod = isPhase2 ? ecalPh2::Samp_Period : ecalPh1::Samp_Period;
0185     time[inputCh] = jitter[inputCh] * sampPeriod + timeCalib + timeOffset;
0186
0187     // NB: calculate the "flagBits extra"  --> not really "flags", but actually an encoded version of energy uncertainty, time unc., ...
0188
0189     //
0190     // extra packing ...
0191     //
0192     uint32_t extravar = 0;
0193
0194     // chi2
0195     // truncate
0196     auto const chi2 = chi2_in[inputCh] > 64 ? 64 : chi2_in[inputCh];
0197     // use 7 bits
0198     uint32_t const rawChi2 = lround(chi2 / 64.f * ((1 << 7) - 1));
0199     extravar = ecal::reconstruction::rechitSetMasked(extravar, rawChi2, 0, 7);
0200
0201     // energy uncertainty (amplitudeError is currently not set in the portable uncalibrated rec hit producer)
0202     auto const energyError = amplitudeError[inputCh] * adc2gev_to_use * intercalib * lasercalib;
0203     uint32_t rawEnergy = 0;
0204     if (energyError > 0.001) {
0205       uint16_t const exponent = static_cast<uint16_t>(ecal::reconstruction::rechitGetPower10(energyError));
0206       uint16_t significand = lround(energyError * ip10[exponent]);
0207       // use 13 bits (3 exponent, 10 significand)
0208       rawEnergy = exponent << 10 | significand;
0209     }
0210     extravar = ecal::reconstruction::rechitSetMasked(extravar, rawEnergy, 8, 13);
0211
0212     // time uncertainty directly from uncalib rechit (the jitter error is currently not stored in aux in the portable uncalibrated rec hit producer)
0213     uint8_t const timeErrBits = aux[inputCh] & 0xFF;
0214     extravar = ecal::reconstruction::rechitSetMasked(extravar, timeErrBits, 24, 8);
0215
0216     //
0217     // set output extra variable
0218     //
0219     extra[inputCh] = extravar;
0220
0221     //
0222     // additional flags setting
0223     // using correctly the flags as calculated at the UncalibRecHit stage
0224     //
0225     bool good = true;
0226     if (checkUncalibRecHitFlag(flags_in[inputCh], EcalUncalibratedRecHit::Flags::kLeadingEdgeRecovered)) {
0227       setFlag(flagBits[inputCh], EcalRecHit::Flags::kLeadingEdgeRecovered);
0228       good = false;
0229     }
0230     if (checkUncalibRecHitFlag(flags_in[inputCh], EcalUncalibratedRecHit::Flags::kSaturated)) {
0231       // leading edge recovery failed - still keep the information
0232       // about the saturation and do not flag as dead
0233       setFlag(flagBits[inputCh], EcalRecHit::Flags::kSaturated);
0234       good = false;
0235     }
0236     if (checkUncalibRecHitFlag(flags_in[inputCh], EcalUncalibratedRecHit::Flags::kOutOfTime)) {
0237       setFlag(flagBits[inputCh], EcalRecHit::Flags::kOutOfTime);
0238       good = false;
0239     }
0240     if (checkUncalibRecHitFlag(flags_in[inputCh], EcalUncalibratedRecHit::Flags::kPoorReco)) {
0241       setFlag(flagBits[inputCh], EcalRecHit::Flags::kPoorReco);
0242       good = false;
0243     }
0244     if (checkUncalibRecHitFlag(flags_in[inputCh], EcalUncalibratedRecHit::Flags::kHasSwitchToGain6)) {
0245       setFlag(flagBits[inputCh], EcalRecHit::Flags::kHasSwitchToGain6);
0246     }
0247     if (checkUncalibRecHitFlag(flags_in[inputCh], EcalUncalibratedRecHit::Flags::kHasSwitchToGain1)) {
0248       setFlag(flagBits[inputCh], EcalRecHit::Flags::kHasSwitchToGain1);
0249     }
0250
0251     if (good) {
0252       setFlag(flagBits[inputCh], EcalRecHit::Flags::kGood);
0253     }
0254
0255     if (lasercalib < laserMIN || lasercalib > laserMAX) {
0256       setFlag(flagBits[inputCh], EcalRecHit::Flags::kPoorCalib);
0257     }
0258
0259     // recover, killing, and other stuff
0260     //
0261     // Structure:
0262     //  EB
0263     //  EE
0264     //
0265     //  - single MVA
0266     //  - democratic sharing
0267     //  - kill all the other cases
0268
0269     // recoverable channel status codes
0270     if (dbChStatus == EcalChannelStatusCode::Code::kFixedG0 ||
0271         dbChStatus == EcalChannelStatusCode::Code::kNonRespondingIsolated ||
0272         dbChStatus == EcalChannelStatusCode::Code::kDeadVFE) {
0273       bool is_Single = false;
0274       bool is_FE = false;
0275       bool is_VFE = false;
0276
0277       if (dbChStatus == EcalChannelStatusCode::Code::kDeadVFE) {
0278         is_VFE = true;
0279       } else if (dbChStatus == EcalChannelStatusCode::Code::kDeadFE) {
0280         is_FE = true;
0281       } else {
0282         is_Single = true;
0283       }
0284
0285       // EB
0286       if (isBarrel) {
0287         if (is_Single || is_FE || is_VFE) {
0288           if (is_Single && (recoverIsolatedChannels || !killDeadChannels)) {
0289             // single MVA
0290             // TODO
0291           } else if (is_FE && (recoverFE || !killDeadChannels)) {
0292             // democratic sharing
0293             // TODO
0294           } else {
0295             // kill all the other cases
0296             energy[inputCh] = 0.;
0297             // TODO: set flags
0298           }
0299         }
0300       } else {  // EE
0301         if (is_Single || is_FE || is_VFE) {
0302           if (is_Single && (recoverIsolatedChannels || !killDeadChannels)) {
0303             // single MVA
0304             // TODO
0305           } else if (is_FE && (recoverFE || !killDeadChannels)) {
0306             // democratic sharing
0307             // TODO
0308           } else {
0309             // kill all the other cases
0310             energy[inputCh] = 0.;
0311             // TODO: set flags
0312           }
0313         }
0314       }
0315     }
0316   }
0317
0318 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::ecal::rechit