Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-10-10 23:56:21 

0001 /** \class EcalRecHitWorkerRecover
0002   *  Algorithms to recover dead channels
0003   *
0004   */
0005 
0006 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h"
0007 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
0008 #include "CalibCalorimetry/EcalTPGTools/interface/EcalTPGScale.h"
0009 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
0010 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0011 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
0012 #include "CondFormats/DataRecord/interface/EcalTimeCalibConstantsRcd.h"
0013 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0014 #include "CondFormats/EcalObjects/interface/EcalTimeCalibConstants.h"
0015 #include "DataFormats/EcalDetId/interface/EcalScDetId.h"
0016 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Utilities/interface/ESGetToken.h"
0022 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0023 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0024 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0025 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0026 #include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
0027 #include "Geometry/EcalMapping/interface/EcalMappingRcd.h"
0028 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0029 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0030 #include "Geometry/Records/interface/EcalBarrelGeometryRecord.h"
0031 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0032 #include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h"
0033 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalRecHitSimpleAlgo.h"
0034 #include "RecoLocalCalo/EcalRecProducers/interface/EcalRecHitWorkerBaseClass.h"
0035 
0036 #include <vector>
0037 
0038 class EcalRecHitWorkerRecover : public EcalRecHitWorkerBaseClass {
0039 public:
0040   EcalRecHitWorkerRecover(const edm::ParameterSet&, edm::ConsumesCollector& c);
0041 
0042   void set(const edm::EventSetup& es) override;
0043   bool run(const edm::Event& evt, const EcalUncalibratedRecHit& uncalibRH, EcalRecHitCollection& result) override;
0044 
0045 protected:
0046   void insertRecHit(const EcalRecHit& hit, EcalRecHitCollection& collection);
0047   float recCheckCalib(float energy, int ieta);
0048   bool alreadyInserted(const DetId& id);
0049   float estimateEnergy(int ieta,
0050                        EcalRecHitCollection* hits,
0051                        const std::set<DetId>& sId,
0052                        const std::vector<DetId>& vId,
0053                        const EcalTPGScale& tpgscale);
0054   bool checkChannelStatus(const DetId& id, const std::vector<int>& statusestoexclude);
0055 
0056   edm::ESHandle<EcalLaserDbService> laser;
0057   edm::ESGetToken<EcalLaserDbService, EcalLaserDbRecord> laserToken_;
0058 
0059   // isolated dead channels
0060   edm::ESHandle<CaloTopology> caloTopology_;
0061   edm::ESHandle<CaloGeometry> caloGeometry_;
0062   edm::ESHandle<EcalChannelStatus> chStatus_;
0063   edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopologyToken_;
0064   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0065   edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> chStatusToken_;
0066 
0067   double singleRecoveryThreshold_;
0068   double sum8RecoveryThreshold_;
0069   std::string singleRecoveryMethod_;
0070   bool killDeadChannels_;
0071 
0072   bool recoverEBIsolatedChannels_;
0073   bool recoverEEIsolatedChannels_;
0074   bool recoverEBVFE_;
0075   bool recoverEEVFE_;
0076   bool recoverEBFE_;
0077   bool recoverEEFE_;
0078 
0079   // list of channel statuses for which recovery in EE should
0080   // not be attempted
0081   std::vector<int> dbStatusToBeExcludedEE_;
0082   std::vector<int> dbStatusToBeExcludedEB_;
0083 
0084   const edm::EventSetup* eventSetup_ = nullptr;
0085   // dead FE
0086   EcalTPGScale::Tokens ecalScaleTokens_;
0087   edm::EDGetTokenT<EcalTrigPrimDigiCollection> tpDigiToken_;
0088   edm::ESHandle<EcalElectronicsMapping> pEcalMapping_;
0089   const EcalElectronicsMapping* ecalMapping_;
0090   double logWarningEtThreshold_EB_FE_;
0091   double logWarningEtThreshold_EE_FE_;
0092 
0093   edm::ESHandle<EcalTrigTowerConstituentsMap> ttMap_;
0094 
0095   edm::ESHandle<CaloSubdetectorGeometry> pEBGeom_;
0096   const CaloSubdetectorGeometry* ebGeom_;
0097   const CaloGeometry* geo_;
0098   edm::ESGetToken<EcalElectronicsMapping, EcalMappingRcd> pEcalMappingToken_;
0099   edm::ESGetToken<EcalTrigTowerConstituentsMap, IdealGeometryRecord> ttMapToken_;
0100   edm::ESGetToken<CaloSubdetectorGeometry, EcalBarrelGeometryRecord> pEBGeomToken_;
0101   std::unique_ptr<EcalRecHitSimpleAlgo> rechitMaker_;
0102 
0103   std::set<DetId> recoveredDetIds_EB_;
0104   std::set<DetId> recoveredDetIds_EE_;
0105 
0106   EcalTPGScale::Tokens tpgscaleTokens_;
0107 
0108   EcalDeadChannelRecoveryAlgos<EBDetId> ebDeadChannelCorrector;
0109   EcalDeadChannelRecoveryAlgos<EEDetId> eeDeadChannelCorrector;
0110 };
0111 
0112 EcalRecHitWorkerRecover::EcalRecHitWorkerRecover(const edm::ParameterSet& ps, edm::ConsumesCollector& c)
0113     : EcalRecHitWorkerBaseClass(ps, c), ecalScaleTokens_(c), tpgscaleTokens_(c) {
0114   rechitMaker_ = std::make_unique<EcalRecHitSimpleAlgo>();
0115   // isolated channel recovery
0116   singleRecoveryMethod_ = ps.getParameter<std::string>("singleChannelRecoveryMethod");
0117   singleRecoveryThreshold_ = ps.getParameter<double>("singleChannelRecoveryThreshold");
0118   sum8RecoveryThreshold_ = ps.getParameter<double>("sum8ChannelRecoveryThreshold");
0119   killDeadChannels_ = ps.getParameter<bool>("killDeadChannels");
0120   recoverEBIsolatedChannels_ = ps.getParameter<bool>("recoverEBIsolatedChannels");
0121   recoverEEIsolatedChannels_ = ps.getParameter<bool>("recoverEEIsolatedChannels");
0122   recoverEBVFE_ = ps.getParameter<bool>("recoverEBVFE");
0123   recoverEEVFE_ = ps.getParameter<bool>("recoverEEVFE");
0124   recoverEBFE_ = ps.getParameter<bool>("recoverEBFE");
0125   recoverEEFE_ = ps.getParameter<bool>("recoverEEFE");
0126   laserToken_ = c.esConsumes<EcalLaserDbService, EcalLaserDbRecord>();
0127   caloTopologyToken_ = c.esConsumes<CaloTopology, CaloTopologyRecord>();
0128   pEcalMappingToken_ = c.esConsumes<EcalElectronicsMapping, EcalMappingRcd>();
0129   pEBGeomToken_ = c.esConsumes<CaloSubdetectorGeometry, EcalBarrelGeometryRecord>(edm::ESInputTag("", "EcalBarrel"));
0130   caloGeometryToken_ = c.esConsumes<CaloGeometry, CaloGeometryRecord>();
0131   chStatusToken_ = c.esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0132   ttMapToken_ = c.esConsumes<EcalTrigTowerConstituentsMap, IdealGeometryRecord>();
0133 
0134   dbStatusToBeExcludedEE_ = ps.getParameter<std::vector<int> >("dbStatusToBeExcludedEE");
0135   dbStatusToBeExcludedEB_ = ps.getParameter<std::vector<int> >("dbStatusToBeExcludedEB");
0136 
0137   logWarningEtThreshold_EB_FE_ = ps.getParameter<double>("logWarningEtThreshold_EB_FE");
0138   logWarningEtThreshold_EE_FE_ = ps.getParameter<double>("logWarningEtThreshold_EE_FE");
0139 
0140   tpDigiToken_ =
0141       c.consumes<EcalTrigPrimDigiCollection>(ps.getParameter<edm::InputTag>("triggerPrimitiveDigiCollection"));
0142 
0143   if (recoverEBIsolatedChannels_ && singleRecoveryMethod_ == "BDTG")
0144     ebDeadChannelCorrector.setParameters(ps);
0145 }
0146 
0147 void EcalRecHitWorkerRecover::set(const edm::EventSetup& es) {
0148   laser = es.getHandle(laserToken_);
0149   caloTopology_ = es.getHandle(caloTopologyToken_);
0150   pEcalMapping_ = es.getHandle(pEcalMappingToken_);
0151   ecalMapping_ = pEcalMapping_.product();
0152   // geometry...
0153   pEBGeom_ = es.getHandle(pEBGeomToken_);
0154   caloGeometry_ = es.getHandle(caloGeometryToken_);
0155   chStatus_ = es.getHandle(chStatusToken_);
0156   geo_ = caloGeometry_.product();
0157   ebGeom_ = pEBGeom_.product();
0158   ttMap_ = es.getHandle(ttMapToken_);
0159   recoveredDetIds_EB_.clear();
0160   recoveredDetIds_EE_.clear();
0161   eventSetup_ = &es;
0162 }
0163 
0164 bool EcalRecHitWorkerRecover::run(const edm::Event& evt,
0165                                   const EcalUncalibratedRecHit& uncalibRH,
0166                                   EcalRecHitCollection& result) {
0167   DetId detId = uncalibRH.id();
0168   uint32_t flags = (0xF & uncalibRH.flags());
0169 
0170   // get laser coefficient
0171   //float lasercalib = laser->getLaserCorrection( detId, evt.time());
0172 
0173   // killDeadChannels_ = true, means explicitely kill dead channels even if the recovered energies are computed in the code
0174   // if you don't want to store the recovered energies in the rechit you can produce LogWarnings if logWarningEtThreshold_EB(EE)_FE>0
0175   // logWarningEtThreshold_EB(EE)_FE_<0 will not compute the recovered energies at all (faster)
0176 
0177   if (killDeadChannels_) {
0178     if ((flags == EcalRecHitWorkerRecover::EB_single && !recoverEBIsolatedChannels_) ||
0179         (flags == EcalRecHitWorkerRecover::EE_single && !recoverEEIsolatedChannels_) ||
0180         (flags == EcalRecHitWorkerRecover::EB_VFE && !recoverEBVFE_) ||
0181         (flags == EcalRecHitWorkerRecover::EE_VFE && !recoverEEVFE_)) {
0182       EcalRecHit hit(detId, 0., 0., EcalRecHit::kDead);
0183       hit.setFlag(EcalRecHit::kDead);
0184       insertRecHit(hit, result);  // insert trivial rechit with kDead flag
0185       return true;
0186     }
0187     if (flags == EcalRecHitWorkerRecover::EB_FE && !recoverEBFE_) {
0188       EcalTrigTowerDetId ttDetId(((EBDetId)detId).tower());
0189       std::vector<DetId> vid = ttMap_->constituentsOf(ttDetId);
0190       for (std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) {
0191         EcalRecHit hit((*dit), 0., 0., EcalRecHit::kDead);
0192         hit.setFlag(EcalRecHit::kDead);
0193         insertRecHit(hit, result);  // insert trivial rechit with kDead flag
0194       }
0195       if (logWarningEtThreshold_EB_FE_ < 0)
0196         return true;  // if you don't want log warning just return true
0197     }
0198     if (flags == EcalRecHitWorkerRecover::EE_FE && !recoverEEFE_) {
0199       EEDetId id(detId);
0200       EcalScDetId sc(1 + (id.ix() - 1) / 5, 1 + (id.iy() - 1) / 5, id.zside());
0201       std::vector<DetId> eeC;
0202       for (int dx = 1; dx <= 5; ++dx) {
0203         for (int dy = 1; dy <= 5; ++dy) {
0204           int ix = (sc.ix() - 1) * 5 + dx;
0205           int iy = (sc.iy() - 1) * 5 + dy;
0206           int iz = sc.zside();
0207           if (EEDetId::validDetId(ix, iy, iz)) {
0208             eeC.push_back(EEDetId(ix, iy, iz));
0209           }
0210         }
0211       }
0212       for (size_t i = 0; i < eeC.size(); ++i) {
0213         EcalRecHit hit(eeC[i], 0., 0., EcalRecHit::kDead);
0214         hit.setFlag(EcalRecHit::kDead);
0215         insertRecHit(hit, result);  // insert trivial rechit with kDead flag
0216       }
0217       if (logWarningEtThreshold_EE_FE_ < 0)
0218         return true;  // if you don't want log warning just return true
0219     }
0220   }
0221 
0222   if (flags == EcalRecHitWorkerRecover::EB_single) {
0223     // recover as single dead channel
0224     ebDeadChannelCorrector.setCaloTopology(caloTopology_.product());
0225 
0226     // channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE
0227     bool AcceptRecHit = true;
0228     float ebEn = ebDeadChannelCorrector.correct(
0229         detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, &AcceptRecHit);
0230     EcalRecHit hit(detId, ebEn, 0., EcalRecHit::kDead);
0231 
0232     if (hit.energy() != 0 and AcceptRecHit == true) {
0233       hit.setFlag(EcalRecHit::kNeighboursRecovered);
0234     } else {
0235       // recovery failed
0236       hit.setFlag(EcalRecHit::kDead);
0237     }
0238     insertRecHit(hit, result);
0239 
0240   } else if (flags == EcalRecHitWorkerRecover::EE_single) {
0241     // recover as single dead channel
0242     eeDeadChannelCorrector.setCaloTopology(caloTopology_.product());
0243 
0244     // channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE
0245     bool AcceptRecHit = true;
0246     float eeEn = eeDeadChannelCorrector.correct(
0247         detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, &AcceptRecHit);
0248     EcalRecHit hit(detId, eeEn, 0., EcalRecHit::kDead);
0249     if (hit.energy() != 0 and AcceptRecHit == true) {
0250       hit.setFlag(EcalRecHit::kNeighboursRecovered);
0251     } else {
0252       // recovery failed
0253       hit.setFlag(EcalRecHit::kDead);
0254     }
0255     insertRecHit(hit, result);
0256 
0257   } else if (flags == EcalRecHitWorkerRecover::EB_VFE) {
0258     // recover as dead VFE
0259     EcalRecHit hit(detId, 0., 0.);
0260     hit.setFlag(EcalRecHit::kDead);
0261     // recovery not implemented
0262     insertRecHit(hit, result);
0263   } else if (flags == EcalRecHitWorkerRecover::EB_FE) {
0264     // recover as dead TT
0265 
0266     EcalTrigTowerDetId ttDetId(((EBDetId)detId).tower());
0267     edm::Handle<EcalTrigPrimDigiCollection> pTPDigis;
0268     evt.getByToken(tpDigiToken_, pTPDigis);
0269     const EcalTrigPrimDigiCollection* tpDigis = nullptr;
0270     tpDigis = pTPDigis.product();
0271 
0272     EcalTrigPrimDigiCollection::const_iterator tp = tpDigis->find(ttDetId);
0273     // recover the whole trigger tower
0274     if (tp != tpDigis->end()) {
0275       EcalTPGScale ecalScale{ecalScaleTokens_, *eventSetup_};
0276       //std::vector<DetId> vid = ecalMapping_->dccTowerConstituents( ecalMapping_->DCCid( ttDetId ), ecalMapping_->iTT( ttDetId ) );
0277       std::vector<DetId> vid = ttMap_->constituentsOf(ttDetId);
0278       float tpEt = ecalScale.getTPGInGeV(tp->compressedEt(), tp->id());
0279       float tpEtThreshEB = logWarningEtThreshold_EB_FE_;
0280       if (tpEt > tpEtThreshEB) {
0281         edm::LogWarning("EnergyInDeadEB_FE") << "TP energy in the dead TT = " << tpEt << " at " << ttDetId;
0282       }
0283       if (!killDeadChannels_ || recoverEBFE_) {
0284         // democratic energy sharing
0285 
0286         for (std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) {
0287           if (alreadyInserted(*dit))
0288             continue;
0289           float theta = ebGeom_->getGeometry(*dit)->getPosition().theta();
0290           float tpEt = ecalScale.getTPGInGeV(tp->compressedEt(), tp->id());
0291           if (checkChannelStatus(*dit, dbStatusToBeExcludedEB_)) {
0292             EcalRecHit hit(*dit, tpEt / ((float)vid.size()) / sin(theta), 0.);
0293             hit.setFlag(EcalRecHit::kTowerRecovered);
0294             if (tp->compressedEt() == 0xFF)
0295               hit.setFlag(EcalRecHit::kTPSaturated);
0296             if (tp->sFGVB())
0297               hit.setFlag(EcalRecHit::kL1SpikeFlag);
0298             insertRecHit(hit, result);
0299           }
0300         }
0301       } else {
0302         // tp not found => recovery failed
0303         std::vector<DetId> vid = ttMap_->constituentsOf(ttDetId);
0304         for (std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) {
0305           if (alreadyInserted(*dit))
0306             continue;
0307           EcalRecHit hit(*dit, 0., 0.);
0308           hit.setFlag(EcalRecHit::kDead);
0309           insertRecHit(hit, result);
0310         }
0311       }
0312     }
0313   } else if (flags == EcalRecHitWorkerRecover::EE_FE) {
0314     // Structure for recovery:
0315     // ** SC --> EEDetId constituents (eeC) --> associated Trigger Towers (aTT) --> EEDetId constituents (aTTC)
0316     // ** energy for a SC EEDetId = [ sum_aTT(energy) - sum_aTTC(energy) ] / N_eeC
0317     // .. i.e. the total energy of the TTs covering the SC minus
0318     // .. the energy of the recHits in the TTs but not in the SC
0319     //std::vector<DetId> vid = ecalMapping_->dccTowerConstituents( ecalMapping_->DCCid( ttDetId ), ecalMapping_->iTT( ttDetId ) );
0320     // due to lack of implementation of the EcalTrigTowerDetId ix,iy methods in EE we compute Et recovered energies (in EB we compute E)
0321 
0322     EEDetId eeId(detId);
0323     EcalScDetId sc((eeId.ix() - 1) / 5 + 1, (eeId.iy() - 1) / 5 + 1, eeId.zside());
0324     std::set<DetId> eeC;
0325     for (int dx = 1; dx <= 5; ++dx) {
0326       for (int dy = 1; dy <= 5; ++dy) {
0327         int ix = (sc.ix() - 1) * 5 + dx;
0328         int iy = (sc.iy() - 1) * 5 + dy;
0329         int iz = sc.zside();
0330         if (EEDetId::validDetId(ix, iy, iz)) {
0331           EEDetId id(ix, iy, iz);
0332           if (checkChannelStatus(id, dbStatusToBeExcludedEE_)) {
0333             eeC.insert(id);
0334           }  // check status
0335         }
0336       }
0337     }
0338 
0339     edm::Handle<EcalTrigPrimDigiCollection> pTPDigis;
0340     evt.getByToken(tpDigiToken_, pTPDigis);
0341     const EcalTrigPrimDigiCollection* tpDigis = nullptr;
0342     tpDigis = pTPDigis.product();
0343 
0344     // associated trigger towers
0345     std::set<EcalTrigTowerDetId> aTT;
0346     for (std::set<DetId>::const_iterator it = eeC.begin(); it != eeC.end(); ++it) {
0347       aTT.insert(ttMap_->towerOf(*it));
0348     }
0349 
0350     EcalTPGScale tpgscale(tpgscaleTokens_, *eventSetup_);
0351     EcalTPGScale ecalScale(ecalScaleTokens_, *eventSetup_);
0352     // associated trigger towers: total energy
0353     float totE = 0;
0354     // associated trigger towers: EEDetId constituents
0355     std::set<DetId> aTTC;
0356     bool atLeastOneTPSaturated = false;
0357     for (std::set<EcalTrigTowerDetId>::const_iterator it = aTT.begin(); it != aTT.end(); ++it) {
0358       // add the energy of this trigger tower
0359       EcalTrigPrimDigiCollection::const_iterator itTP = tpDigis->find(*it);
0360       if (itTP != tpDigis->end()) {
0361         std::vector<DetId> v = ttMap_->constituentsOf(*it);
0362 
0363         // from the constituents, remove dead channels
0364         std::vector<DetId>::iterator ttcons = v.begin();
0365         while (ttcons != v.end()) {
0366           if (!checkChannelStatus(*ttcons, dbStatusToBeExcludedEE_)) {
0367             ttcons = v.erase(ttcons);
0368           } else {
0369             ++ttcons;
0370           }
0371         }  // while
0372 
0373         if (itTP->compressedEt() == 0xFF) {  // In the case of a saturated trigger tower, a fraction
0374           atLeastOneTPSaturated =
0375               true;  //of the saturated energy is put in: number of xtals in dead region/total xtals in TT *63.75
0376 
0377           //Alternative recovery algorithm that I will now investigate.
0378           //Estimate energy sums the energy in the working channels, then decides how much energy
0379           //to put here depending on that. Duncan 20101203
0380 
0381           totE += estimateEnergy(itTP->id().ietaAbs(), &result, eeC, v, tpgscale);
0382 
0383           /* 
0384                          These commented out lines use
0385                          64GeV*fraction of the TT overlapping the dead FE
0386                         
0387                       int count = 0;
0388                       for (std::vector<DetId>::const_iterator idsit = v.begin(); idsit != v.end(); ++ idsit){
0389                       std::set<DetId>::const_iterator itFind = eeC.find(*idsit);
0390                       if (itFind != eeC.end())
0391                       ++count;
0392                       }
0393                       //std::cout << count << ", " << v.size() << std::endl;
0394                       totE+=((float)count/(float)v.size())* ((it->ietaAbs()>26)?2*ecalScale_.getTPGInGeV( itTP->compressedEt(), itTP->id() ):ecalScale_.getTPGInGeV( itTP->compressedEt(), itTP->id() ));*/
0395         } else {
0396           totE += ((it->ietaAbs() > 26) ? 2 : 1) * ecalScale.getTPGInGeV(itTP->compressedEt(), itTP->id());
0397         }
0398 
0399         // get the trigger tower constituents
0400 
0401         if (itTP->compressedEt() == 0) {  // If there's no energy in TT, the constituents are removed from the recovery.
0402           for (size_t i = 0; i < v.size(); ++i)
0403             eeC.erase(v[i]);
0404         } else if (itTP->compressedEt() != 0xFF) {
0405           //If it's saturated the energy has already been determined, so we do not want to subtract any channels
0406           for (size_t j = 0; j < v.size(); ++j) {
0407             aTTC.insert(v[j]);
0408           }
0409         }
0410       }
0411     }
0412     // remove crystals of dead SC
0413     // (this step is not needed if sure that SC crystals are not
0414     // in the recHit collection)
0415 
0416     for (std::set<DetId>::const_iterator it = eeC.begin(); it != eeC.end(); ++it) {
0417       aTTC.erase(*it);
0418     }
0419     // compute the total energy for the dead SC
0420     const EcalRecHitCollection* hits = &result;
0421     for (std::set<DetId>::const_iterator it = aTTC.begin(); it != aTTC.end(); ++it) {
0422       EcalRecHitCollection::const_iterator jt = hits->find(*it);
0423       if (jt != hits->end()) {
0424         float energy = jt->energy();  // Correct conversion to Et
0425         float eta = geo_->getPosition(jt->id()).eta();
0426         float pf = 1.0 / cosh(eta);
0427         // use Et instead of E, consistent with the Et estimation of the associated TT
0428         totE -= energy * pf;
0429       }
0430     }
0431 
0432     float scEt = totE;
0433     float scEtThreshEE = logWarningEtThreshold_EE_FE_;
0434     if (scEt > scEtThreshEE) {
0435       edm::LogWarning("EnergyInDeadEE_FE") << "TP energy in the dead TT = " << scEt << " at " << sc;
0436     }
0437 
0438     // assign the energy to the SC crystals
0439     if (!killDeadChannels_ || recoverEEFE_) {  // if eeC is empty, i.e. there are no hits
0440                                                // in the tower, nothing is returned. No negative values from noise.
0441       for (std::set<DetId>::const_iterator it = eeC.begin(); it != eeC.end(); ++it) {
0442         float eta = geo_->getPosition(*it).eta();  //Convert back to E from Et for the recovered hits
0443         float pf = 1.0 / cosh(eta);
0444         EcalRecHit hit(*it, totE / ((float)eeC.size() * pf), 0);
0445 
0446         if (atLeastOneTPSaturated)
0447           hit.setFlag(EcalRecHit::kTPSaturated);
0448         hit.setFlag(EcalRecHit::kTowerRecovered);
0449         insertRecHit(hit, result);
0450 
0451       }  // for
0452     }    // if
0453   }
0454   return true;
0455 }
0456 
0457 float EcalRecHitWorkerRecover::estimateEnergy(int ieta,
0458                                               EcalRecHitCollection* hits,
0459                                               const std::set<DetId>& sId,
0460                                               const std::vector<DetId>& vId,
0461                                               const EcalTPGScale& tpgscale) {
0462   float xtalE = 0;
0463   int count = 0;
0464   for (std::vector<DetId>::const_iterator vIdit = vId.begin(); vIdit != vId.end(); ++vIdit) {
0465     std::set<DetId>::const_iterator sIdit = sId.find(*vIdit);
0466     if (sIdit == sId.end()) {
0467       float energy = hits->find(*vIdit)->energy();
0468       float eta = geo_->getPosition(*vIdit).eta();
0469       float pf = 1.0 / cosh(eta);
0470       xtalE += energy * pf;
0471       count++;
0472     }
0473   }
0474 
0475   if (count == 0) {  // If there are no overlapping crystals return saturated value.
0476 
0477     double etsat = tpgscale.getTPGInGeV(0xFF,
0478                                         ttMap_->towerOf(*vId.begin()));  // get saturation value for the first
0479                                                                          // constituent, for the others it's the same
0480 
0481     return etsat / cosh(ieta) * (ieta > 26 ? 2 : 1);  // account for duplicated TT in EE for ieta>26
0482   } else
0483     return xtalE * ((vId.size() / (float)count) - 1) * (ieta > 26 ? 2 : 1);
0484 }
0485 
0486 void EcalRecHitWorkerRecover::insertRecHit(const EcalRecHit& hit, EcalRecHitCollection& collection) {
0487   // skip already inserted DetId's and raise a log warning
0488   if (alreadyInserted(hit.id())) {
0489     edm::LogWarning("EcalRecHitWorkerRecover") << "DetId already recovered! Skipping...";
0490     return;
0491   }
0492   EcalRecHitCollection::iterator it = collection.find(hit.id());
0493   if (it == collection.end()) {
0494     // insert the hit in the collection
0495     collection.push_back(hit);
0496   } else {
0497     // overwrite existing recHit
0498     *it = hit;
0499   }
0500   if (hit.id().subdetId() == EcalBarrel) {
0501     recoveredDetIds_EB_.insert(hit.id());
0502   } else if (hit.id().subdetId() == EcalEndcap) {
0503     recoveredDetIds_EE_.insert(hit.id());
0504   } else {
0505     edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << hit.id().rawId();
0506   }
0507 }
0508 
0509 bool EcalRecHitWorkerRecover::alreadyInserted(const DetId& id) {
0510   bool res = false;
0511   if (id.subdetId() == EcalBarrel) {
0512     res = (recoveredDetIds_EB_.find(id) != recoveredDetIds_EB_.end());
0513   } else if (id.subdetId() == EcalEndcap) {
0514     res = (recoveredDetIds_EE_.find(id) != recoveredDetIds_EE_.end());
0515   } else {
0516     edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << id.rawId();
0517   }
0518   return res;
0519 }
0520 
0521 // In the future, this will be used to calibrate the TT energy. There is a dependance on
0522 // eta at lower energies that can be corrected for here after more validation.
0523 float EcalRecHitWorkerRecover::recCheckCalib(float eTT, int ieta) { return eTT; }
0524 
0525 // return false is the channel has status  in the list of statusestoexclude
0526 // true otherwise (channel ok)
0527 // Careful: this function works on raw (encoded) channel statuses
0528 bool EcalRecHitWorkerRecover::checkChannelStatus(const DetId& id, const std::vector<int>& statusestoexclude) {
0529   if (!chStatus_.isValid())
0530     edm::LogError("ObjectNotFound") << "Channel Status not set";
0531 
0532   EcalChannelStatus::const_iterator chIt = chStatus_->find(id);
0533   uint16_t dbStatus = 0;
0534   if (chIt != chStatus_->end()) {
0535     dbStatus = chIt->getEncodedStatusCode();
0536   } else {
0537     edm::LogError("ObjectNotFound") << "No channel status found for xtal " << id.rawId()
0538                                     << "! something wrong with EcalChannelStatus in your DB? ";
0539   }
0540 
0541   for (std::vector<int>::const_iterator status = statusestoexclude.begin(); status != statusestoexclude.end();
0542        ++status) {
0543     if (*status == dbStatus)
0544       return false;
0545   }
0546 
0547   return true;
0548 }
0549 
0550 #include "FWCore/Framework/interface/MakerMacros.h"
0551 #include "RecoLocalCalo/EcalRecProducers/interface/EcalRecHitWorkerFactory.h"
0552 DEFINE_EDM_PLUGIN(EcalRecHitWorkerFactory, EcalRecHitWorkerRecover, "EcalRecHitWorkerRecover");
[ 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!