Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:23:51 

0001 #include "SimCalorimetry/HcalTrigPrimAlgos/interface/HcalTriggerPrimitiveAlgo.h"
0002 
0003 #include "CalibFormats/CaloObjects/interface/IntegerCaloSamples.h"
0004 #include "CondFormats/HcalObjects/interface/HcalTPParameters.h"
0005 #include "CondFormats/HcalObjects/interface/HcalTPChannelParameters.h"
0006 
0007 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0008 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
0009 #include "DataFormats/HcalDetId/interface/HcalElectronicsId.h"
0010 
0011 #include "EventFilter/HcalRawToDigi/interface/HcalDCCHeader.h"
0012 #include "EventFilter/HcalRawToDigi/interface/HcalHTRData.h"
0013 
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 
0017 #include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
0018 
0019 #include <iostream>
0020 
0021 using namespace std;
0022 
0023 HcalTriggerPrimitiveAlgo::HcalTriggerPrimitiveAlgo(bool pf,
0024                                                    const std::vector<double>& w,
0025                                                    int latency,
0026                                                    uint32_t FG_threshold,
0027                                                    const std::vector<uint32_t>& FG_HF_thresholds,
0028                                                    uint32_t ZS_threshold,
0029                                                    int numberOfSamples,
0030                                                    int numberOfPresamples,
0031                                                    int numberOfFilterPresamplesHBQIE11,
0032                                                    int numberOfFilterPresamplesHEQIE11,
0033                                                    int numberOfSamplesHF,
0034                                                    int numberOfPresamplesHF,
0035                                                    bool useTDCInMinBiasBits,
0036                                                    uint32_t minSignalThreshold,
0037                                                    uint32_t PMT_NoiseThreshold)
0038     : incoder_(nullptr),
0039       outcoder_(nullptr),
0040       theThreshold(0),
0041       peakfind_(pf),
0042       weights_(w),
0043       latency_(latency),
0044       FG_threshold_(FG_threshold),
0045       FG_HF_thresholds_(FG_HF_thresholds),
0046       ZS_threshold_(ZS_threshold),
0047       numberOfSamples_(numberOfSamples),
0048       numberOfPresamples_(numberOfPresamples),
0049       numberOfFilterPresamplesHBQIE11_(numberOfFilterPresamplesHBQIE11),
0050       numberOfFilterPresamplesHEQIE11_(numberOfFilterPresamplesHEQIE11),
0051       numberOfSamplesHF_(numberOfSamplesHF),
0052       numberOfPresamplesHF_(numberOfPresamplesHF),
0053       useTDCInMinBiasBits_(useTDCInMinBiasBits),
0054       minSignalThreshold_(minSignalThreshold),
0055       PMT_NoiseThreshold_(PMT_NoiseThreshold),
0056       NCTScaleShift(0),
0057       RCTScaleShift(0),
0058       peak_finder_algorithm_(2),
0059       override_parameters_() {
0060   //No peak finding setting (for Fastsim)
0061   if (!peakfind_) {
0062     numberOfSamples_ = 1;
0063     numberOfPresamples_ = 0;
0064     numberOfSamplesHF_ = 1;
0065     numberOfPresamplesHF_ = 0;
0066   }
0067   // Switch to integer for comparisons - remove compiler warning
0068   ZS_threshold_I_ = ZS_threshold_;
0069 }
0070 
0071 HcalTriggerPrimitiveAlgo::~HcalTriggerPrimitiveAlgo() {}
0072 
0073 void HcalTriggerPrimitiveAlgo::setUpgradeFlags(bool hb, bool he, bool hf) {
0074   upgrade_hb_ = hb;
0075   upgrade_he_ = he;
0076   upgrade_hf_ = hf;
0077 }
0078 
0079 void HcalTriggerPrimitiveAlgo::setFixSaturationFlag(bool fix_saturation) { fix_saturation_ = fix_saturation; }
0080 
0081 void HcalTriggerPrimitiveAlgo::overrideParameters(const edm::ParameterSet& ps) {
0082   override_parameters_ = ps;
0083 
0084   if (override_parameters_.exists("ADCThresholdHF")) {
0085     override_adc_hf_ = true;
0086     override_adc_hf_value_ = override_parameters_.getParameter<uint32_t>("ADCThresholdHF");
0087   }
0088   if (override_parameters_.exists("TDCMaskHF")) {
0089     override_tdc_hf_ = true;
0090     override_tdc_hf_value_ = override_parameters_.getParameter<unsigned long long>("TDCMaskHF");
0091   }
0092 }
0093 
0094 void HcalTriggerPrimitiveAlgo::addSignal(const HBHEDataFrame& frame) {
0095   // TODO: Need to add support for seperate 28, 29 in HE
0096   //Hack for 300_pre10, should be removed.
0097   if (frame.id().depth() == 5)
0098     return;
0099 
0100   std::vector<HcalTrigTowerDetId> ids = theTrigTowerGeometry->towerIds(frame.id());
0101   assert(ids.size() == 1 || ids.size() == 2);
0102   IntegerCaloSamples samples1(ids[0], int(frame.size()));
0103 
0104   samples1.setPresamples(frame.presamples());
0105   incoder_->adc2Linear(frame, samples1);
0106 
0107   std::vector<bool> msb;
0108   incoder_->lookupMSB(frame, msb);
0109 
0110   if (ids.size() == 2) {
0111     // make a second trigprim for the other one, and split the energy
0112     IntegerCaloSamples samples2(ids[1], samples1.size());
0113     for (int i = 0; i < samples1.size(); ++i) {
0114       samples1[i] = uint32_t(samples1[i] * 0.5);
0115       samples2[i] = samples1[i];
0116     }
0117     samples2.setPresamples(frame.presamples());
0118     addSignal(samples2);
0119     addFG(ids[1], msb);
0120   }
0121   addSignal(samples1);
0122   addFG(ids[0], msb);
0123 }
0124 
0125 void HcalTriggerPrimitiveAlgo::addSignal(const HFDataFrame& frame) {
0126   if (frame.id().depth() == 1 || frame.id().depth() == 2) {
0127     std::vector<HcalTrigTowerDetId> ids = theTrigTowerGeometry->towerIds(frame.id());
0128     std::vector<HcalTrigTowerDetId>::const_iterator it;
0129     for (it = ids.begin(); it != ids.end(); ++it) {
0130       HcalTrigTowerDetId trig_tower_id = *it;
0131       IntegerCaloSamples samples(trig_tower_id, frame.size());
0132       samples.setPresamples(frame.presamples());
0133       incoder_->adc2Linear(frame, samples);
0134 
0135       // Don't add to final collection yet
0136       // HF PMT veto sum is calculated in analyzerHF()
0137       IntegerCaloSamples zero_samples(trig_tower_id, frame.size());
0138       zero_samples.setPresamples(frame.presamples());
0139       addSignal(zero_samples);
0140 
0141       // Pre-LS1 Configuration
0142       if (trig_tower_id.version() == 0) {
0143         // Mask off depths: fgid is the same for both depths
0144         uint32_t fgid = (frame.id().maskDepth());
0145 
0146         if (theTowerMapFGSum.find(trig_tower_id) == theTowerMapFGSum.end()) {
0147           SumFGContainer sumFG;
0148           theTowerMapFGSum.insert(std::pair<HcalTrigTowerDetId, SumFGContainer>(trig_tower_id, sumFG));
0149         }
0150 
0151         SumFGContainer& sumFG = theTowerMapFGSum[trig_tower_id];
0152         SumFGContainer::iterator sumFGItr;
0153         for (sumFGItr = sumFG.begin(); sumFGItr != sumFG.end(); ++sumFGItr) {
0154           if (sumFGItr->id() == fgid) {
0155             break;
0156           }
0157         }
0158         // If find
0159         if (sumFGItr != sumFG.end()) {
0160           for (int i = 0; i < samples.size(); ++i) {
0161             (*sumFGItr)[i] += samples[i];
0162           }
0163         } else {
0164           //Copy samples (change to fgid)
0165           IntegerCaloSamples sumFGSamples(DetId(fgid), samples.size());
0166           sumFGSamples.setPresamples(samples.presamples());
0167           for (int i = 0; i < samples.size(); ++i) {
0168             sumFGSamples[i] = samples[i];
0169           }
0170           sumFG.push_back(sumFGSamples);
0171         }
0172 
0173         // set veto to true if Long or Short less than threshold
0174         if (HF_Veto.find(fgid) == HF_Veto.end()) {
0175           vector<bool> vetoBits(samples.size(), false);
0176           HF_Veto[fgid] = vetoBits;
0177         }
0178         for (int i = 0; i < samples.size(); ++i) {
0179           if (samples[i] < minSignalThreshold_) {
0180             HF_Veto[fgid][i] = true;
0181           }
0182         }
0183       }
0184       // HF 1x1
0185       else if (trig_tower_id.version() == 1) {
0186         uint32_t fgid = (frame.id().maskDepth());
0187         HFDetails& details = theHFDetailMap[trig_tower_id][fgid];
0188         // Check the frame type to determine long vs short
0189         if (frame.id().depth() == 1) {  // Long
0190           details.long_fiber = samples;
0191           details.LongDigi = frame;
0192         } else if (frame.id().depth() == 2) {  // Short
0193           details.short_fiber = samples;
0194           details.ShortDigi = frame;
0195         } else {
0196           // Neither long nor short... So we have no idea what to do
0197           edm::LogWarning("HcalTPAlgo") << "Unable to figure out what to do with data frame for " << frame.id();
0198           return;
0199         }
0200       }
0201       // Uh oh, we are in a bad/unknown state! Things will start crashing.
0202       else {
0203         return;
0204       }
0205     }
0206   }
0207 }
0208 
0209 void HcalTriggerPrimitiveAlgo::addSignal(const QIE10DataFrame& frame) {
0210   HcalDetId detId = frame.detid();
0211   // prevent QIE10 calibration channels from entering TP emulation
0212   if (detId.subdet() != HcalForward)
0213     return;
0214 
0215   auto ids = theTrigTowerGeometry->towerIds(frame.id());
0216   for (const auto& id : ids) {
0217     if (id.version() == 0) {
0218       edm::LogError("HcalTPAlgo") << "Encountered QIE10 data frame mapped to TP version 0:" << id;
0219       continue;
0220     }
0221 
0222     int nsamples = frame.samples();
0223 
0224     IntegerCaloSamples samples(id, nsamples);
0225     samples.setPresamples(frame.presamples());
0226     incoder_->adc2Linear(frame, samples);
0227 
0228     // Don't add to final collection yet
0229     // HF PMT veto sum is calculated in analyzerHF()
0230     IntegerCaloSamples zero_samples(id, nsamples);
0231     zero_samples.setPresamples(frame.presamples());
0232     addSignal(zero_samples);
0233 
0234     auto fid = HcalDetId(frame.id());
0235     auto& details = theHFUpgradeDetailMap[id][fid.maskDepth()];
0236     auto& detail = details[fid.depth() - 1];
0237     detail.samples = samples;
0238     detail.digi = frame;
0239     detail.validity.resize(nsamples);
0240     detail.passTDC.resize(nsamples);
0241     incoder_->lookupMSB(frame, detail.fgbits);
0242     for (int idx = 0; idx < nsamples; ++idx) {
0243       detail.validity[idx] = validChannel(frame, idx);
0244       detail.passTDC[idx] = passTDC(frame, idx);
0245     }
0246   }
0247 }
0248 
0249 void HcalTriggerPrimitiveAlgo::addSignal(const QIE11DataFrame& frame) {
0250   HcalDetId detId(frame.id());
0251   // prevent QIE11 calibration channels from entering TP emulation
0252   if (detId.subdet() != HcalEndcap && detId.subdet() != HcalBarrel)
0253     return;
0254 
0255   std::vector<HcalTrigTowerDetId> ids = theTrigTowerGeometry->towerIds(detId);
0256   assert(ids.size() == 1 || ids.size() == 2);
0257   IntegerCaloSamples samples1(ids[0], int(frame.samples()));
0258 
0259   samples1.setPresamples(frame.presamples());
0260   incoder_->adc2Linear(frame, samples1);
0261 
0262   std::vector<std::bitset<2>> msb(frame.samples(), 0);
0263   incoder_->lookupMSB(frame, msb);
0264 
0265   if (ids.size() == 2) {
0266     // make a second trigprim for the other one, and share the energy
0267     IntegerCaloSamples samples2(ids[1], samples1.size());
0268     for (int i = 0; i < samples1.size(); ++i) {
0269       samples1[i] = uint32_t(samples1[i]);
0270       samples2[i] = samples1[i];
0271     }
0272     samples2.setPresamples(frame.presamples());
0273     addSignal(samples2);
0274     addUpgradeFG(ids[1], detId.depth(), msb);
0275     addUpgradeTDCFG(ids[1], frame);
0276   }
0277   addSignal(samples1);
0278   addUpgradeFG(ids[0], detId.depth(), msb);
0279   addUpgradeTDCFG(ids[0], frame);
0280 }
0281 
0282 void HcalTriggerPrimitiveAlgo::addSignal(const IntegerCaloSamples& samples) {
0283   HcalTrigTowerDetId id(samples.id());
0284   SumMap::iterator itr = theSumMap.find(id);
0285 
0286   if (itr == theSumMap.end()) {
0287     theSumMap.insert(std::make_pair(id, samples));
0288   } else {
0289     // wish CaloSamples had a +=
0290     for (int i = 0; i < samples.size(); ++i) {
0291       (itr->second)[i] += samples[i];
0292     }
0293   }
0294 
0295   // if fix_saturation == true, keep track of tower with saturated input LUT
0296   if (fix_saturation_) {
0297     SatMap::iterator itr_sat = theSatMap.find(id);
0298 
0299     assert((itr == theSumMap.end()) == (itr_sat == theSatMap.end()));
0300 
0301     if (itr_sat == theSatMap.end()) {
0302       vector<bool> check_sat;
0303       for (int i = 0; i < samples.size(); ++i) {
0304         if (!(samples[i] < QIE11_LINEARIZATION_ET)) {
0305           check_sat.push_back(true);
0306         } else
0307           check_sat.push_back(false);
0308       }
0309       theSatMap.insert(std::make_pair(id, check_sat));
0310     } else {
0311       for (int i = 0; i < samples.size(); ++i) {
0312         if (!(samples[i] < QIE11_LINEARIZATION_ET))
0313           (itr_sat->second)[i] = true;
0314       }
0315     }
0316   }
0317 }
0318 
0319 void HcalTriggerPrimitiveAlgo::analyze(IntegerCaloSamples& samples, HcalTriggerPrimitiveDigi& result) {
0320   int shrink = weights_.size() - 1;
0321   std::vector<bool>& msb = fgMap_[samples.id()];
0322   IntegerCaloSamples sum(samples.id(), samples.size());
0323 
0324   //slide algo window
0325   for (int ibin = 0; ibin < int(samples.size()) - shrink; ++ibin) {
0326     int algosumvalue = 0;
0327     for (unsigned int i = 0; i < weights_.size(); i++) {
0328       //add up value * scale factor
0329       algosumvalue += int(samples[ibin + i] * weights_[i]);
0330     }
0331     if (algosumvalue < 0)
0332       sum[ibin] = 0;  // low-side
0333                       //high-side
0334     //else if (algosumvalue>QIE8_LINEARIZATION_ET) sum[ibin]=QIE8_LINEARIZATION_ET;
0335     else
0336       sum[ibin] = algosumvalue;  //assign value to sum[]
0337   }
0338 
0339   // Align digis and TP
0340   int dgPresamples = samples.presamples();
0341   int tpPresamples = numberOfPresamples_;
0342   int shift = dgPresamples - tpPresamples;
0343   int dgSamples = samples.size();
0344   int tpSamples = numberOfSamples_;
0345   if (peakfind_) {
0346     if ((shift < shrink) || (shift + tpSamples + shrink > dgSamples - (peak_finder_algorithm_ - 1))) {
0347       edm::LogInfo("HcalTriggerPrimitiveAlgo::analyze")
0348           << "TP presample or size from the configuration file is out of the accessible range. Using digi values from "
0349              "data instead...";
0350       shift = shrink;
0351       tpPresamples = dgPresamples - shrink;
0352       tpSamples = dgSamples - (peak_finder_algorithm_ - 1) - shrink - shift;
0353     }
0354   }
0355 
0356   std::vector<int> finegrain(tpSamples, false);
0357 
0358   IntegerCaloSamples output(samples.id(), tpSamples);
0359   output.setPresamples(tpPresamples);
0360 
0361   for (int ibin = 0; ibin < tpSamples; ++ibin) {
0362     // ibin - index for output TP
0363     // idx - index for samples + shift
0364     int idx = ibin + shift;
0365 
0366     //Peak finding
0367     if (peakfind_) {
0368       bool isPeak = false;
0369       switch (peak_finder_algorithm_) {
0370         case 1:
0371           isPeak = (samples[idx] > samples[idx - 1] && samples[idx] >= samples[idx + 1] && samples[idx] > theThreshold);
0372           break;
0373         case 2:
0374           isPeak = (sum[idx] > sum[idx - 1] && sum[idx] >= sum[idx + 1] && sum[idx] > theThreshold);
0375           break;
0376         default:
0377           break;
0378       }
0379 
0380       if (isPeak) {
0381         output[ibin] = std::min<unsigned int>(sum[idx], QIE8_LINEARIZATION_ET);
0382         finegrain[ibin] = msb[idx];
0383       }
0384       // Not a peak
0385       else
0386         output[ibin] = 0;
0387     } else {  // No peak finding, just output running sum
0388       output[ibin] = std::min<unsigned int>(sum[idx], QIE8_LINEARIZATION_ET);
0389       finegrain[ibin] = msb[idx];
0390     }
0391 
0392     // Only Pegged for 1-TS algo.
0393     if (peak_finder_algorithm_ == 1) {
0394       if (samples[idx] >= QIE8_LINEARIZATION_ET)
0395         output[ibin] = QIE8_LINEARIZATION_ET;
0396     }
0397   }
0398   outcoder_->compress(output, finegrain, result);
0399 }
0400 
0401 void HcalTriggerPrimitiveAlgo::analyzeQIE11(IntegerCaloSamples& samples,
0402                                             vector<bool> sample_saturation,
0403                                             HcalTriggerPrimitiveDigi& result,
0404                                             const HcalFinegrainBit& fg_algo) {
0405   HcalDetId detId(samples.id());
0406 
0407   // Get the |ieta| for current sample
0408   int theIeta = detId.ietaAbs();
0409 
0410   unsigned int dgSamples = samples.size();
0411   unsigned int dgPresamples = samples.presamples();
0412 
0413   unsigned int tpSamples = numberOfSamples_;
0414   unsigned int tpPresamples = numberOfPresamples_;
0415 
0416   unsigned int filterSamples = weightsQIE11_[theIeta].size();
0417   unsigned int filterPresamples = theIeta > theTrigTowerGeometry->topology().lastHBRing()
0418                                       ? numberOfFilterPresamplesHEQIE11_
0419                                       : numberOfFilterPresamplesHBQIE11_;
0420 
0421   unsigned int shift = dgPresamples - tpPresamples;
0422 
0423   // shrink keeps the FIR filter from going off the end of the 8TS vector
0424   unsigned int shrink = filterSamples - 1;
0425 
0426   auto& msb = fgUpgradeMap_[samples.id()];
0427   auto& timingTDC = fgUpgradeTDCMap_[samples.id()];
0428   IntegerCaloSamples sum(samples.id(), samples.size());
0429 
0430   std::vector<HcalTrigTowerDetId> ids = theTrigTowerGeometry->towerIds(detId);
0431 
0432   // keep track of tower with saturated energy and force the total TP saturated
0433   bool force_saturation[samples.size()];
0434   for (int i = 0; i < samples.size(); i++) {
0435     force_saturation[i] = false;
0436   }
0437 
0438   //slide algo window
0439   for (unsigned int ibin = 0; ibin < dgSamples - shrink; ++ibin) {
0440     int algosumvalue = 0;
0441     bool check_sat = false;
0442     //TP energy calculation for PFA2
0443     if (weightsQIE11_[theIeta][0] == 255) {
0444       for (unsigned int i = 0; i < filterSamples; i++) {
0445         //add up value * scale factor
0446         // In addition, divide by two in the 10 degree phi segmentation region
0447         // to mimic 5 degree segmentation for the trigger
0448         unsigned int sample = samples[ibin + i];
0449 
0450         if (fix_saturation_ && (sample_saturation.size() > ibin + i))
0451           check_sat = (check_sat | sample_saturation[ibin + i] | (sample > QIE11_MAX_LINEARIZATION_ET));
0452 
0453         if (sample > QIE11_MAX_LINEARIZATION_ET)
0454           sample = QIE11_MAX_LINEARIZATION_ET;
0455 
0456         // Usually use a segmentation factor of 1.0 but for ieta >= 21 use 2
0457         int segmentationFactor = 1;
0458         if (ids.size() == 2) {
0459           segmentationFactor = 2;
0460         }
0461 
0462         algosumvalue += int(sample / segmentationFactor);
0463       }
0464       if (algosumvalue < 0)
0465         sum[ibin] = 0;  // low-side
0466                         //high-side
0467       //else if (algosumvalue>QIE11_LINEARIZATION_ET) sum[ibin]=QIE11_LINEARIZATION_ET;
0468       else
0469         sum[ibin] = algosumvalue;  //assign value to sum[]
0470 
0471       if (check_sat)
0472         force_saturation[ibin] = true;
0473       //TP energy calculation for PFA1' and PFA1
0474     } else {
0475       //add up value * scale factor
0476       // In addition, divide by two in the 10 degree phi segmentation region
0477       // to mimic 5 degree segmentation for the trigger
0478       int sampleTS = samples[ibin + 1];
0479       int sampleTSminus1 = samples[ibin];
0480 
0481       if (fix_saturation_ && (sample_saturation.size() > ibin + 1))
0482         check_sat = (sample_saturation[ibin + 1] || (sampleTS >= QIE11_MAX_LINEARIZATION_ET) ||
0483                      sample_saturation[ibin] || (sampleTSminus1 >= QIE11_MAX_LINEARIZATION_ET));
0484 
0485       if (sampleTS > QIE11_MAX_LINEARIZATION_ET)
0486         sampleTS = QIE11_MAX_LINEARIZATION_ET;
0487 
0488       if (sampleTSminus1 > QIE11_MAX_LINEARIZATION_ET)
0489         sampleTSminus1 = QIE11_MAX_LINEARIZATION_ET;
0490 
0491       // Usually use a segmentation factor of 1.0 but for ieta >= 21 use 2
0492       int segmentationFactor = 1;
0493       if (ids.size() == 2) {
0494         segmentationFactor = 2;
0495       }
0496 
0497       // Based on the |ieta| of the sample, retrieve the correct region weight
0498       int theWeight = weightsQIE11_[theIeta][0];
0499 
0500       algosumvalue = ((sampleTS << 8) - (sampleTSminus1 * theWeight)) / 256 / segmentationFactor;
0501 
0502       if (algosumvalue < 0)
0503         sum[ibin] = 0;  // low-side
0504                         //high-side
0505       //else if (algosumvalue>QIE11_LINEARIZATION_ET) sum[ibin]=QIE11_LINEARIZATION_ET;
0506       else
0507         sum[ibin] = algosumvalue;  //assign value to sum[]
0508 
0509       if (check_sat)
0510         force_saturation[ibin] = true;
0511     }
0512   }
0513 
0514   std::vector<int> finegrain(tpSamples, false);
0515 
0516   IntegerCaloSamples output(samples.id(), tpSamples);
0517   output.setPresamples(tpPresamples);
0518 
0519   for (unsigned int ibin = 0; ibin < tpSamples; ++ibin) {
0520     // ibin - index for output TP
0521     // idx - index for samples + shift - filterPresamples
0522     int idx = ibin + shift - filterPresamples;
0523 
0524     // When idx is <= 0 peakfind would compare out-of-bounds of the vector. Avoid this ambiguity
0525     if (idx <= 0) {
0526       output[ibin] = 0;
0527       continue;
0528     }
0529 
0530     //Only run the peak-finder when the PFA2 FIR filter is running, which corresponds to weights = 1
0531     if (weightsQIE11_[theIeta][0] == 255) {
0532       bool isPeak = (sum[idx] > sum[idx - 1] && sum[idx] >= sum[idx + 1] && sum[idx] > theThreshold);
0533       if (isPeak) {
0534         output[ibin] = std::min<unsigned int>(sum[idx], QIE11_MAX_LINEARIZATION_ET);
0535 
0536         if (fix_saturation_ && force_saturation[idx] && ids.size() == 2)
0537           output[ibin] = QIE11_MAX_LINEARIZATION_ET / 2;
0538         else if (fix_saturation_ && force_saturation[idx])
0539           output[ibin] = QIE11_MAX_LINEARIZATION_ET;
0540 
0541       } else {
0542         // Not a peak
0543         output[ibin] = 0;
0544       }
0545     } else {
0546       output[ibin] = std::min<unsigned int>(sum[idx], QIE11_MAX_LINEARIZATION_ET);
0547 
0548       if (fix_saturation_ && force_saturation[idx] && ids.size() == 2)
0549         output[ibin] = QIE11_MAX_LINEARIZATION_ET / 2;
0550       else if (fix_saturation_ && force_saturation[idx])
0551         output[ibin] = QIE11_MAX_LINEARIZATION_ET;
0552     }
0553     // peak-finding is not applied for FG bits
0554     // compute(msb) returns two bits (MIP). compute(timingTDC,ids) returns 6 bits (1 depth, 1 prompt, 1 delayed 01, 1 delayed 10, 2 reserved)
0555     finegrain[ibin] = fg_algo.compute(timingTDC[idx + filterPresamples], ids[0]).to_ulong() |
0556                       fg_algo.compute(msb[idx + filterPresamples]).to_ulong() << 4;
0557     if (ibin == tpPresamples && (idx + filterPresamples) != dgPresamples)
0558       edm::LogError("HcalTriggerPritimveAlgo")
0559           << "TP SOI (tpPresamples = " << tpPresamples
0560           << ") is not aligned with digi SOI (dgPresamples = " << dgPresamples << ")";
0561   }
0562   outcoder_->compress(output, finegrain, result);
0563 }
0564 
0565 void HcalTriggerPrimitiveAlgo::analyzeHF(IntegerCaloSamples& samples,
0566                                          HcalTriggerPrimitiveDigi& result,
0567                                          const int hf_lumi_shift) {
0568   HcalTrigTowerDetId detId(samples.id());
0569 
0570   // Align digis and TP
0571   int dgPresamples = samples.presamples();
0572   int tpPresamples = numberOfPresamplesHF_;
0573   int shift = dgPresamples - tpPresamples;
0574   int dgSamples = samples.size();
0575   int tpSamples = numberOfSamplesHF_;
0576   if (shift < 0 || shift + tpSamples > dgSamples) {
0577     edm::LogInfo("HcalTriggerPrimitiveAlgo::analyzeHF")
0578         << "TP presample or size from the configuration file is out of the accessible range. Using digi values from "
0579            "data instead...";
0580     tpPresamples = dgPresamples;
0581     shift = 0;
0582     tpSamples = dgSamples;
0583   }
0584 
0585   std::vector<int> finegrain(tpSamples, false);
0586 
0587   TowerMapFGSum::const_iterator tower2fg = theTowerMapFGSum.find(detId);
0588   assert(tower2fg != theTowerMapFGSum.end());
0589 
0590   const SumFGContainer& sumFG = tower2fg->second;
0591   // Loop over all L+S pairs that mapped from samples.id()
0592   // Note: 1 samples.id() = 6 x (L+S) without noZS
0593   for (SumFGContainer::const_iterator sumFGItr = sumFG.begin(); sumFGItr != sumFG.end(); ++sumFGItr) {
0594     const std::vector<bool>& veto = HF_Veto[sumFGItr->id().rawId()];
0595     for (int ibin = 0; ibin < tpSamples; ++ibin) {
0596       int idx = ibin + shift;
0597       // if not vetod, add L+S to total sum and calculate FG
0598       bool vetoed = idx < int(veto.size()) && veto[idx];
0599       if (!(vetoed && (*sumFGItr)[idx] > PMT_NoiseThreshold_)) {
0600         samples[idx] += (*sumFGItr)[idx];
0601         finegrain[ibin] = (finegrain[ibin] || (*sumFGItr)[idx] >= FG_threshold_);
0602       }
0603     }
0604   }
0605 
0606   IntegerCaloSamples output(samples.id(), tpSamples);
0607   output.setPresamples(tpPresamples);
0608 
0609   for (int ibin = 0; ibin < tpSamples; ++ibin) {
0610     int idx = ibin + shift;
0611     output[ibin] = samples[idx] >> hf_lumi_shift;
0612     static const int MAX_OUTPUT = QIE8_LINEARIZATION_ET;  // QIE8_LINEARIZATION_ET = 1023
0613     if (output[ibin] > MAX_OUTPUT)
0614       output[ibin] = MAX_OUTPUT;
0615   }
0616   outcoder_->compress(output, finegrain, result);
0617 }
0618 
0619 void HcalTriggerPrimitiveAlgo::analyzeHF2016(const IntegerCaloSamples& samples,
0620                                              HcalTriggerPrimitiveDigi& result,
0621                                              const int hf_lumi_shift,
0622                                              const HcalFeatureBit* embit) {
0623   // Align digis and TP
0624   const int SHIFT = samples.presamples() - numberOfPresamplesHF_;
0625   assert(SHIFT >= 0);
0626   assert((SHIFT + numberOfSamplesHF_) <= samples.size());
0627 
0628   // Try to find the HFDetails from the map corresponding to our samples
0629   const HcalTrigTowerDetId detId(samples.id());
0630   HFDetailMap::const_iterator it = theHFDetailMap.find(detId);
0631   // Missing values will give an empty digi
0632   if (it == theHFDetailMap.end()) {
0633     return;
0634   }
0635 
0636   std::vector<std::bitset<2>> finegrain(numberOfSamplesHF_, false);
0637 
0638   // Set up out output of IntergerCaloSamples
0639   IntegerCaloSamples output(samples.id(), numberOfSamplesHF_);
0640   output.setPresamples(numberOfPresamplesHF_);
0641 
0642   for (const auto& item : it->second) {
0643     auto& details = item.second;
0644     for (int ibin = 0; ibin < numberOfSamplesHF_; ++ibin) {
0645       const int IDX = ibin + SHIFT;
0646       int long_fiber_val = 0;
0647       if (IDX < details.long_fiber.size()) {
0648         long_fiber_val = details.long_fiber[IDX];
0649       }
0650       int short_fiber_val = 0;
0651       if (IDX < details.short_fiber.size()) {
0652         short_fiber_val = details.short_fiber[IDX];
0653       }
0654       output[ibin] += (long_fiber_val + short_fiber_val);
0655 
0656       uint32_t ADCLong = details.LongDigi[ibin].adc();
0657       uint32_t ADCShort = details.ShortDigi[ibin].adc();
0658 
0659       if (details.LongDigi.id().ietaAbs() >= FIRST_FINEGRAIN_TOWER) {
0660         finegrain[ibin][1] = (ADCLong > FG_HF_thresholds_[0] || ADCShort > FG_HF_thresholds_[0]);
0661 
0662         if (embit != nullptr)
0663           finegrain[ibin][0] = embit->fineGrainbit(details.ShortDigi, details.LongDigi, ibin);
0664       }
0665     }
0666   }
0667 
0668   for (int bin = 0; bin < numberOfSamplesHF_; ++bin) {
0669     static const unsigned int MAX_OUTPUT = QIE8_LINEARIZATION_ET;  // QIE8_LINEARIZATION_ET = 1023
0670     output[bin] = min({MAX_OUTPUT, output[bin] >> hf_lumi_shift});
0671   }
0672 
0673   std::vector<int> finegrain_converted;
0674   finegrain_converted.reserve(finegrain.size());
0675   for (const auto& fg : finegrain)
0676     finegrain_converted.push_back(fg.to_ulong());
0677   outcoder_->compress(output, finegrain_converted, result);
0678 }
0679 
0680 bool HcalTriggerPrimitiveAlgo::passTDC(const QIE10DataFrame& digi, int ts) const {
0681   auto parameters = conditions_->getHcalTPParameters();
0682   auto adc_threshold = parameters->getADCThresholdHF();
0683   auto tdc_mask = parameters->getTDCMaskHF();
0684 
0685   if (override_adc_hf_)
0686     adc_threshold = override_adc_hf_value_;
0687   if (override_tdc_hf_)
0688     tdc_mask = override_tdc_hf_value_;
0689 
0690   if (digi[ts].adc() < adc_threshold)
0691     return true;
0692 
0693   return (1ul << digi[ts].le_tdc()) & tdc_mask;
0694 }
0695 
0696 bool HcalTriggerPrimitiveAlgo::validChannel(const QIE10DataFrame& digi, int ts) const {
0697   // channels with invalid data should not contribute to the sum
0698   if (digi.linkError() || ts >= digi.samples() || !digi[ts].ok())
0699     return false;
0700 
0701   auto mask = conditions_->getHcalTPChannelParameter(HcalDetId(digi.id()))->getMask();
0702   if (mask)
0703     return false;
0704 
0705   return true;
0706 }
0707 
0708 void HcalTriggerPrimitiveAlgo::analyzeHFQIE10(const IntegerCaloSamples& samples,
0709                                               HcalTriggerPrimitiveDigi& result,
0710                                               const int hf_lumi_shift,
0711                                               const HcalFeatureBit* embit) {
0712   // Align digis and TP
0713   const int shift = samples.presamples() - numberOfPresamplesHF_;
0714   assert(shift >= 0);
0715   assert((shift + numberOfSamplesHF_) <= samples.size());
0716   assert(hf_lumi_shift >= 2);
0717 
0718   // Try to find the HFDetails from the map corresponding to our samples
0719   const HcalTrigTowerDetId detId(samples.id());
0720   auto it = theHFUpgradeDetailMap.find(detId);
0721   // Missing values will give an empty digi
0722   if (it == theHFUpgradeDetailMap.end()) {
0723     return;
0724   }
0725 
0726   std::vector<std::bitset<2>> finegrain(numberOfSamplesHF_, false);
0727 
0728   // Set up out output of IntergerCaloSamples
0729   IntegerCaloSamples output(samples.id(), numberOfSamplesHF_);
0730   output.setPresamples(numberOfPresamplesHF_);
0731 
0732   for (const auto& item : it->second) {
0733     auto& details = item.second;
0734     for (int ibin = 0; ibin < numberOfSamplesHF_; ++ibin) {
0735       const int idx = ibin + shift;
0736 
0737       int long_fiber_val = 0;
0738       int long_fiber_count = 0;
0739       int short_fiber_val = 0;
0740       int short_fiber_count = 0;
0741 
0742       bool saturated = false;
0743 
0744       for (auto i : {0, 2}) {
0745         if (idx < details[i].samples.size() and details[i].validity[idx] and details[i].passTDC[idx]) {
0746           long_fiber_val += details[i].samples[idx];
0747           saturated = saturated || (details[i].samples[idx] == QIE10_LINEARIZATION_ET);
0748           ++long_fiber_count;
0749         }
0750       }
0751       for (auto i : {1, 3}) {
0752         if (idx < details[i].samples.size() and details[i].validity[idx] and details[i].passTDC[idx]) {
0753           short_fiber_val += details[i].samples[idx];
0754           saturated = saturated || (details[i].samples[idx] == QIE10_LINEARIZATION_ET);
0755           ++short_fiber_count;
0756         }
0757       }
0758 
0759       if (saturated) {
0760         output[ibin] = QIE10_MAX_LINEARIZATION_ET;
0761       } else {
0762         // For details of the energy handling, see:
0763         // https://cms-docdb.cern.ch/cgi-bin/DocDB/ShowDocument?docid=12306
0764         // If both readouts are valid, average of the two energies is taken
0765         // division by 2 is compensated by adjusting the total scale shift in the end
0766         if (long_fiber_count == 2)
0767           long_fiber_val >>= 1;
0768         if (short_fiber_count == 2)
0769           short_fiber_val >>= 1;
0770 
0771         auto sum = long_fiber_val + short_fiber_val;
0772         // Similar to above, if both channels are valid,
0773         // average of the two energies is calculated
0774         // division by 2 here is also compensated by adjusting the total scale shift in the end
0775         if (long_fiber_count > 0 and short_fiber_count > 0)
0776           sum >>= 1;
0777 
0778         output[ibin] += sum;
0779       }
0780 
0781       for (const auto& detail : details) {
0782         if (idx < int(detail.digi.size()) and detail.validity[idx] and
0783             HcalDetId(detail.digi.id()).ietaAbs() >= FIRST_FINEGRAIN_TOWER) {
0784           if (useTDCInMinBiasBits_ && !detail.passTDC[idx])
0785             continue;
0786           finegrain[ibin][1] = finegrain[ibin][1] or detail.fgbits[idx][0];
0787           // what is commonly called the "second" HF min-bias bit is
0788           // actually the 0-th bit, which can also be used instead for the EM bit
0789           // (called finegrain[ibin][0] below) in non-HI running
0790           finegrain[ibin][0] = finegrain[ibin][0] or detail.fgbits[idx][1];
0791         }
0792       }
0793       // the EM bit is only used if the "second" FG bit is disabled
0794       if (embit != nullptr and FG_HF_thresholds_.at(1) != 255) {
0795         finegrain[ibin][0] = embit->fineGrainbit(details[1].digi,
0796                                                  details[3].digi,
0797                                                  details[0].digi,
0798                                                  details[2].digi,
0799                                                  details[1].validity[idx],
0800                                                  details[3].validity[idx],
0801                                                  details[0].validity[idx],
0802                                                  details[2].validity[idx],
0803                                                  idx);
0804       }
0805     }
0806   }
0807 
0808   for (int bin = 0; bin < numberOfSamplesHF_; ++bin) {
0809     output[bin] = min({(unsigned int)QIE10_MAX_LINEARIZATION_ET, output[bin] >> (hf_lumi_shift - 2)});
0810   }
0811   std::vector<int> finegrain_converted;
0812   finegrain_converted.reserve(finegrain.size());
0813   for (const auto& fg : finegrain)
0814     finegrain_converted.push_back(fg.to_ulong());
0815   outcoder_->compress(output, finegrain_converted, result);
0816 }
0817 
0818 void HcalTriggerPrimitiveAlgo::runZS(HcalTrigPrimDigiCollection& result) {
0819   for (HcalTrigPrimDigiCollection::iterator tp = result.begin(); tp != result.end(); ++tp) {
0820     bool ZS = true;
0821     for (int i = 0; i < tp->size(); ++i) {
0822       if (tp->sample(i).compressedEt() > ZS_threshold_I_) {
0823         ZS = false;
0824         break;
0825       }
0826     }
0827     if (ZS)
0828       tp->setZSInfo(false, true);
0829     else
0830       tp->setZSInfo(true, false);
0831   }
0832 }
0833 
0834 void HcalTriggerPrimitiveAlgo::runFEFormatError(const FEDRawDataCollection* rawraw,
0835                                                 const HcalElectronicsMap* emap,
0836                                                 HcalTrigPrimDigiCollection& result) {
0837   std::set<uint32_t> FrontEndErrors;
0838 
0839   for (int i = FEDNumbering::MINHCALFEDID; i <= FEDNumbering::MAXHCALFEDID; ++i) {
0840     const FEDRawData& raw = rawraw->FEDData(i);
0841     if (raw.size() < 12)
0842       continue;
0843     const HcalDCCHeader* dccHeader = (const HcalDCCHeader*)(raw.data());
0844     if (!dccHeader)
0845       continue;
0846     HcalHTRData htr;
0847     for (int spigot = 0; spigot < HcalDCCHeader::SPIGOT_COUNT; spigot++) {
0848       if (!dccHeader->getSpigotPresent(spigot))
0849         continue;
0850       dccHeader->getSpigotData(spigot, htr, raw.size());
0851       int dccid = dccHeader->getSourceId();
0852       int errWord = htr.getErrorsWord() & 0x1FFFF;
0853       bool HTRError = (!htr.check() || htr.isHistogramEvent() || (errWord & 0x800) != 0);
0854 
0855       if (HTRError) {
0856         bool valid = false;
0857         for (int fchan = 0; fchan < 3 && !valid; fchan++) {
0858           for (int fib = 0; fib < 9 && !valid; fib++) {
0859             HcalElectronicsId eid(fchan, fib, spigot, dccid - FEDNumbering::MINHCALFEDID);
0860             eid.setHTR(htr.readoutVMECrateId(), htr.htrSlot(), htr.htrTopBottom());
0861             DetId detId = emap->lookup(eid);
0862             if (detId.null())
0863               continue;
0864             HcalSubdetector subdet = (HcalSubdetector(detId.subdetId()));
0865             if (detId.det() != 4 || (subdet != HcalBarrel && subdet != HcalEndcap && subdet != HcalForward))
0866               continue;
0867             std::vector<HcalTrigTowerDetId> ids = theTrigTowerGeometry->towerIds(detId);
0868             for (std::vector<HcalTrigTowerDetId>::const_iterator triggerId = ids.begin(); triggerId != ids.end();
0869                  ++triggerId) {
0870               FrontEndErrors.insert(triggerId->rawId());
0871             }
0872             //valid = true;
0873           }
0874         }
0875       }
0876     }
0877   }
0878 
0879   // Loop over TP collection
0880   // Set TP to zero if there is FE Format Error
0881   HcalTriggerPrimitiveSample zeroSample(0);
0882   for (HcalTrigPrimDigiCollection::iterator tp = result.begin(); tp != result.end(); ++tp) {
0883     if (FrontEndErrors.find(tp->id().rawId()) != FrontEndErrors.end()) {
0884       for (int i = 0; i < tp->size(); ++i)
0885         tp->setSample(i, zeroSample);
0886     }
0887   }
0888 }
0889 
0890 void HcalTriggerPrimitiveAlgo::addFG(const HcalTrigTowerDetId& id, std::vector<bool>& msb) {
0891   FGbitMap::iterator itr = fgMap_.find(id);
0892   if (itr != fgMap_.end()) {
0893     std::vector<bool>& _msb = itr->second;
0894     for (size_t i = 0; i < msb.size(); ++i)
0895       _msb[i] = _msb[i] || msb[i];
0896   } else
0897     fgMap_[id] = msb;
0898 }
0899 
0900 bool HcalTriggerPrimitiveAlgo::validUpgradeFG(const HcalTrigTowerDetId& id, int depth) const {
0901   if (depth > LAST_FINEGRAIN_DEPTH)
0902     return false;
0903   if (id.ietaAbs() > LAST_FINEGRAIN_TOWER)
0904     return false;
0905   if (id.ietaAbs() == HBHE_OVERLAP_TOWER and not upgrade_hb_)
0906     return false;
0907   return true;
0908 }
0909 
0910 bool HcalTriggerPrimitiveAlgo::needLegacyFG(const HcalTrigTowerDetId& id) const {
0911   // This tower (ietaAbs == 16) does not accept upgraded FG bits,
0912   // but needs pseudo legacy ones to ensure that the tower is processed
0913   // even when the QIE8 depths in front of it do not have energy deposits.
0914   if (id.ietaAbs() == HBHE_OVERLAP_TOWER and not upgrade_hb_)
0915     return true;
0916   return false;
0917 }
0918 
0919 bool HcalTriggerPrimitiveAlgo::needUpgradeID(const HcalTrigTowerDetId& id, int depth) const {
0920   // Depth 7 for TT 26, 27, and 28 is not considered a fine grain depth.
0921   // However, the trigger tower for these ieta should still be added to the fgUpgradeMap_
0922   // Otherwise, depth 7-only signal will not be analyzed.
0923   unsigned int aieta = id.ietaAbs();
0924   if (aieta >= FIRST_DEPTH7_TOWER and aieta <= LAST_FINEGRAIN_TOWER and depth > LAST_FINEGRAIN_DEPTH)
0925     return true;
0926   return false;
0927 }
0928 
0929 void HcalTriggerPrimitiveAlgo::addUpgradeFG(const HcalTrigTowerDetId& id,
0930                                             int depth,
0931                                             const std::vector<std::bitset<2>>& bits) {
0932   if (not validUpgradeFG(id, depth)) {
0933     if (needLegacyFG(id)) {
0934       std::vector<bool> pseudo(bits.size(), false);
0935       addFG(id, pseudo);
0936     } else if (needUpgradeID(id, depth)) {
0937       // If the tower id is not in the map yet
0938       // then for safety's sake add it, otherwise, no need
0939       // Likewise, we're here with non-fg depth 7 so the bits are not to be added
0940       auto it = fgUpgradeMap_.find(id);
0941       if (it == fgUpgradeMap_.end()) {
0942         FGUpgradeContainer element;
0943         element.resize(bits.size());
0944         fgUpgradeMap_.insert(std::make_pair(id, element));
0945       }
0946     }
0947 
0948     return;
0949   }
0950 
0951   auto it = fgUpgradeMap_.find(id);
0952   if (it == fgUpgradeMap_.end()) {
0953     FGUpgradeContainer element;
0954     element.resize(bits.size());
0955     it = fgUpgradeMap_.insert(std::make_pair(id, element)).first;
0956   }
0957   for (unsigned int i = 0; i < bits.size(); ++i) {
0958     it->second[i][0][depth - 1] = bits[i][0];
0959     it->second[i][1][depth - 1] = bits[i][1];
0960   }
0961 }
0962 
0963 void HcalTriggerPrimitiveAlgo::addUpgradeTDCFG(const HcalTrigTowerDetId& id, const QIE11DataFrame& frame) {
0964   HcalDetId detId(frame.id());
0965   if (detId.subdet() != HcalEndcap && detId.subdet() != HcalBarrel)
0966     return;
0967 
0968   std::vector<HcalTrigTowerDetId> ids = theTrigTowerGeometry->towerIds(detId);
0969   assert(ids.size() == 1 || ids.size() == 2);
0970   IntegerCaloSamples samples1(ids[0], int(frame.samples()));
0971   samples1.setPresamples(frame.presamples());
0972   incoder_->adc2Linear(frame, samples1);                                  // use linearization LUT
0973   std::vector<unsigned short> bits12_15 = incoder_->group0FGbits(frame);  // get 4 energy bits (12-15) from group 0 LUT
0974 
0975   bool is_compressed = false;
0976   if (detId.subdet() == HcalBarrel) {
0977     is_compressed = (frame.flavor() == 3);
0978     // 0 if frame.flavor is 0 (uncompressed), 1 if frame.flavor is 3 (compressed)
0979   }
0980 
0981   auto it = fgUpgradeTDCMap_.find(id);
0982   if (it == fgUpgradeTDCMap_.end()) {
0983     FGUpgradeTDCContainer element;
0984     element.resize(frame.samples());
0985     it = fgUpgradeTDCMap_.insert(std::make_pair(id, element)).first;
0986   }
0987   for (int i = 0; i < frame.samples(); i++) {
0988     it->second[i][detId.depth() - 1] =
0989         std::make_pair(std::make_pair(bits12_15[i], is_compressed), std::make_pair(frame[i].tdc(), samples1[i]));
0990   }
0991 }
0992 
0993 void HcalTriggerPrimitiveAlgo::setWeightsQIE11(const edm::ParameterSet& weightsQIE11) {
0994   // Names are just abs(ieta) for HBHE
0995   std::vector<std::string> ietaStrs = weightsQIE11.getParameterNames();
0996   for (auto& ietaStr : ietaStrs) {
0997     // Strip off "ieta" part of key and just use integer value in map
0998     auto const& v = weightsQIE11.getParameter<std::vector<int>>(ietaStr);
0999     weightsQIE11_[std::stoi(ietaStr.substr(4))] = {{v[0], v[1]}};
1000   }
1001 }
1002 
1003 void HcalTriggerPrimitiveAlgo::setWeightQIE11(int aieta, int weight) {
1004   // Simple map of |ieta| in HBHE to weight
1005   // Only one weight for SOI-1 TS
1006   weightsQIE11_[aieta] = {{weight, 255}};
1007 }
1008 
1009 void HcalTriggerPrimitiveAlgo::setPeakFinderAlgorithm(int algo) {
1010   if (algo <= 0 || algo > 2)
1011     throw cms::Exception("ERROR: Only algo 1 & 2 are supported.") << std::endl;
1012   peak_finder_algorithm_ = algo;
1013 }
1014 
1015 void HcalTriggerPrimitiveAlgo::setNCTScaleShift(int shift) { NCTScaleShift = shift; }
1016 
1017 void HcalTriggerPrimitiveAlgo::setRCTScaleShift(int shift) { RCTScaleShift = shift; }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!