EgammaCoreTools/src/EcalClusterLazyTools.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002
0003 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
0004 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0005
0006 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0007 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0008
0009 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0010 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0011 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0012 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0013 #include "RecoCaloTools/Navigation/interface/EcalPreshowerNavigator.h"
0014
0015 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0016 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0017
0018 EcalClusterLazyToolsBase::EcalClusterLazyToolsBase(const edm::Event &ev,
0019                                                    ESData const &esData,
0020                                                    edm::EDGetTokenT<EcalRecHitCollection> token1,
0021                                                    edm::EDGetTokenT<EcalRecHitCollection> token2,
0022                                                    std::optional<edm::EDGetTokenT<EcalRecHitCollection>> token3)
0023     : geometry_(&esData.caloGeometry),
0024       topology_(&esData.caloTopology),
0025       ebRecHits_(&edm::get(ev, token1)),
0026       eeRecHits_(&edm::get(ev, token2)) {
0027   if (token3) {
0028     if (geometry_->getSubdetectorGeometry(DetId::Ecal, EcalPreshower)) {
0029       ecalPS_topology_ = std::make_unique<EcalPreshowerTopology>();
0030     } else {
0031       edm::LogInfo("subdetector geometry not available") << "EcalPreshower geometry is missing" << std::endl;
0032     }
0033
0034     getESRecHits(ev, *token3);
0035   }
0036
0037   // get IC's
0038   ical = &esData.ecalIntercalibConstants;
0039   icalMap = &ical->getMap();
0040   // get ADCtoGeV
0041   agc = &esData.ecalADCToGeV;
0042   // transp corrections
0043   laser = &esData.ecalLaserDbService;
0044 }
0045
0046 void EcalClusterLazyToolsBase::getESRecHits(const edm::Event &ev,
0047                                             edm::EDGetTokenT<EcalRecHitCollection> const &esRecHitsToken) {
0048   auto pESRecHits = ev.getHandle(esRecHitsToken);
0049   esRecHits_ = pESRecHits.product();
0050   // make the map of rechits
0051   rechits_map_.clear();
0052   if (pESRecHits.isValid()) {
0053     for (auto it = pESRecHits->begin(); it != pESRecHits->end(); ++it) {
0054       // remove bad ES rechits
0055       std::vector<int> badf = {
0056           EcalRecHit::ESFlags::kESDead,  // 1
0057           EcalRecHit::ESFlags::kESTwoGoodRatios,
0058           EcalRecHit::ESFlags::kESBadRatioFor12,  // 5
0059           EcalRecHit::ESFlags::kESBadRatioFor23Upper,
0060           EcalRecHit::ESFlags::kESBadRatioFor23Lower,
0061           EcalRecHit::ESFlags::kESTS1Largest,
0062           EcalRecHit::ESFlags::kESTS3Largest,
0063           EcalRecHit::ESFlags::kESTS3Negative,  // 10
0064           EcalRecHit::ESFlags::kESTS13Sigmas,   // 14
0065       };
0066
0067       if (it->checkFlags(badf))
0068         continue;
0069
0070       //Make the map of DetID, EcalRecHit pairs
0071       rechits_map_.insert(std::make_pair(it->id(), *it));
0072     }
0073   }
0074 }
0075
0076 const EcalRecHitCollection *EcalClusterLazyToolsBase::getEcalRecHitCollection(const reco::BasicCluster &cluster) const {
0077   if (cluster.size() == 0) {
0078     throw cms::Exception("InvalidCluster") << "The cluster has no crystals!";
0079   }
0080   DetId id = (cluster.hitsAndFractions()[0]).first;  // size is by definition > 0 -- FIXME??
0081   const EcalRecHitCollection *recHits = nullptr;
0082   if (id.subdetId() == EcalBarrel) {
0083     recHits = ebRecHits_;
0084   } else if (id.subdetId() == EcalEndcap) {
0085     recHits = eeRecHits_;
0086   } else {
0087     throw cms::Exception("InvalidSubdetector")
0088         << "The subdetId() " << id.subdetId() << " does not correspond to EcalBarrel neither EcalEndcap";
0089   }
0090   return recHits;
0091 }
0092
0093 // get time of basic cluster seed crystal
0094 float EcalClusterLazyToolsBase::BasicClusterSeedTime(const reco::BasicCluster &cluster) {
0095   const EcalRecHitCollection *recHits = getEcalRecHitCollection(cluster);
0096
0097   DetId id = cluster.seed();
0098   auto theSeedHit = recHits->find(id);
0099   //  std::cout << "the seed of the BC has time: "
0100   //<< (*theSeedHit).time()
0101   //<< "and energy: " << (*theSeedHit).energy() << " collection size: " << recHits->size()
0102   //<< "\n" <<std::endl; // GF debug
0103
0104   return (*theSeedHit).time();
0105 }
0106
0107 // error-weighted average of time from constituents of basic cluster
0108 float EcalClusterLazyToolsBase::BasicClusterTime(const reco::BasicCluster &cluster, const edm::Event &ev) {
0109   auto clusterComponents = (cluster).hitsAndFractions();
0110   //std::cout << "BC has this many components: " << clusterComponents.size() << std::endl; // GF debug
0111
0112   const EcalRecHitCollection *recHits = getEcalRecHitCollection(cluster);
0113   //std::cout << "BasicClusterClusterTime - rechits are this many: " << recHits->size() << std::endl; // GF debug
0114
0115   float weightedTsum = 0;
0116   float sumOfWeights = 0;
0117
0118   for (auto detitr = clusterComponents.begin(); detitr != clusterComponents.end(); detitr++) {
0119     //      EcalRecHitCollection::const_iterator theSeedHit = recHits->find (id); // trash this
0120     auto oneHit = recHits->find((detitr->first));
0121
0122     // in order to get back the ADC counts from the recHit energy, three ingredients are necessary:
0123     // 1) get laser correction coefficient
0124     float lasercalib = 1.;
0125     lasercalib = laser->getLaserCorrection(detitr->first, ev.time());
0126     // 2) get intercalibration
0127     auto icalit = icalMap->find(detitr->first);
0128     EcalIntercalibConstant icalconst = 1.;
0129     if (icalit != icalMap->end()) {
0130       icalconst = (*icalit);
0131       // std::cout << "icalconst set to: " << icalconst << std::endl;
0132     } else {
0133       edm::LogError("EcalClusterLazyTools")
0134           << "No intercalib const found for xtal " << (detitr->first).rawId() << "bailing out";
0135       assert(false);
0136     }
0137     // 3) get adc2GeV
0138     float adcToGeV = 1.;
0139     if ((detitr->first).subdetId() == EcalBarrel)
0140       adcToGeV = float(agc->getEBValue());
0141     else if ((detitr->first).subdetId() == EcalEndcap)
0142       adcToGeV = float(agc->getEEValue());
0143     float adc = 2.;
0144     if (icalconst > 0 && lasercalib > 0 && adcToGeV > 0)
0145       adc = (*oneHit).energy() / (icalconst * lasercalib * adcToGeV);
0146
0147     // don't consider recHits with too little amplitude; take sigma_noise_total into account
0148     if ((detitr->first).subdetId() == EcalBarrel && adc < (1.1 * 20))
0149       continue;
0150     if ((detitr->first).subdetId() == EcalEndcap && adc < (2.2 * 20))
0151       continue;
0152
0153     // count only on rechits whose error is trusted by the method (ratio)
0154     if (!(*oneHit).isTimeErrorValid())
0155       continue;
0156
0157     float timeError = (*oneHit).timeError();
0158     // the constant used to build timeError is largely over-estimated ; remove in quadrature 0.6 and add 0.15 back.
0159     // could be prettier if value of constant term was changed at recHit production level
0160     if (timeError > 0.6)
0161       timeError = sqrt(timeError * timeError - 0.6 * 0.6 + 0.15 * 0.15);
0162     else
0163       timeError = sqrt(timeError * timeError + 0.15 * 0.15);
0164
0165     // do the error weighting
0166     weightedTsum += (*oneHit).time() / (timeError * timeError);
0167     sumOfWeights += 1. / (timeError * timeError);
0168   }
0169
0170   // what if no crytal is available for weighted average?
0171   if (sumOfWeights == 0)
0172     return -999;
0173   else
0174     return (weightedTsum / sumOfWeights);
0175 }
0176
0177 // get BasicClusterSeedTime of the seed basic cluser of the supercluster
0178 float EcalClusterLazyToolsBase::SuperClusterSeedTime(const reco::SuperCluster &cluster) {
0179   return BasicClusterSeedTime((*cluster.seed()));
0180 }
0181
0182 // get Preshower effective sigmaIRIR
0183 float EcalClusterLazyToolsBase::eseffsirir(const reco::SuperCluster &cluster) {
0184   if (!(fabs(cluster.eta()) > 1.6 && fabs(cluster.eta()) < 3.))
0185     return 0.;
0186
0187   if (!ecalPS_topology_)
0188     return 0.;
0189
0190   std::vector<float> phoESHitsIXIX =
0191       getESHits(cluster.x(), cluster.y(), cluster.z(), rechits_map_, geometry_, ecalPS_topology_.get(), 0, 1);
0192   std::vector<float> phoESHitsIYIY =
0193       getESHits(cluster.x(), cluster.y(), cluster.z(), rechits_map_, geometry_, ecalPS_topology_.get(), 0, 2);
0194   float phoESShapeIXIX = getESShape(phoESHitsIXIX);
0195   float phoESShapeIYIY = getESShape(phoESHitsIYIY);
0196
0197   return sqrt(phoESShapeIXIX * phoESShapeIXIX + phoESShapeIYIY * phoESShapeIYIY);
0198 }
0199
0200 // get Preshower effective sigmaIXIX
0201 float EcalClusterLazyToolsBase::eseffsixix(const reco::SuperCluster &cluster) {
0202   if (!(fabs(cluster.eta()) > 1.6 && fabs(cluster.eta()) < 3.))
0203     return 0.;
0204
0205   if (!ecalPS_topology_)
0206     return 0.;
0207
0208   std::vector<float> phoESHitsIXIX =
0209       getESHits(cluster.x(), cluster.y(), cluster.z(), rechits_map_, geometry_, ecalPS_topology_.get(), 0, 1);
0210   float phoESShapeIXIX = getESShape(phoESHitsIXIX);
0211
0212   return phoESShapeIXIX;
0213 }
0214
0215 // get Preshower effective sigmaIYIY
0216 float EcalClusterLazyToolsBase::eseffsiyiy(const reco::SuperCluster &cluster) {
0217   if (!(fabs(cluster.eta()) > 1.6 && fabs(cluster.eta()) < 3.))
0218     return 0.;
0219
0220   if (!ecalPS_topology_)
0221     return 0.;
0222
0223   std::vector<float> phoESHitsIYIY =
0224       getESHits(cluster.x(), cluster.y(), cluster.z(), rechits_map_, geometry_, ecalPS_topology_.get(), 0, 2);
0225   float phoESShapeIYIY = getESShape(phoESHitsIYIY);
0226
0227   return phoESShapeIYIY;
0228 }
0229
0230 // get Preshower Rechits
0231 std::vector<float> EcalClusterLazyToolsBase::getESHits(double X,
0232                                                        double Y,
0233                                                        double Z,
0234                                                        const std::map<DetId, EcalRecHit> &_rechits_map,
0235                                                        const CaloGeometry *geometry,
0236                                                        CaloSubdetectorTopology const *topology_p,
0237                                                        int row,
0238                                                        int plane) {
0239   std::map<DetId, EcalRecHit> rechits_map = _rechits_map;
0240   std::vector<float> esHits;
0241
0242   const GlobalPoint point(X, Y, Z);
0243
0244   const CaloSubdetectorGeometry *geometry_p = geometry->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0245
0246   DetId esId = (dynamic_cast<const EcalPreshowerGeometry *>(geometry_p))->getClosestCellInPlane(point, plane);
0247   ESDetId esDetId = (esId == DetId(0)) ? ESDetId(0) : ESDetId(esId);
0248
0249   std::map<DetId, EcalRecHit>::iterator it;
0250   ESDetId next;
0251   ESDetId strip;
0252   strip = esDetId;
0253
0254   EcalPreshowerNavigator theESNav(strip, topology_p);
0255   theESNav.setHome(strip);
0256
0257   if (row == 1) {
0258     if (plane == 1 && strip != ESDetId(0))
0259       strip = theESNav.north();
0260     if (plane == 2 && strip != ESDetId(0))
0261       strip = theESNav.east();
0262   } else if (row == -1) {
0263     if (plane == 1 && strip != ESDetId(0))
0264       strip = theESNav.south();
0265     if (plane == 2 && strip != ESDetId(0))
0266       strip = theESNav.west();
0267   }
0268
0269   if (strip == ESDetId(0)) {
0270     for (int i = 0; i < 31; ++i)
0271       esHits.push_back(0);
0272   } else {
0273     it = rechits_map.find(strip);
0274     if (it != rechits_map.end() && it->second.energy() > 1.0e-10)
0275       esHits.push_back(it->second.energy());
0276     else
0277       esHits.push_back(0);
0278     //cout<<"center : "<<strip<<" "<<it->second.energy()<<endl;
0279
0280     // Front Plane
0281     if (plane == 1) {
0282       // east road
0283       for (int i = 0; i < 15; ++i) {
0284         next = theESNav.east();
0285         if (next != ESDetId(0)) {
0286           it = rechits_map.find(next);
0287           if (it != rechits_map.end() && it->second.energy() > 1.0e-10)
0288             esHits.push_back(it->second.energy());
0289           else
0290             esHits.push_back(0);
0291           //cout<<"east "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
0292         } else {
0293           for (int j = i; j < 15; j++)
0294             esHits.push_back(0);
0295           break;
0296           //cout<<"east "<<i<<" : "<<next<<" "<<0<<endl;
0297         }
0298       }
0299
0300       // west road
0301       theESNav.setHome(strip);
0302       theESNav.home();
0303       for (int i = 0; i < 15; ++i) {
0304         next = theESNav.west();
0305         if (next != ESDetId(0)) {
0306           it = rechits_map.find(next);
0307           if (it != rechits_map.end() && it->second.energy() > 1.0e-10)
0308             esHits.push_back(it->second.energy());
0309           else
0310             esHits.push_back(0);
0311           //cout<<"west "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
0312         } else {
0313           for (int j = i; j < 15; j++)
0314             esHits.push_back(0);
0315           break;
0316           //cout<<"west "<<i<<" : "<<next<<" "<<0<<endl;
0317         }
0318       }
0319     }  // End of Front Plane
0320
0321     // Rear Plane
0322     if (plane == 2) {
0323       // north road
0324       for (int i = 0; i < 15; ++i) {
0325         next = theESNav.north();
0326         if (next != ESDetId(0)) {
0327           it = rechits_map.find(next);
0328           if (it != rechits_map.end() && it->second.energy() > 1.0e-10)
0329             esHits.push_back(it->second.energy());
0330           else
0331             esHits.push_back(0);
0332           //cout<<"north "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
0333         } else {
0334           for (int j = i; j < 15; j++)
0335             esHits.push_back(0);
0336           break;
0337           //cout<<"north "<<i<<" : "<<next<<" "<<0<<endl;
0338         }
0339       }
0340
0341       // south road
0342       theESNav.setHome(strip);
0343       theESNav.home();
0344       for (int i = 0; i < 15; ++i) {
0345         next = theESNav.south();
0346         if (next != ESDetId(0)) {
0347           it = rechits_map.find(next);
0348           if (it != rechits_map.end() && it->second.energy() > 1.0e-10)
0349             esHits.push_back(it->second.energy());
0350           else
0351             esHits.push_back(0);
0352           //cout<<"south "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
0353         } else {
0354           for (int j = i; j < 15; j++)
0355             esHits.push_back(0);
0356           break;
0357           //cout<<"south "<<i<<" : "<<next<<" "<<0<<endl;
0358         }
0359       }
0360     }  // End of Rear Plane
0361   }    // Fill ES RecHits
0362
0363   return esHits;
0364 }
0365
0366 // get Preshower hit shape
0367 float EcalClusterLazyToolsBase::getESShape(const std::vector<float> &ESHits0) {
0368   const int nBIN = 21;
0369   float esRH[nBIN];
0370   for (int idx = 0; idx < nBIN; idx++) {
0371     esRH[idx] = 0.;
0372   }
0373
0374   for (int ibin = 0; ibin < ((nBIN + 1) / 2); ibin++) {
0375     if (ibin == 0) {
0376       esRH[(nBIN - 1) / 2] = ESHits0[ibin];
0377     } else {
0378       esRH[(nBIN - 1) / 2 + ibin] = ESHits0[ibin];
0379       esRH[(nBIN - 1) / 2 - ibin] = ESHits0[ibin + 15];
0380     }
0381   }
0382
0383   // ---- Effective Energy Deposit Width ---- //
0384   double EffWidthSigmaISIS = 0.;
0385   double totalEnergyISIS = 0.;
0386   double EffStatsISIS = 0.;
0387   for (int id_X = 0; id_X < 21; id_X++) {
0388     totalEnergyISIS += esRH[id_X];
0389     EffStatsISIS += esRH[id_X] * (id_X - 10) * (id_X - 10);
0390   }
0391   EffWidthSigmaISIS = (totalEnergyISIS > 0.) ? sqrt(fabs(EffStatsISIS / totalEnergyISIS)) : 0.;
0392
0393   return EffWidthSigmaISIS;
0394 }