EcalSelectiveReadoutAlgos/src/EcalSelectiveReadoutSuppressor.cc

0001 #include "SimCalorimetry/EcalSelectiveReadoutAlgos/interface/EcalSelectiveReadoutSuppressor.h"
0002 #include "SimCalorimetry/EcalSelectiveReadoutAlgos/interface/EcalSelectiveReadout.h"
0003 #include "DataFormats/EcalDigi/interface/EEDataFrame.h"
0004 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 // Geometry
0008 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0009 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0010 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0011 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0012 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0013
0014 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0015
0016 //exceptions:
0017 #include "FWCore/Utilities/interface/Exception.h"
0018
0019 #include <cmath>
0020 #include <iostream>
0021 #include <limits>
0022 #include <memory>
0023
0024 using namespace boost;
0025 using namespace std;
0026
0027 const int EcalSelectiveReadoutSuppressor::nFIRTaps = 6;
0028
0029 EcalSelectiveReadoutSuppressor::EcalSelectiveReadoutSuppressor(const edm::ParameterSet& params,
0030                                                                edm::ConsumesCollector iC)
0031     : ttThresOnCompressedEt_(false), ievt_(0), geoToken_(iC.esConsumes()) {
0032   int defTtf = params.getParameter<int>("defaultTtf");
0033   if (defTtf < 0 || defTtf > 7) {
0034     throw cms::Exception("InvalidParameter") << "Value of EcalSelectiveReadoutProducer module parameter defaultTtf, "
0035                                              << defaultTtf_ << ", is out of the valid range 0..7\n";
0036   } else {
0037     defaultTtf_ = (EcalSelectiveReadout::ttFlag_t)defTtf;
0038   }
0039
0040   trigPrimBypass_ = params.getParameter<bool>("trigPrimBypass");
0041   trigPrimBypassMode_ = params.getParameter<int>("trigPrimBypassMode");
0042   trigPrimBypassWithPeakFinder_ = params.getParameter<bool>("trigPrimBypassWithPeakFinder");
0043   trigPrimBypassLTH_ = params.getParameter<double>("trigPrimBypassLTH");
0044   trigPrimBypassHTH_ = params.getParameter<double>("trigPrimBypassHTH");
0045   if (trigPrimBypass_) {
0046     edm::LogWarning("Digitization") << "Beware a simplified trigger primitive "
0047                                        "computation is used for the ECAL selective readout";
0048     if (trigPrimBypassMode_ != 0 && trigPrimBypassMode_ != 1) {
0049       throw cms::Exception("InvalidParameter")
0050           << "Invalid value for EcalSelectiveReadoutProducer parameter 'trigPrimBypassMode_'."
0051              " Valid values are 0 and 1.\n";
0052     }
0053     ttThresOnCompressedEt_ = (trigPrimBypassMode_ == 1);
0054   }
0055 }
0056
0057 void EcalSelectiveReadoutSuppressor::setSettings(const EcalSRSettings* settings) {
0058   firstFIRSample = settings->ecalDccZs1stSample_[0];
0059   weights = settings->dccNormalizedWeights_[0];
0060   symetricZS = settings->symetricZS_[0];
0061   actions_ = settings->actions_;
0062
0063   //online configuration has only 4 actions flags, the 4 'forced' flags being the same with the force
0064   //bit set to 1. Extends the actions vector for case of online-type configuration:
0065   if (actions_.size() == 4) {
0066     for (int i = 0; i < 4; ++i) {
0067       actions_.push_back(actions_[i] | 0x4);
0068     }
0069   }
0070
0071   bool actionValid = actions_.size() == 8;
0072   for (size_t i = 0; i < actions_.size(); ++i) {
0073     if (actions_[i] < 0 || actions_[i] > 7)
0074       actionValid = false;
0075   }
0076
0077   if (!actionValid) {
0078     throw cms::Exception("InvalidParameter")
0079         << "EcalSelectiveReadoutProducer module parameter 'actions' is "
0080            "not valid. It must be a vector of 8 integer values comprised between 0 and 7\n";
0081   }
0082
0083   double adcToGeV = settings->ebDccAdcToGeV_;
0084   thrUnit[BARREL] = adcToGeV / 4.;  //unit=1/4th ADC count
0085
0086   adcToGeV = settings->eeDccAdcToGeV_;
0087   thrUnit[ENDCAP] = adcToGeV / 4.;  //unit=1/4th ADC count
0088   ecalSelectiveReadout = std::make_unique<EcalSelectiveReadout>(settings->deltaEta_[0], settings->deltaPhi_[0]);
0089   const int eb = 0;
0090   const int ee = 1;
0091   initCellThresholds(settings->srpLowInterestChannelZS_[eb],
0092                      settings->srpLowInterestChannelZS_[ee],
0093                      settings->srpHighInterestChannelZS_[eb],
0094                      settings->srpHighInterestChannelZS_[ee]);
0095 }
0096
0097 void EcalSelectiveReadoutSuppressor::setTriggerMap(const EcalTrigTowerConstituentsMap* map) {
0098   theTriggerMap = map;
0099   ecalSelectiveReadout->setTriggerMap(map);
0100 }
0101
0102 void EcalSelectiveReadoutSuppressor::setElecMap(const EcalElectronicsMapping* map) {
0103   ecalSelectiveReadout->setElecMap(map);
0104 }
0105
0106 void EcalSelectiveReadoutSuppressor::setGeometry(const CaloGeometry* caloGeometry) {
0107 #ifndef ECALSELECTIVEREADOUT_NOGEOM
0108   ecalSelectiveReadout->setGeometry(caloGeometry);
0109 #endif
0110 }
0111
0112 void EcalSelectiveReadoutSuppressor::initCellThresholds(double barrelLowInterest,
0113                                                         double endcapLowInterest,
0114                                                         double barrelHighInterest,
0115                                                         double endcapHighInterest) {
0116   //center, neighbour and single RUs are grouped into a single
0117   //'high interest' group
0118   int lowInterestThr[2];  //index for BARREL/ENDCAP
0119   //  int lowInterestSrFlag[2];
0120   int highInterestThr[2];
0121   // int highInterestSrFlag[2];
0122
0123   lowInterestThr[BARREL] = internalThreshold(barrelLowInterest, BARREL);
0124   //  lowInterestSrFlag[BARREL] = thr2Srf(lowInterestThr[BARREL],
0125   //                      EcalSrFlag::SRF_ZS1);
0126
0127   highInterestThr[BARREL] = internalThreshold(barrelHighInterest, BARREL);
0128   //  highInterestSrFlag[BARREL] = thr2Srf(highInterestThr[BARREL],
0129   //                       EcalSrFlag::SRF_ZS2);
0130
0131   lowInterestThr[ENDCAP] = internalThreshold(endcapLowInterest, ENDCAP);
0132   //lowInterestSrFlag[ENDCAP] = thr2Srf(lowInterestThr[ENDCAP],
0133   //                  EcalSrFlag::SRF_ZS1);
0134
0135   highInterestThr[ENDCAP] = internalThreshold(endcapHighInterest, ENDCAP);
0136   //  highInterestSrFlag[ENDCAP] = thr2Srf(highInterestThr[ENDCAP],
0137   //                   EcalSrFlag::SRF_ZS2);
0138
0139   const int FORCED_MASK = EcalSelectiveReadout::FORCED_MASK;
0140
0141   for (int iSubDet = 0; iSubDet < 2; ++iSubDet) {
0142     //low interest
0143     //zsThreshold[iSubDet][0] = lowInterestThr[iSubDet];
0144     //srFlags[iSubDet][0] = lowInterestSrFlag[iSubDet];
0145     //srFlags[iSubDet][0 + FORCED_MASK] = FORCED_MASK | lowInterestSrFlag[iSubDet];
0146
0147     //single->high interest
0148     //zsThreshold[iSubDet][1] = highInterestThr[iSubDet];
0149     //srFlags[iSubDet][1] = highInterestSrFlag[iSubDet];
0150     //srFlags[iSubDet][1 +  FORCED_MASK] = FORCED_MASK | highInterestSrFlag[iSubDet];
0151
0152     //neighbour->high interest
0153     //zsThreshold[iSubDet][2] = highInterestThr[iSubDet];
0154     //srFlags[iSubDet][2] = highInterestSrFlag[iSubDet];
0155     //srFlags[iSubDet][2 + FORCED_MASK] = FORCED_MASK | highInterestSrFlag[iSubDet];
0156
0157     //center->high interest
0158     //zsThreshold[iSubDet][3] = highInterestThr[iSubDet];
0159     //srFlags[iSubDet][3] = highInterestSrFlag[iSubDet];
0160     //srFlags[iSubDet][3 + FORCED_MASK] = FORCED_MASK | highInterestSrFlag[iSubDet];
0161     for (size_t i = 0; i < 8; ++i) {
0162       srFlags[iSubDet][i] = actions_[i];
0163       if ((actions_[i] & ~FORCED_MASK) == 0)
0164         zsThreshold[iSubDet][i] = numeric_limits<int>::max();
0165       else if ((actions_[i] & ~FORCED_MASK) == 1)
0166         zsThreshold[iSubDet][i] = lowInterestThr[iSubDet];
0167       else if ((actions_[i] & ~FORCED_MASK) == 2)
0168         zsThreshold[iSubDet][i] = highInterestThr[iSubDet];
0169       else
0170         zsThreshold[iSubDet][i] = numeric_limits<int>::min();
0171     }
0172
0173     //     for(size_t i = 0; i < 8; ++i){
0174     //       cout << "zsThreshold[" << iSubDet << "]["  << i << "] = " << zsThreshold[iSubDet][i] << endl;
0175     //     }
0176   }
0177 }
0178
0179 // int EcalSelectiveReadoutSuppressor::thr2Srf(int thr, int zsFlag) const{
0180 //   if(thr==numeric_limits<int>::max()){
0181 //     return EcalSrFlag::SRF_SUPPRESS;
0182 //   }
0183 //   if(thr==numeric_limits<int>::min()){
0184 //     return EcalSrFlag::SRF_FULL;
0185 //   }
0186 //   return zsFlag;
0187 // }
0188
0189 int EcalSelectiveReadoutSuppressor::internalThreshold(double thresholdInGeV, int iSubDet) const {
0190   double thr_ = thresholdInGeV / thrUnit[iSubDet];
0191   //treating over- and underflows, threshold is coded on 11+1 signed bits
0192   //an underflow threshold is considered here as if NoRO DCC switch is on
0193   //an overflow threshold is considered here as if ForcedRO DCC switch in on
0194   //Beware that conparison must be done on a double type, because conversion
0195   //cast to an int of a double higher than MAX_INT is undefined.
0196   int thr;
0197   if (thr_ >= 0x7FF + .5) {
0198     thr = numeric_limits<int>::max();
0199   } else if (thr_ <= -0x7FF - .5) {
0200     thr = numeric_limits<int>::min();
0201   } else {
0202     thr = lround(thr_);
0203   }
0204   return thr;
0205 }
0206
0207 //This implementation  assumes that int is coded on at least 28-bits,
0208 //which in pratice should be always true.
0209 bool EcalSelectiveReadoutSuppressor::accept(edm::DataFrame const& frame, int thr) {
0210   //FIR filter weights:
0211   const vector<int>& w = getFIRWeigths();
0212
0213   //accumulator used to compute weighted sum of samples
0214   int acc = 0;
0215   bool gain12saturated = false;
0216   const int gain12 = 0x01;
0217   const int lastFIRSample = firstFIRSample + nFIRTaps - 1;
0218   //LogDebug("DccFir") << "DCC FIR operation: ";
0219   int iWeight = 0;
0220   for (int iSample = firstFIRSample - 1; iSample < lastFIRSample; ++iSample, ++iWeight) {
0221     if (iSample >= 0 && iSample < (int)frame.size()) {
0222       EcalMGPASample sample(frame[iSample]);
0223       if (sample.gainId() != gain12)
0224         gain12saturated = true;
0225       //LogTrace("DccFir") << (iSample>=firstFIRSample?"+":"") << sample.adc()
0226       //         << "*(" << w[iWeight] << ")";
0227       acc += sample.adc() * w[iWeight];
0228     } else {
0229       edm::LogWarning("DccFir") << __FILE__ << ":" << __LINE__
0230                                 << ": Not enough samples in data frame or 'ecalDccZs1stSample' module "
0231                                    "parameter is not valid...";
0232     }
0233   }
0234
0235   if (symetricZS) {  //cut on absolute value
0236     if (acc < 0)
0237       acc = -acc;
0238   }
0239
0240   //LogTrace("DccFir") << "\n";
0241   //discards the 8 LSBs
0242   //(result of shift operator on negative numbers depends on compiler
0243   //implementation, therefore the value is offset to make sure it is positive
0244   //before performing the bit shift).
0245   acc = ((acc + (1 << 30)) >> 8) - (1 << (30 - 8));
0246
0247   //ZS passed if weigthed sum acc above ZS threshold or if
0248   //one sample has a lower gain than gain 12 (that is gain 12 output
0249   //is saturated)
0250
0251   const bool result = (acc >= thr) || gain12saturated;
0252
0253   //LogTrace("DccFir") << "acc: " << acc << "\n"
0254   //                   << "threshold: " << thr << " ("
0255   //                   << thr*thrUnit[((EcalDataFrame&)frame).id().subdetId()==EcalBarrel?0:1]
0256   //                   << "GeV)\n"
0257   //                   << "saturated: " << (gain12saturated?"yes":"no") << "\n"
0258   //                   << "ZS passed: " << (result?"yes":"no")
0259   //                   << (symetricZS?" (symetric cut)":"") << "\n";
0260
0261   return result;
0262 }
0263
0264 void EcalSelectiveReadoutSuppressor::run(const edm::EventSetup& eventSetup,
0265                                          const EcalTrigPrimDigiCollection& trigPrims,
0266                                          EBDigiCollection& barrelDigis,
0267                                          EEDigiCollection& endcapDigis) {
0268   EBDigiCollection selectedBarrelDigis;
0269   EEDigiCollection selectedEndcapDigis;
0270
0271   run(eventSetup, trigPrims, barrelDigis, endcapDigis, &selectedBarrelDigis, &selectedEndcapDigis, nullptr, nullptr);
0272
0273   //replaces the input with the suppressed version
0274   barrelDigis.swap(selectedBarrelDigis);
0275   endcapDigis.swap(selectedEndcapDigis);
0276 }
0277
0278 void EcalSelectiveReadoutSuppressor::run(const edm::EventSetup& eventSetup,
0279                                          const EcalTrigPrimDigiCollection& trigPrims,
0280                                          const EBDigiCollection& barrelDigis,
0281                                          const EEDigiCollection& endcapDigis,
0282                                          EBDigiCollection* selectedBarrelDigis,
0283                                          EEDigiCollection* selectedEndcapDigis,
0284                                          EBSrFlagCollection* ebSrFlags,
0285                                          EESrFlagCollection* eeSrFlags) {
0286   ++ievt_;
0287   if (!trigPrimBypass_ || ttThresOnCompressedEt_) {  //uses output of TPG emulator
0288     setTtFlags(trigPrims);
0289   } else {  //debug mode, run w/o TP digis
0290     setTtFlags(eventSetup, barrelDigis, endcapDigis);
0291   }
0292
0293   ecalSelectiveReadout->runSelectiveReadout0(ttFlags);
0294
0295   if (selectedBarrelDigis) {
0296     selectedBarrelDigis->reserve(barrelDigis.size() / 20);
0297
0298     // do barrel first
0299     for (EBDigiCollection::const_iterator digiItr = barrelDigis.begin(); digiItr != barrelDigis.end(); ++digiItr) {
0300       int interestLevel = ecalSelectiveReadout->getCrystalInterest(EBDigiCollection::DetId(digiItr->id())) &
0301                           ~EcalSelectiveReadout::FORCED_MASK;
0302       if (accept(*digiItr, zsThreshold[BARREL][interestLevel])) {
0303         selectedBarrelDigis->push_back(digiItr->id(), digiItr->begin());
0304       }
0305     }
0306   }
0307
0308   // and endcaps
0309   if (selectedEndcapDigis) {
0310     selectedEndcapDigis->reserve(endcapDigis.size() / 20);
0311     for (EEDigiCollection::const_iterator digiItr = endcapDigis.begin(); digiItr != endcapDigis.end(); ++digiItr) {
0312       int interestLevel = ecalSelectiveReadout->getCrystalInterest(EEDigiCollection::DetId(digiItr->id())) &
0313                           ~EcalSelectiveReadout::FORCED_MASK;
0314       if (accept(*digiItr, zsThreshold[ENDCAP][interestLevel])) {
0315         selectedEndcapDigis->push_back(digiItr->id(), digiItr->begin());
0316       }
0317     }
0318   }
0319
0320   if (ievt_ <= 10) {
0321     int neb = (selectedBarrelDigis ? selectedBarrelDigis->size() : 0);
0322     if (selectedEndcapDigis)
0323       LogDebug("EcalSelectiveReadout")
0324           //       << __FILE__ << ":" << __LINE__ << ": "
0325           << "Number of EB digis passing the SR: " << neb << " / " << barrelDigis.size() << "\n";
0326     if (selectedEndcapDigis)
0327       LogDebug("EcalSelectiveReadout")
0328           //       << __FILE__ << ":" << __LINE__ << ": "
0329           << "\nNumber of EE digis passing the SR: " << selectedEndcapDigis->size() << " / " << endcapDigis.size()
0330           << "\n";
0331   }
0332
0333   if (ebSrFlags)
0334     ebSrFlags->reserve(34 * 72);
0335   if (eeSrFlags)
0336     eeSrFlags->reserve(624);
0337   //SR flags:
0338   for (int iZ = -1; iZ <= 1; iZ += 2) {  //-1=>EE-, EB-, +1=>EE+, EB+
0339     //barrel:
0340     for (unsigned iEta = 1; iEta <= nBarrelTriggerTowersInEta / 2; ++iEta) {
0341       for (unsigned iPhi = 1; iPhi <= nTriggerTowersInPhi; ++iPhi) {
0342         const EcalTrigTowerDetId id(iZ, EcalBarrel, iEta, iPhi);
0343         EcalSelectiveReadout::towerInterest_t interest = ecalSelectiveReadout->getTowerInterest(id);
0344         if (interest < 0) {
0345           throw cms::Exception("EcalSelectiveReadout") << __FILE__ << ":" << __LINE__ << ": "
0346                                                        << "unknown SR flag. for "
0347                                                        << " TT " << id << ". Most probably a bug.";
0348         }
0349         int flag;
0350         //  if(interest==EcalSelectiveReadout::FORCED_RO){
0351         //  flag = EcalSrFlag::SRF_FORCED_MASK | EcalSrFlag::SRF_FULL;
0352         //} else{
0353         flag = srFlags[BARREL][interest];
0354         //}
0355         if (ebSrFlags)
0356           ebSrFlags->push_back(EBSrFlag(id, flag));
0357       }  //next iPhi
0358     }    //next barrel iEta
0359
0360     //endcap:
0361     EcalScDetId id;
0362     for (int iX = 1; iX <= 20; ++iX) {
0363       for (int iY = 1; iY <= 20; ++iY) {
0364         if (EcalScDetId::validDetId(iX, iY, iZ))
0365           id = EcalScDetId(iX, iY, iZ);
0366         else
0367           continue;
0368
0369         EcalSelectiveReadout::towerInterest_t interest = ecalSelectiveReadout->getSuperCrystalInterest(id);
0370
0371         if (interest >= 0) {  //negative no SC at (iX,iY) coordinates
0372           int flag;
0373           //      if(interest==EcalSelectiveReadout::FORCED_RO){
0374           //  flag = EcalSrFlag::SRF_FORCED_MASK | EcalSrFlag::SRF_FULL;
0375           //} else{
0376           flag = srFlags[ENDCAP][interest];
0377           //}
0378           if (eeSrFlags)
0379             eeSrFlags->push_back(EESrFlag(id, flag));
0380         } else if (iX < 9 || iX > 12 || iY < 9 || iY > 12) {  //not an inner partial SC
0381           edm::LogError("EcalSelectiveReadout") << __FILE__ << ":" << __LINE__ << ": "
0382                                                 << "negative interest in EE for SC " << id << "\n";
0383         }
0384       }  //next iY
0385     }    //next iX
0386   }
0387 }
0388
0389 void EcalSelectiveReadoutSuppressor::setTtFlags(const EcalTrigPrimDigiCollection& trigPrims) {
0390   for (size_t iEta0 = 0; iEta0 < nTriggerTowersInEta; ++iEta0) {
0391     for (size_t iPhi0 = 0; iPhi0 < nTriggerTowersInPhi; ++iPhi0) {
0392       ttFlags[iEta0][iPhi0] = defaultTtf_;
0393     }
0394   }
0395   for (EcalTrigPrimDigiCollection::const_iterator trigPrim = trigPrims.begin(); trigPrim != trigPrims.end();
0396        ++trigPrim) {
0397     int iEta = trigPrim->id().ieta();
0398     unsigned int iEta0;
0399     if (iEta < 0) {  //z- half ECAL: transforming ranges -28;-1 => 0;27
0400       iEta0 = iEta + nTriggerTowersInEta / 2;
0401     } else {  //z+ halfECAL: transforming ranges 1;28 => 28;55
0402       iEta0 = iEta + nTriggerTowersInEta / 2 - 1;
0403     }
0404
0405     unsigned int iPhi0 = trigPrim->id().iphi() - 1;
0406
0407     if (!ttThresOnCompressedEt_) {
0408       ttFlags[iEta0][iPhi0] = (EcalSelectiveReadout::ttFlag_t)trigPrim->ttFlag();
0409     } else {
0410       int compressedEt = trigPrim->compressedEt();
0411       if (compressedEt < trigPrimBypassLTH_) {
0412         ttFlags[iEta0][iPhi0] = EcalSelectiveReadout::TTF_LOW_INTEREST;
0413       } else if (compressedEt < trigPrimBypassHTH_) {
0414         ttFlags[iEta0][iPhi0] = EcalSelectiveReadout::TTF_MID_INTEREST;
0415       } else {
0416         ttFlags[iEta0][iPhi0] = EcalSelectiveReadout::TTF_HIGH_INTEREST;
0417       }
0418     }
0419   }
0420 }
0421
0422 vector<int> EcalSelectiveReadoutSuppressor::getFIRWeigths() {
0423   if (firWeights.empty()) {
0424     firWeights = vector<int>(nFIRTaps, 0);  //default weight: 0;
0425     const static int maxWeight = 0xEFF;     //weights coded on 11+1 signed bits
0426     for (unsigned i = 0; i < min((size_t)nFIRTaps, weights.size()); ++i) {
0427       firWeights[i] = lround(weights[i] * (1 << 10));
0428       if (abs(firWeights[i]) > maxWeight) {  //overflow
0429         firWeights[i] = firWeights[i] < 0 ? -maxWeight : maxWeight;
0430       }
0431     }
0432   }
0433   return firWeights;
0434 }
0435
0436 void EcalSelectiveReadoutSuppressor::setTtFlags(const edm::EventSetup& es,
0437                                                 const EBDigiCollection& ebDigis,
0438                                                 const EEDigiCollection& eeDigis) {
0439   double trigPrim[nTriggerTowersInEta][nTriggerTowersInPhi];
0440
0441   //ecal geometry:
0442   //  static const CaloSubdetectorGeometry* eeGeometry = 0;
0443   //  static const CaloSubdetectorGeometry* ebGeometry = 0;
0444   const CaloSubdetectorGeometry* eeGeometry = nullptr;
0445   const CaloSubdetectorGeometry* ebGeometry = nullptr;
0446   //  if(eeGeometry==0 || ebGeometry==0){
0447   auto const& geo = es.getData(geoToken_);
0448   eeGeometry = geo.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0449   ebGeometry = geo.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0450   //  }
0451
0452   //init trigPrim array:
0453   bzero(trigPrim, sizeof(trigPrim));
0454
0455   for (EBDigiCollection::const_iterator it = ebDigis.begin(); it != ebDigis.end(); ++it) {
0456     EBDataFrame frame(*it);
0457     const EcalTrigTowerDetId& ttId = theTriggerMap->towerOf(frame.id());
0458     //      edm:::LogDebug("TT") << __FILE__ << ":" << __LINE__ << ": "
0459     //       <<  ((EBDetId&)frame.id()).ieta()
0460     //       << "," << ((EBDetId&)frame.id()).iphi()
0461     //       << " -> " << ttId.ieta() << "," << ttId.iphi() << "\n";
0462     const int iTTEta0 = iTTEta2cIndex(ttId.ieta());
0463     const int iTTPhi0 = iTTPhi2cIndex(ttId.iphi());
0464     double theta = ebGeometry->getGeometry(frame.id())->getPosition().theta();
0465     double e = frame2Energy(frame);
0466     if (!trigPrimBypassWithPeakFinder_ || ((frame2Energy(frame, -1) < e) && (frame2Energy(frame, 1) < e))) {
0467       trigPrim[iTTEta0][iTTPhi0] += e * sin(theta);
0468     }
0469   }
0470
0471   for (EEDigiCollection::const_iterator it = eeDigis.begin(); it != eeDigis.end(); ++it) {
0472     EEDataFrame frame(*it);
0473     const EcalTrigTowerDetId& ttId = theTriggerMap->towerOf(frame.id());
0474     const int iTTEta0 = iTTEta2cIndex(ttId.ieta());
0475     const int iTTPhi0 = iTTPhi2cIndex(ttId.iphi());
0476     //     cout << __FILE__ << ":" << __LINE__ << ": EE xtal->TT "
0477     //   <<  ((EEDetId&)frame.id()).ix()
0478     //   << "," << ((EEDetId&)frame.id()).iy()
0479     //   << " -> " << ttId.ieta() << "," << ttId.iphi() << "\n";
0480     double theta = eeGeometry->getGeometry(frame.id())->getPosition().theta();
0481     double e = frame2Energy(frame);
0482     if (!trigPrimBypassWithPeakFinder_ || ((frame2Energy(frame, -1) < e) && (frame2Energy(frame, 1) < e))) {
0483       trigPrim[iTTEta0][iTTPhi0] += e * sin(theta);
0484     }
0485   }
0486
0487   //dealing with pseudo-TT in two inner EE eta-ring:
0488   int innerTTEtas[] = {0, 1, 54, 55};
0489   for (unsigned iRing = 0; iRing < sizeof(innerTTEtas) / sizeof(innerTTEtas[0]); ++iRing) {
0490     int iTTEta0 = innerTTEtas[iRing];
0491     //this detector eta-section is divided in only 36 phi bins
0492     //For this eta regions,
0493     //current tower eta numbering scheme is inconsistent. For geometry
0494     //version 133:
0495     //- TT are numbered from 0 to 72 for 36 bins
0496     //- some TT have an even index, some an odd index
0497     //For geometry version 125, there are 72 phi bins.
0498     //The code below should handle both geometry definition.
0499     //If there are 72 input trigger primitives for each inner eta-ring,
0500     //then the average of the trigger primitive of the two pseudo-TT of
0501     //a pair (nEta, nEta+1) is taken as Et of both pseudo TTs.
0502     //If there are only 36 input TTs for each inner eta ring, then half
0503     //of the present primitive of a pseudo TT pair is used as Et of both
0504     //pseudo TTs.
0505
0506     for (unsigned iTTPhi0 = 0; iTTPhi0 < nTriggerTowersInPhi - 1; iTTPhi0 += 2) {
0507       double et = .5 * (trigPrim[iTTEta0][iTTPhi0] + trigPrim[iTTEta0][iTTPhi0 + 1]);
0508       //divides the TT into 2 phi bins in order to match with 72 phi-bins SRP
0509       //scheme or average the Et on the two pseudo TTs if the TT is already
0510       //divided into two trigger primitives.
0511       trigPrim[iTTEta0][iTTPhi0] = et;
0512       trigPrim[iTTEta0][iTTPhi0 + 1] = et;
0513     }
0514   }
0515
0516   for (unsigned iTTEta0 = 0; iTTEta0 < nTriggerTowersInEta; ++iTTEta0) {
0517     for (unsigned iTTPhi0 = 0; iTTPhi0 < nTriggerTowersInPhi; ++iTTPhi0) {
0518       if (trigPrim[iTTEta0][iTTPhi0] > trigPrimBypassHTH_) {
0519         ttFlags[iTTEta0][iTTPhi0] = EcalSelectiveReadout::TTF_HIGH_INTEREST;
0520       } else if (trigPrim[iTTEta0][iTTPhi0] > trigPrimBypassLTH_) {
0521         ttFlags[iTTEta0][iTTPhi0] = EcalSelectiveReadout::TTF_MID_INTEREST;
0522       } else {
0523         ttFlags[iTTEta0][iTTPhi0] = EcalSelectiveReadout::TTF_LOW_INTEREST;
0524       }
0525
0526       // cout /*LogDebug("TT")*/
0527       //    << "ttFlags[" << iTTEta0 << "][" << iTTPhi0 << "] = "
0528       //    << ttFlags[iTTEta0][iTTPhi0] << "\n";
0529     }
0530   }
0531 }
0532
0533 template <class T>
0534 double EcalSelectiveReadoutSuppressor::frame2Energy(const T& frame, int offset) const {
0535   //we have to start by 0 in order to handle offset=-1
0536   //(however Fenix FIR has AFAK only 5 taps)
0537   double weights[] = {0., -1 / 3., -1 / 3., -1 / 3., 0., 1.};
0538
0539   double adc2GeV = 0.;
0540   if (typeid(frame) == typeid(EBDataFrame)) {
0541     adc2GeV = 0.035;
0542   } else if (typeid(frame) == typeid(EEDataFrame)) {
0543     adc2GeV = 0.060;
0544   } else {  //T is an invalid type!
0545     //TODO: replace message by a cms exception
0546     cerr << "Severe error. " << __FILE__ << ":" << __LINE__ << ": "
0547          << "this is a bug. Please report it.\n";
0548   }
0549
0550   double acc = 0;
0551
0552   const int n = min<int>(frame.size(), sizeof(weights) / sizeof(weights[0]));
0553
0554   double gainInv[] = {12., 1., 6., 12.};
0555
0556   //cout << __PRETTY_FUNCTION__ << ": ";
0557   for (int i = offset; i < n; ++i) {
0558     int iframe = i + offset;
0559     if (iframe >= 0 && iframe < frame.size()) {
0560       acc += weights[i] * frame[iframe].adc() * gainInv[frame[iframe].gainId()] * adc2GeV;
0561       //cout << (iframe>offset?"+":"")
0562       //     << frame[iframe].adc() << "*" << gainInv[frame[iframe].gainId()]
0563       //     << "*" << adc2GeV << "*(" << weights[i] << ")";
0564     }
0565   }
0566   //cout << "\n";
0567   return acc;
0568 }
0569
0570 void EcalSelectiveReadoutSuppressor::printTTFlags(ostream& os, int iEvent, bool withHeader) const {
0571   const char tccFlagMarker[] = {'?', '.', 'S', '?', 'C', 'E', 'E', 'E', 'E'};
0572   const int nEta = EcalSelectiveReadout::nTriggerTowersInEta;
0573   const int nPhi = EcalSelectiveReadout::nTriggerTowersInPhi;
0574
0575   if (withHeader) {
0576     os << "# TCC flag map\n#\n"
0577           "# +-->Phi            "
0578        << tccFlagMarker[1 + 0]
0579        << ": 000 (low interest)\n"
0580           "# |                  "
0581        << tccFlagMarker[1 + 1]
0582        << ": 001 (mid interest)\n"
0583           "# |                  "
0584        << tccFlagMarker[1 + 2]
0585        << ": 010 (not valid)\n"
0586           "# V Eta              "
0587        << tccFlagMarker[1 + 3]
0588        << ": 011 (high interest)\n"
0589           "#                    "
0590        << tccFlagMarker[1 + 4] << ": 1xx forced readout (Hw error)\n";
0591   }
0592
0593   if (iEvent >= 0) {
0594     os << "#\n#Event " << iEvent << "\n";
0595   }
0596
0597   for (int iEta = 0; iEta < nEta; ++iEta) {
0598     for (int iPhi = 0; iPhi < nPhi; ++iPhi) {
0599       os << tccFlagMarker[ttFlags[iEta][iPhi] + 1];
0600     }
0601     os << "\n";
0602   }
0603 }