EgammaClusterProducers/src/PreshowerClusterProducer.cc

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_;  // internal counter of events
0059
0060   //clustering parameters:
0061   edm::EDGetTokenT<EcalRecHitCollection> preshHitsToken_;               // name of module/plugin/producer producing hits
0062   edm::EDGetTokenT<reco::SuperClusterCollection> endcapSClusterToken_;  // ditto SuperClusters
0063
0064   // name out output collections
0065   std::string preshClusterCollectionX_;
0066   std::string preshClusterCollectionY_;
0067
0068   int preshNclust_;
0069   float preshClustECut;
0070   double etThresh_;
0071
0072   // association parameters:
0073   std::string assocSClusterCollection_;  // name of super cluster output collection
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;  // algorithm doing the real work
0100                                      // The set of used DetID's
0101   //std::set<DetId> used_strips;
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   // use configuration file to setup input/output collection names
0112   preshHitsToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("preshRecHitProducer"));
0113
0114   // Name of a SuperClusterCollection to make associations:
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   // Output collections:
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   // get the ECAL geometry:
0154   edm::ESHandle<CaloGeometry> geoHandle = es.getHandle(caloGeometryToken_);
0155
0156   // retrieve ES-EE intercalibration constants and channel status
0157   set(es);
0158   const ESChannelStatus* channelStatus = esChannelStatus_.product();
0159
0160   const CaloSubdetectorGeometry* geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0161
0162   // create unique_ptr to a PreshowerClusterCollection
0163   auto clusters_p1 = std::make_unique<reco::PreshowerClusterCollection>();
0164   auto clusters_p2 = std::make_unique<reco::PreshowerClusterCollection>();
0165   // create new collection of corrected super clusters
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   // fetch the product (pSuperClusters)
0173   evt.getByToken(endcapSClusterToken_, pSuperClusters);
0174   const reco::SuperClusterCollection* SClusts = pSuperClusters.product();
0175
0176   // fetch the product (RecHits)
0177   evt.getByToken(preshHitsToken_, pRecHits);
0178   // pointer to the object in the product
0179   const EcalRecHitCollection* rechits = pRecHits.product();  // EcalRecHitCollection hit_collection = *rhcHandle;
0180
0181   LogTrace("EcalClusters") << "PreshowerClusterProducerInfo: ### Total # of preshower RecHits: " << rechits->size();
0182
0183   // make the map of rechits:
0184   std::map<DetId, EcalRecHit> rechits_map;
0185   EcalRecHitCollection::const_iterator it;
0186   for (it = rechits->begin(); it != rechits->end(); it++) {
0187     // remove bad ES rechits
0188     if (it->recoFlag() == 1 || it->recoFlag() == 14 || (it->recoFlag() <= 10 && it->recoFlag() >= 5))
0189       continue;
0190     //Make the map of DetID, EcalRecHit pairs
0191     rechits_map.insert(std::make_pair(it->id(), *it));
0192   }
0193   // The set of used DetID's for a given event:
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;  // output collection of corrected PCs
0200   reco::SuperClusterCollection new_SC;                    // output collection of corrected SCs
0201
0202   //make cycle over super clusters
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;  // 0: dead channel; 1: active channel
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         // Get strip position at intersection point of the line EE - Vertex:
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         // Get a vector of ES clusters (found by the PreshSeeded algorithm) associated with a given EE basic cluster.
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         }  // end of cycle over ES clusters
0264       }
0265
0266       new_BC.push_back(*bc_iter);
0267       nBC++;
0268     }  // end of cycle over BCs
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     // update energy of the SuperCluster
0276     if (e1 + e2 > 1.0e-10) {
0277       e1 = e1 / mip_;  // GeV to #MIPs
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     //corrected Energy
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   }  // end of cycle over SCs
0312
0313   // copy the preshower clusters into collections and put in the Event:
0314   clusters_p1->assign(clusters1.begin(), clusters1.end());
0315   clusters_p2->assign(clusters2.begin(), clusters2.end());
0316   // put collection of preshower clusters to the event
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   // put collection of corrected super clusters to the event
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   // both planes work
0346   gamma0_ = (ESGain == 1) ? 0.02 : esEEInterCalib->getGammaHigh0();
0347   alpha0_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow0() : esEEInterCalib->getAlphaHigh0();
0348
0349   // only first plane works
0350   gamma1_ = (ESGain == 1) ? (0.02 * esEEInterCalib->getGammaLow1()) : esEEInterCalib->getGammaHigh1();
0351   alpha1_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow1() : esEEInterCalib->getAlphaHigh1();
0352
0353   // only second plane works
0354   gamma2_ = (ESGain == 1) ? (0.02 * esEEInterCalib->getGammaLow2()) : esEEInterCalib->getGammaHigh2();
0355   alpha2_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow2() : esEEInterCalib->getAlphaHigh2();
0356
0357   // both planes do not work
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   // |eta| < 1.9
0365   aEta_[0] = esMissingECalib->getConstAEta0();
0366   bEta_[0] = esMissingECalib->getConstBEta0();
0367
0368   // 1.9 < |eta| < 2.1
0369   aEta_[1] = esMissingECalib->getConstAEta1();
0370   bEta_[1] = esMissingECalib->getConstBEta1();
0371
0372   // 2.1 < |eta| < 2.3
0373   aEta_[2] = esMissingECalib->getConstAEta2();
0374   bEta_[2] = esMissingECalib->getConstBEta2();
0375
0376   // 2.3 < |eta| < 2.5
0377   aEta_[3] = esMissingECalib->getConstAEta3();
0378   bEta_[3] = esMissingECalib->getConstBEta3();
0379 }