File indexing completed on 2024-04-06 12:24:41
0001 #include "CondFormats/DataRecord/interface/ESChannelStatusRcd.h"
0002 #include "CondFormats/DataRecord/interface/ESEEIntercalibConstantsRcd.h"
0003 #include "CondFormats/DataRecord/interface/ESGainRcd.h"
0004 #include "CondFormats/DataRecord/interface/ESMIPToGeVConstantRcd.h"
0005 #include "CondFormats/DataRecord/interface/ESMissingEnergyCalibrationRcd.h"
0006 #include "CondFormats/ESObjects/interface/ESChannelStatus.h"
0007 #include "CondFormats/ESObjects/interface/ESEEIntercalibConstants.h"
0008 #include "CondFormats/ESObjects/interface/ESGain.h"
0009 #include "CondFormats/ESObjects/interface/ESMIPToGeVConstant.h"
0010 #include "CondFormats/ESObjects/interface/ESMissingEnergyCalibration.h"
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0013 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0014 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0015 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0016 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0017 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0018 #include "DataFormats/EgammaReco/interface/PreshowerCluster.h"
0019 #include "DataFormats/EgammaReco/interface/PreshowerClusterFwd.h"
0020 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0021 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0022 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0023 #include "DataFormats/Math/interface/Point3D.h"
0024 #include "FWCore/Framework/interface/ESHandle.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/Framework/interface/stream/EDProducer.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/ESGetToken.h"
0032 #include "FWCore/Utilities/interface/Exception.h"
0033 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0034 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0035 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0036 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0037 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0038 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0039 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0040 #include "RecoEcal/EgammaClusterAlgos/interface/PreshowerClusterAlgo.h"
0041
0042 #include <fstream>
0043 #include <memory>
0044 #include <vector>
0045
0046 class PreshowerClusterProducer : public edm::stream::EDProducer<> {
0047 public:
0048 typedef math::XYZPoint Point;
0049
0050 explicit PreshowerClusterProducer(const edm::ParameterSet& ps);
0051
0052 ~PreshowerClusterProducer() override;
0053
0054 void produce(edm::Event& evt, const edm::EventSetup& es) override;
0055 void set(const edm::EventSetup& es);
0056
0057 private:
0058 int nEvt_;
0059
0060
0061 edm::EDGetTokenT<EcalRecHitCollection> preshHitsToken_;
0062 edm::EDGetTokenT<reco::SuperClusterCollection> endcapSClusterToken_;
0063
0064
0065 std::string preshClusterCollectionX_;
0066 std::string preshClusterCollectionY_;
0067
0068 int preshNclust_;
0069 float preshClustECut;
0070 double etThresh_;
0071
0072
0073 std::string assocSClusterCollection_;
0074
0075 edm::ESHandle<ESGain> esgain_;
0076 edm::ESHandle<ESMIPToGeVConstant> esMIPToGeV_;
0077 edm::ESHandle<ESEEIntercalibConstants> esEEInterCalib_;
0078 edm::ESHandle<ESMissingEnergyCalibration> esMissingECalib_;
0079 edm::ESHandle<ESChannelStatus> esChannelStatus_;
0080 edm::ESGetToken<ESGain, ESGainRcd> esGainToken_;
0081 edm::ESGetToken<ESMIPToGeVConstant, ESMIPToGeVConstantRcd> esMIPToGeVToken_;
0082 edm::ESGetToken<ESEEIntercalibConstants, ESEEIntercalibConstantsRcd> esEEInterCalibToken_;
0083 edm::ESGetToken<ESMissingEnergyCalibration, ESMissingEnergyCalibrationRcd> esMissingECalibToken_;
0084 edm::ESGetToken<ESChannelStatus, ESChannelStatusRcd> esChannelStatusToken_;
0085 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0086
0087 double mip_;
0088 double gamma0_;
0089 double gamma1_;
0090 double gamma2_;
0091 double gamma3_;
0092 double alpha0_;
0093 double alpha1_;
0094 double alpha2_;
0095 double alpha3_;
0096 double aEta_[4];
0097 double bEta_[4];
0098
0099 PreshowerClusterAlgo* presh_algo;
0100
0101
0102 };
0103
0104 #include "FWCore/Framework/interface/MakerMacros.h"
0105 DEFINE_FWK_MODULE(PreshowerClusterProducer);
0106
0107 using namespace edm;
0108 using namespace std;
0109
0110 PreshowerClusterProducer::PreshowerClusterProducer(const edm::ParameterSet& ps) {
0111
0112 preshHitsToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("preshRecHitProducer"));
0113
0114
0115 endcapSClusterToken_ =
0116 consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("endcapSClusterProducer"));
0117
0118 esGainToken_ = esConsumes<ESGain, ESGainRcd>();
0119 esMIPToGeVToken_ = esConsumes<ESMIPToGeVConstant, ESMIPToGeVConstantRcd>();
0120 esEEInterCalibToken_ = esConsumes<ESEEIntercalibConstants, ESEEIntercalibConstantsRcd>();
0121 esMissingECalibToken_ = esConsumes<ESMissingEnergyCalibration, ESMissingEnergyCalibrationRcd>();
0122 esChannelStatusToken_ = esConsumes<ESChannelStatus, ESChannelStatusRcd>();
0123 caloGeometryToken_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0124
0125
0126 preshClusterCollectionX_ = ps.getParameter<std::string>("preshClusterCollectionX");
0127 preshClusterCollectionY_ = ps.getParameter<std::string>("preshClusterCollectionY");
0128 preshNclust_ = ps.getParameter<int>("preshNclust");
0129
0130 etThresh_ = ps.getParameter<double>("etThresh");
0131
0132 assocSClusterCollection_ = ps.getParameter<std::string>("assocSClusterCollection");
0133
0134 produces<reco::PreshowerClusterCollection>(preshClusterCollectionX_);
0135 produces<reco::PreshowerClusterCollection>(preshClusterCollectionY_);
0136 produces<reco::SuperClusterCollection>(assocSClusterCollection_);
0137
0138 float preshStripECut = ps.getParameter<double>("preshStripEnergyCut");
0139 int preshSeededNst = ps.getParameter<int>("preshSeededNstrip");
0140 preshClustECut = ps.getParameter<double>("preshClusterEnergyCut");
0141
0142 presh_algo = new PreshowerClusterAlgo(preshStripECut, preshClustECut, preshSeededNst);
0143
0144 nEvt_ = 0;
0145 }
0146
0147 PreshowerClusterProducer::~PreshowerClusterProducer() { delete presh_algo; }
0148
0149 void PreshowerClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0150 edm::Handle<EcalRecHitCollection> pRecHits;
0151 edm::Handle<reco::SuperClusterCollection> pSuperClusters;
0152
0153
0154 edm::ESHandle<CaloGeometry> geoHandle = es.getHandle(caloGeometryToken_);
0155
0156
0157 set(es);
0158 const ESChannelStatus* channelStatus = esChannelStatus_.product();
0159
0160 const CaloSubdetectorGeometry* geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0161
0162
0163 auto clusters_p1 = std::make_unique<reco::PreshowerClusterCollection>();
0164 auto clusters_p2 = std::make_unique<reco::PreshowerClusterCollection>();
0165
0166 auto superclusters_p = std::make_unique<reco::SuperClusterCollection>();
0167
0168 std::unique_ptr<CaloSubdetectorTopology> topology_p;
0169 if (geometry_p)
0170 topology_p = std::make_unique<EcalPreshowerTopology>();
0171
0172
0173 evt.getByToken(endcapSClusterToken_, pSuperClusters);
0174 const reco::SuperClusterCollection* SClusts = pSuperClusters.product();
0175
0176
0177 evt.getByToken(preshHitsToken_, pRecHits);
0178
0179 const EcalRecHitCollection* rechits = pRecHits.product();
0180
0181 LogTrace("EcalClusters") << "PreshowerClusterProducerInfo: ### Total # of preshower RecHits: " << rechits->size();
0182
0183
0184 std::map<DetId, EcalRecHit> rechits_map;
0185 EcalRecHitCollection::const_iterator it;
0186 for (it = rechits->begin(); it != rechits->end(); it++) {
0187
0188 if (it->recoFlag() == 1 || it->recoFlag() == 14 || (it->recoFlag() <= 10 && it->recoFlag() >= 5))
0189 continue;
0190
0191 rechits_map.insert(std::make_pair(it->id(), *it));
0192 }
0193
0194 std::set<DetId> used_strips;
0195 used_strips.clear();
0196 LogTrace("EcalClusters") << "PreshowerClusterProducerInfo: ### rechits_map of size " << rechits_map.size()
0197 << " was created!";
0198
0199 reco::PreshowerClusterCollection clusters1, clusters2;
0200 reco::SuperClusterCollection new_SC;
0201
0202
0203 reco::SuperClusterCollection::const_iterator it_super;
0204 int isc = 0;
0205 for (it_super = SClusts->begin(); it_super != SClusts->end(); ++it_super) {
0206 float e1 = 0;
0207 float e2 = 0;
0208 float deltaE = 0;
0209
0210 reco::CaloClusterPtrVector new_BC;
0211 ++isc;
0212 LogTrace("EcalClusters") << " superE = " << it_super->energy() << " superETA = " << it_super->eta()
0213 << " superPHI = " << it_super->phi();
0214
0215 int nBC = 0;
0216 int condP1 = 1;
0217 int condP2 = 1;
0218 reco::CaloCluster_iterator bc_iter = it_super->clustersBegin();
0219 for (; bc_iter != it_super->clustersEnd(); ++bc_iter) {
0220 if (geometry_p) {
0221
0222 double X = (*bc_iter)->x();
0223 double Y = (*bc_iter)->y();
0224 double Z = (*bc_iter)->z();
0225 const GlobalPoint point(X, Y, Z);
0226
0227 DetId tmp1 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 1);
0228 DetId tmp2 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 2);
0229 ESDetId strip1 = (tmp1 == DetId(0)) ? ESDetId(0) : ESDetId(tmp1);
0230 ESDetId strip2 = (tmp2 == DetId(0)) ? ESDetId(0) : ESDetId(tmp2);
0231
0232 if (nBC == 0) {
0233 if (strip1 != ESDetId(0) && strip2 != ESDetId(0)) {
0234 ESChannelStatusMap::const_iterator status_p1 = channelStatus->getMap().find(strip1);
0235 ESChannelStatusMap::const_iterator status_p2 = channelStatus->getMap().find(strip2);
0236 if (status_p1->getStatusCode() == 1)
0237 condP1 = 0;
0238 if (status_p2->getStatusCode() == 1)
0239 condP2 = 0;
0240 } else if (strip1 == ESDetId(0))
0241 condP1 = 0;
0242 else if (strip2 == ESDetId(0))
0243 condP2 = 0;
0244 }
0245
0246
0247 for (int i = 0; i < preshNclust_; i++) {
0248 reco::PreshowerCluster cl1 =
0249 presh_algo->makeOneCluster(strip1, &used_strips, &rechits_map, geometry_p, topology_p.get());
0250 cl1.setBCRef(*bc_iter);
0251 if (cl1.energy() > preshClustECut) {
0252 clusters1.push_back(cl1);
0253 e1 += cl1.energy();
0254 }
0255 reco::PreshowerCluster cl2 =
0256 presh_algo->makeOneCluster(strip2, &used_strips, &rechits_map, geometry_p, topology_p.get());
0257 cl2.setBCRef(*bc_iter);
0258
0259 if (cl2.energy() > preshClustECut) {
0260 clusters2.push_back(cl2);
0261 e2 += cl2.energy();
0262 }
0263 }
0264 }
0265
0266 new_BC.push_back(*bc_iter);
0267 nBC++;
0268 }
0269
0270 LogTrace("EcalClusters") << " For SC #" << isc - 1 << ", containing " << it_super->clustersSize()
0271 << " basic clusters, PreshowerClusterAlgo made " << clusters1.size() << " in X plane and "
0272 << clusters2.size() << " in Y plane "
0273 << " preshower clusters ";
0274
0275
0276 if (e1 + e2 > 1.0e-10) {
0277 e1 = e1 / mip_;
0278 e2 = e2 / mip_;
0279
0280 if (condP1 == 1 && condP2 == 1) {
0281 deltaE = gamma0_ * (e1 + alpha0_ * e2);
0282 } else if (condP1 == 1 && condP2 == 0) {
0283 deltaE = gamma1_ * (e1 + alpha1_ * e2);
0284 } else if (condP1 == 0 && condP2 == 1) {
0285 deltaE = gamma2_ * (e1 + alpha2_ * e2);
0286 } else if (condP1 == 0 && condP2 == 0) {
0287 deltaE = gamma3_ * (e1 + alpha3_ * e2);
0288 }
0289 }
0290
0291
0292 float E = it_super->energy() + deltaE;
0293
0294 LogTrace("EcalClusters") << " Creating corrected SC ";
0295
0296 reco::SuperCluster sc(E, it_super->position(), it_super->seed(), new_BC, deltaE);
0297 if (condP1 == 1 && condP2 == 1)
0298 sc.setPreshowerPlanesStatus(0);
0299 else if (condP1 == 1 && condP2 == 0)
0300 sc.setPreshowerPlanesStatus(1);
0301 else if (condP1 == 0 && condP2 == 1)
0302 sc.setPreshowerPlanesStatus(2);
0303 else if (condP1 == 0 && condP2 == 0)
0304 sc.setPreshowerPlanesStatus(3);
0305
0306 if (sc.energy() * sin(sc.position().theta()) > etThresh_)
0307 new_SC.push_back(sc);
0308 LogTrace("EcalClusters") << " SuperClusters energies: new E = " << sc.energy()
0309 << " vs. old E =" << it_super->energy();
0310
0311 }
0312
0313
0314 clusters_p1->assign(clusters1.begin(), clusters1.end());
0315 clusters_p2->assign(clusters2.begin(), clusters2.end());
0316
0317 evt.put(std::move(clusters_p1), preshClusterCollectionX_);
0318 evt.put(std::move(clusters_p2), preshClusterCollectionY_);
0319 LogTrace("EcalClusters") << "Preshower clusters added to the event";
0320
0321
0322 superclusters_p->assign(new_SC.begin(), new_SC.end());
0323 evt.put(std::move(superclusters_p), assocSClusterCollection_);
0324 LogTrace("EcalClusters") << "Corrected SClusters added to the event";
0325
0326 nEvt_++;
0327 }
0328
0329 void PreshowerClusterProducer::set(const edm::EventSetup& es) {
0330 esgain_ = es.getHandle(esGainToken_);
0331 const ESGain* gain = esgain_.product();
0332
0333 double ESGain = gain->getESGain();
0334
0335 esMIPToGeV_ = es.getHandle(esMIPToGeVToken_);
0336 const ESMIPToGeVConstant* mipToGeV = esMIPToGeV_.product();
0337
0338 mip_ = (ESGain == 1) ? mipToGeV->getESValueLow() : mipToGeV->getESValueHigh();
0339
0340 esChannelStatus_ = es.getHandle(esChannelStatusToken_);
0341
0342 esEEInterCalib_ = es.getHandle(esEEInterCalibToken_);
0343 const ESEEIntercalibConstants* esEEInterCalib = esEEInterCalib_.product();
0344
0345
0346 gamma0_ = (ESGain == 1) ? 0.02 : esEEInterCalib->getGammaHigh0();
0347 alpha0_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow0() : esEEInterCalib->getAlphaHigh0();
0348
0349
0350 gamma1_ = (ESGain == 1) ? (0.02 * esEEInterCalib->getGammaLow1()) : esEEInterCalib->getGammaHigh1();
0351 alpha1_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow1() : esEEInterCalib->getAlphaHigh1();
0352
0353
0354 gamma2_ = (ESGain == 1) ? (0.02 * esEEInterCalib->getGammaLow2()) : esEEInterCalib->getGammaHigh2();
0355 alpha2_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow2() : esEEInterCalib->getAlphaHigh2();
0356
0357
0358 gamma3_ = (ESGain == 1) ? 0.02 : esEEInterCalib->getGammaHigh3();
0359 alpha3_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow3() : esEEInterCalib->getAlphaHigh3();
0360
0361 esMissingECalib_ = es.getHandle(esMissingECalibToken_);
0362 const ESMissingEnergyCalibration* esMissingECalib = esMissingECalib_.product();
0363
0364
0365 aEta_[0] = esMissingECalib->getConstAEta0();
0366 bEta_[0] = esMissingECalib->getConstBEta0();
0367
0368
0369 aEta_[1] = esMissingECalib->getConstAEta1();
0370 bEta_[1] = esMissingECalib->getConstBEta1();
0371
0372
0373 aEta_[2] = esMissingECalib->getConstAEta2();
0374 bEta_[2] = esMissingECalib->getConstBEta2();
0375
0376
0377 aEta_[3] = esMissingECalib->getConstAEta3();
0378 bEta_[3] = esMissingECalib->getConstBEta3();
0379 }