EgammaClusterProducers/src/PreshowerPhiClusterProducer.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/deltaPhi.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/Exception.h"
0032 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0033 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0034 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0035 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0036 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0037 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0038 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0039 #include "RecoEcal/EgammaClusterAlgos/interface/PreshowerPhiClusterAlgo.h"
0040
0041 #include <fstream>
0042 #include <memory>
0043 #include <vector>
0044
0045 class PreshowerPhiClusterProducer : public edm::stream::EDProducer<> {
0046 public:
0047   typedef math::XYZPoint Point;
0048
0049   explicit PreshowerPhiClusterProducer(const edm::ParameterSet& ps);
0050
0051   ~PreshowerPhiClusterProducer() override;
0052
0053   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0054   void set(const edm::EventSetup& es);
0055
0056 private:
0057   int nEvt_;  // internal counter of events
0058
0059   //clustering parameters:
0060   edm::EDGetTokenT<EcalRecHitCollection> preshHitToken_;                // name of module/plugin/producer
0061                                                                         // 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   // association parameters:
0069   std::string assocSClusterCollection_;  // name of super cluster output collection
0070
0071   edm::ESHandle<ESGain> esgain_;
0072   edm::ESHandle<ESMIPToGeVConstant> esMIPToGeV_;
0073   edm::ESHandle<ESEEIntercalibConstants> esEEInterCalib_;
0074   edm::ESHandle<ESMissingEnergyCalibration> esMissingECalib_;
0075   edm::ESHandle<ESChannelStatus> esChannelStatus_;
0076   edm::ESGetToken<ESGain, ESGainRcd> esGainToken_;
0077   edm::ESGetToken<ESMIPToGeVConstant, ESMIPToGeVConstantRcd> esMIPToGeVToken_;
0078   edm::ESGetToken<ESEEIntercalibConstants, ESEEIntercalibConstantsRcd> esEEInterCalibToken_;
0079   edm::ESGetToken<ESMissingEnergyCalibration, ESMissingEnergyCalibrationRcd> esMissingECalibToken_;
0080   edm::ESGetToken<ESChannelStatus, ESChannelStatusRcd> esChannelStatusToken_;
0081   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0082   double mip_;
0083   double gamma0_;
0084   double gamma1_;
0085   double gamma2_;
0086   double gamma3_;
0087   double alpha0_;
0088   double alpha1_;
0089   double alpha2_;
0090   double alpha3_;
0091   double aEta_[4];
0092   double bEta_[4];
0093
0094   double etThresh_;
0095
0096   PreshowerPhiClusterAlgo* presh_algo;  // algorithm doing the real work
0097                                         // The set of used DetID's
0098   //std::set<DetId> used_strips;
0099
0100   float esPhiClusterDeltaEta_;
0101   float esPhiClusterDeltaPhi_;
0102 };
0103
0104 #include "FWCore/Framework/interface/MakerMacros.h"
0105 DEFINE_FWK_MODULE(PreshowerPhiClusterProducer);
0106
0107 using namespace edm;
0108 using namespace std;
0109
0110 PreshowerPhiClusterProducer::PreshowerPhiClusterProducer(const edm::ParameterSet& ps) {
0111   // use configuration file to setup input/output collection names
0112   preshHitToken_ = 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
0129   assocSClusterCollection_ = ps.getParameter<std::string>("assocSClusterCollection");
0130
0131   produces<reco::PreshowerClusterCollection>(preshClusterCollectionX_);
0132   produces<reco::PreshowerClusterCollection>(preshClusterCollectionY_);
0133   produces<reco::SuperClusterCollection>(assocSClusterCollection_);
0134
0135   float esStripECut = ps.getParameter<double>("esStripEnergyCut");
0136   esPhiClusterDeltaEta_ = ps.getParameter<double>("esPhiClusterDeltaEta");
0137   esPhiClusterDeltaPhi_ = ps.getParameter<double>("esPhiClusterDeltaPhi");
0138
0139   etThresh_ = ps.getParameter<double>("etThresh");
0140
0141   presh_algo = new PreshowerPhiClusterAlgo(esStripECut);
0142 }
0143
0144 PreshowerPhiClusterProducer::~PreshowerPhiClusterProducer() { delete presh_algo; }
0145
0146 void PreshowerPhiClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0147   edm::Handle<EcalRecHitCollection> pRecHits;
0148   edm::Handle<reco::SuperClusterCollection> pSuperClusters;
0149
0150   // get the ECAL geometry:
0151   edm::ESHandle<CaloGeometry> geoHandle = es.getHandle(caloGeometryToken_);
0152
0153   // retrieve ES-EE intercalibration constants and channel status
0154   set(es);
0155   const ESChannelStatus* channelStatus = esChannelStatus_.product();
0156
0157   const CaloSubdetectorGeometry* geometry_p = (geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower));
0158
0159   // create unique_ptr to a PreshowerClusterCollection
0160   auto clusters_p1 = std::make_unique<reco::PreshowerClusterCollection>();
0161   auto clusters_p2 = std::make_unique<reco::PreshowerClusterCollection>();
0162   // create new collection of corrected super clusters
0163   auto superclusters_p = std::make_unique<reco::SuperClusterCollection>();
0164
0165   // fetch the product (pSuperClusters)
0166   evt.getByToken(endcapSClusterToken_, pSuperClusters);
0167   const reco::SuperClusterCollection* SClusts = pSuperClusters.product();
0168
0169   // fetch the product (RecHits)
0170   evt.getByToken(preshHitToken_, pRecHits);
0171   // pointer to the object in the product
0172   const EcalRecHitCollection* rechits = pRecHits.product();  // EcalRecHitCollection hit_collection = *rhcHandle;
0173
0174   LogTrace("EcalClusters") << "PreshowerPhiClusterProducerInfo: ### Total # of preshower RecHits: " << rechits->size();
0175
0176   // make the map of rechits:
0177   std::map<DetId, EcalRecHit> rechits_map;
0178   EcalRecHitCollection::const_iterator it;
0179   for (it = rechits->begin(); it != rechits->end(); it++) {
0180     // remove bad ES rechits
0181     if (it->recoFlag() == 1 || it->recoFlag() == 14 || (it->recoFlag() <= 10 && it->recoFlag() >= 5))
0182       continue;
0183     //Make the map of DetID, EcalRecHit pairs
0184     rechits_map.insert(std::make_pair(it->id(), *it));
0185   }
0186   // The set of used DetID's for a given event:
0187   std::set<DetId> used_strips;
0188   used_strips.clear();
0189   LogTrace("EcalClusters") << "PreshowerPhiClusterProducerInfo: ### rechits_map of size " << rechits_map.size()
0190                            << " was created!";
0191
0192   reco::PreshowerClusterCollection clusters1, clusters2;  // output collection of corrected PCs
0193   reco::SuperClusterCollection new_SC;                    // output collection of corrected SCs
0194
0195   //make cycle over super clusters
0196   reco::SuperClusterCollection::const_iterator it_super;
0197   int isc = 0;
0198   for (it_super = SClusts->begin(); it_super != SClusts->end(); ++it_super) {
0199     float e1 = 0;
0200     float e2 = 0;
0201     float deltaE = 0;
0202
0203     reco::CaloClusterPtrVector new_BC;
0204     ++isc;
0205     LogTrace("EcalClusters") << " superE = " << it_super->energy() << " superETA = " << it_super->eta()
0206                              << " superPHI = " << it_super->phi();
0207
0208     //cout<<"=== new SC ==="<<endl;
0209     //cout<<"superE = "<<it_super->energy()<<" superETA = "<<it_super->eta()<<" superPHI = "<<it_super->phi()<<endl;
0210
0211     int nBC = 0;
0212     int condP1 = 1;  // 0: dead channel; 1: active channel
0213     int condP2 = 1;
0214     float maxDeltaPhi = 0;
0215     float minDeltaPhi = 0;
0216     float refPhi = 0;
0217
0218     reco::CaloCluster_iterator bc_iter = it_super->clustersBegin();
0219     for (; bc_iter != it_super->clustersEnd(); ++bc_iter) {
0220       if (nBC == 0) {
0221         refPhi = (*bc_iter)->phi();
0222       } else {
0223         if (reco::deltaPhi((*bc_iter)->phi(), refPhi) > 0 && reco::deltaPhi((*bc_iter)->phi(), refPhi) > maxDeltaPhi)
0224           maxDeltaPhi = reco::deltaPhi((*bc_iter)->phi(), refPhi);
0225         if (reco::deltaPhi((*bc_iter)->phi(), refPhi) < 0 && reco::deltaPhi((*bc_iter)->phi(), refPhi) < minDeltaPhi)
0226           minDeltaPhi = reco::deltaPhi((*bc_iter)->phi(), refPhi);
0227         //cout<<"delta phi : "<<reco::deltaPhi((*bc_iter)->phi(), refPhi)<<endl;
0228       }
0229       //cout<<"BC : "<<nBC<<" "<<(*bc_iter)->energy()<<" "<<(*bc_iter)->eta()<<" "<<(*bc_iter)->phi()<<endl;
0230       nBC++;
0231     }
0232     maxDeltaPhi += esPhiClusterDeltaPhi_;
0233     minDeltaPhi -= esPhiClusterDeltaPhi_;
0234
0235     nBC = 0;
0236     for (bc_iter = it_super->clustersBegin(); bc_iter != it_super->clustersEnd(); ++bc_iter) {
0237       if (geometry_p) {
0238         // Get strip position at intersection point of the line EE - Vertex:
0239         double X = (*bc_iter)->x();
0240         double Y = (*bc_iter)->y();
0241         double Z = (*bc_iter)->z();
0242         const GlobalPoint point(X, Y, Z);
0243
0244         DetId tmp1 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 1);
0245         DetId tmp2 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 2);
0246         ESDetId strip1 = (tmp1 == DetId(0)) ? ESDetId(0) : ESDetId(tmp1);
0247         ESDetId strip2 = (tmp2 == DetId(0)) ? ESDetId(0) : ESDetId(tmp2);
0248
0249         if (nBC == 0) {
0250           if (strip1 != ESDetId(0) && strip2 != ESDetId(0)) {
0251             ESChannelStatusMap::const_iterator status_p1 = channelStatus->getMap().find(strip1);
0252             ESChannelStatusMap::const_iterator status_p2 = channelStatus->getMap().find(strip2);
0253             if (status_p1->getStatusCode() == 1)
0254               condP1 = 0;
0255             if (status_p2->getStatusCode() == 1)
0256               condP2 = 0;
0257           } else if (strip1 == ESDetId(0)) {
0258             condP1 = 0;
0259           } else if (strip2 == ESDetId(0)) {
0260             condP2 = 0;
0261           }
0262
0263           //cout<<"starting cluster : "<<maxDeltaPhi<<" "<<minDeltaPhi<<" "<<esPhiClusterDeltaEta_<<endl;
0264           //cout<<"do plane 1 === "<<strip1.zside()<<" "<<strip1.plane()<<" "<<strip1.six()<<" "<<strip1.siy()<<" "<<strip1.strip()<<endl;
0265           // Get a vector of ES clusters (found by the PreshSeeded algorithm) associated with a given EE basic cluster.
0266           reco::PreshowerCluster cl1 = presh_algo->makeOneCluster(
0267               strip1, &used_strips, &rechits_map, geometry_p, esPhiClusterDeltaEta_, minDeltaPhi, maxDeltaPhi);
0268           cl1.setBCRef(*bc_iter);
0269           clusters1.push_back(cl1);
0270           e1 += cl1.energy();
0271
0272           //cout<<"do plane 2 === "<<strip2.zside()<<" "<<strip2.plane()<<" "<<strip2.six()<<" "<<strip2.siy()<<" "<<strip2.strip()<<endl;
0273           reco::PreshowerCluster cl2 = presh_algo->makeOneCluster(
0274               strip2, &used_strips, &rechits_map, geometry_p, esPhiClusterDeltaEta_, minDeltaPhi, maxDeltaPhi);
0275           cl2.setBCRef(*bc_iter);
0276           clusters2.push_back(cl2);
0277           e2 += cl2.energy();
0278         }
0279       }
0280
0281       new_BC.push_back(*bc_iter);
0282       nBC++;
0283     }  // end of cycle over BCs
0284
0285     LogTrace("EcalClusters") << " For SC #" << isc - 1 << ", containing " << it_super->clustersSize()
0286                              << " basic clusters, PreshowerPhiClusterAlgo made " << clusters1.size()
0287                              << " in X plane and " << clusters2.size() << " in Y plane "
0288                              << " preshower clusters ";
0289
0290     // update energy of the SuperCluster
0291     if (e1 + e2 > 1.0e-10) {
0292       e1 = e1 / mip_;  // GeV to #MIPs
0293       e2 = e2 / mip_;
0294
0295       if (condP1 == 1 && condP2 == 1) {
0296         deltaE = gamma0_ * (e1 + alpha0_ * e2);
0297       } else if (condP1 == 1 && condP2 == 0) {
0298         deltaE = gamma1_ * (e1 + alpha1_ * e2);
0299       } else if (condP1 == 0 && condP2 == 1) {
0300         deltaE = gamma2_ * (e1 + alpha2_ * e2);
0301       } else if (condP1 == 0 && condP2 == 0) {
0302         deltaE = gamma3_ * (e1 + alpha3_ * e2);
0303       }
0304     }
0305
0306     //corrected Energy
0307     float E = it_super->energy() + deltaE;
0308
0309     LogTrace("EcalClusters") << " Creating corrected SC ";
0310
0311     reco::SuperCluster sc(E, it_super->position(), it_super->seed(), new_BC, deltaE);
0312     sc.setPreshowerEnergyPlane1(e1 * mip_);
0313     sc.setPreshowerEnergyPlane2(e2 * mip_);
0314     if (condP1 == 1 && condP2 == 1)
0315       sc.setPreshowerPlanesStatus(0);
0316     else if (condP1 == 1 && condP2 == 0)
0317       sc.setPreshowerPlanesStatus(1);
0318     else if (condP1 == 0 && condP2 == 1)
0319       sc.setPreshowerPlanesStatus(2);
0320     else if (condP1 == 0 && condP2 == 0)
0321       sc.setPreshowerPlanesStatus(3);
0322
0323     if (etThresh_ > 0) {  // calling postion().theta can be expensive
0324       if (sc.energy() * sin(sc.position().theta()) > etThresh_)
0325         new_SC.push_back(sc);
0326     } else {
0327       new_SC.push_back(sc);
0328     }
0329
0330   }  // end of cycle over SCs
0331
0332   // copy the preshower clusters into collections and put in the Event:
0333   clusters_p1->assign(clusters1.begin(), clusters1.end());
0334   clusters_p2->assign(clusters2.begin(), clusters2.end());
0335   // put collection of preshower clusters to the event
0336   evt.put(std::move(clusters_p1), preshClusterCollectionX_);
0337   evt.put(std::move(clusters_p2), preshClusterCollectionY_);
0338   LogTrace("EcalClusters") << "Preshower clusters added to the event";
0339
0340   // put collection of corrected super clusters to the event
0341   superclusters_p->assign(new_SC.begin(), new_SC.end());
0342   evt.put(std::move(superclusters_p), assocSClusterCollection_);
0343   LogTrace("EcalClusters") << "Corrected SClusters added to the event";
0344 }
0345
0346 void PreshowerPhiClusterProducer::set(const edm::EventSetup& es) {
0347   esgain_ = es.getHandle(esGainToken_);
0348   const ESGain* gain = esgain_.product();
0349
0350   double ESGain = gain->getESGain();
0351
0352   esMIPToGeV_ = es.getHandle(esMIPToGeVToken_);
0353   const ESMIPToGeVConstant* mipToGeV = esMIPToGeV_.product();
0354
0355   mip_ = (ESGain == 1) ? mipToGeV->getESValueLow() : mipToGeV->getESValueHigh();
0356
0357   esChannelStatus_ = es.getHandle(esChannelStatusToken_);
0358
0359   esEEInterCalib_ = es.getHandle(esEEInterCalibToken_);
0360   const ESEEIntercalibConstants* esEEInterCalib = esEEInterCalib_.product();
0361
0362   // both planes work
0363   gamma0_ = (ESGain == 1) ? 0.02 : esEEInterCalib->getGammaHigh0();
0364   alpha0_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow0() : esEEInterCalib->getAlphaHigh0();
0365
0366   // only first plane works
0367   gamma1_ = (ESGain == 1) ? (0.02 * esEEInterCalib->getGammaLow1()) : esEEInterCalib->getGammaHigh1();
0368   alpha1_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow1() : esEEInterCalib->getAlphaHigh1();
0369
0370   // only second plane works
0371   gamma2_ = (ESGain == 1) ? (0.02 * esEEInterCalib->getGammaLow2()) : esEEInterCalib->getGammaHigh2();
0372   alpha2_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow2() : esEEInterCalib->getAlphaHigh2();
0373
0374   // both planes do not work
0375   gamma3_ = (ESGain == 1) ? 0.02 : esEEInterCalib->getGammaHigh3();
0376   alpha3_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow3() : esEEInterCalib->getAlphaHigh3();
0377
0378   esMissingECalib_ = es.getHandle(esMissingECalibToken_);
0379   const ESMissingEnergyCalibration* esMissingECalib = esMissingECalib_.product();
0380
0381   // |eta| < 1.9
0382   aEta_[0] = esMissingECalib->getConstAEta0();
0383   bEta_[0] = esMissingECalib->getConstBEta0();
0384
0385   // 1.9 < |eta| < 2.1
0386   aEta_[1] = esMissingECalib->getConstAEta1();
0387   bEta_[1] = esMissingECalib->getConstBEta1();
0388
0389   // 2.1 < |eta| < 2.3
0390   aEta_[2] = esMissingECalib->getConstAEta2();
0391   bEta_[2] = esMissingECalib->getConstBEta2();
0392
0393   // 2.3 < |eta| < 2.5
0394   aEta_[3] = esMissingECalib->getConstAEta3();
0395   bEta_[3] = esMissingECalib->getConstBEta3();
0396 }