File indexing completed on 2023-03-17 11:17:22
0001
0002
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0005 #include "DataFormats/EgammaReco/interface/PreshowerClusterShape.h"
0006 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0007 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0008 #include "DataFormats/Math/interface/Point3D.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/ESGetToken.h"
0017 #include "FWCore/Utilities/interface/Exception.h"
0018 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0019 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0020 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0021 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0022 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0023 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "RecoEcal/EgammaClusterAlgos/interface/EndcapPiZeroDiscriminatorAlgo.h"
0026
0027 #include <fstream>
0028 #include <memory>
0029 #include <sstream>
0030 #include <vector>
0031
0032 class PreshowerClusterShapeProducer : public edm::stream::EDProducer<> {
0033 public:
0034 typedef math::XYZPoint Point;
0035
0036 explicit PreshowerClusterShapeProducer(const edm::ParameterSet& ps);
0037
0038 ~PreshowerClusterShapeProducer() override;
0039
0040 void produce(edm::Event& evt, const edm::EventSetup& es) override;
0041
0042 private:
0043 int nEvt_;
0044
0045
0046
0047 edm::EDGetTokenT<EcalRecHitCollection> preshHitToken_;
0048
0049 edm::EDGetTokenT<reco::SuperClusterCollection> endcapSClusterToken_;
0050
0051 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0052
0053 std::string PreshowerClusterShapeCollectionX_;
0054 std::string PreshowerClusterShapeCollectionY_;
0055
0056 EndcapPiZeroDiscriminatorAlgo* presh_pi0_algo;
0057 };
0058
0059 #include "FWCore/Framework/interface/MakerMacros.h"
0060 DEFINE_FWK_MODULE(PreshowerClusterShapeProducer);
0061
0062 using namespace std;
0063 using namespace reco;
0064 using namespace edm;
0065
0066
0067 PreshowerClusterShapeProducer::PreshowerClusterShapeProducer(const ParameterSet& ps) {
0068
0069
0070 preshHitToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("preshRecHitProducer"));
0071 endcapSClusterToken_ =
0072 consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("endcapSClusterProducer"));
0073 caloGeometryToken_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0074
0075 PreshowerClusterShapeCollectionX_ = ps.getParameter<string>("PreshowerClusterShapeCollectionX");
0076 PreshowerClusterShapeCollectionY_ = ps.getParameter<string>("PreshowerClusterShapeCollectionY");
0077
0078 produces<reco::PreshowerClusterShapeCollection>(PreshowerClusterShapeCollectionX_);
0079 produces<reco::PreshowerClusterShapeCollection>(PreshowerClusterShapeCollectionY_);
0080
0081 float preshStripECut = ps.getParameter<double>("preshStripEnergyCut");
0082 int preshNst = ps.getParameter<int>("preshPi0Nstrip");
0083
0084 string debugString = ps.getParameter<string>("debugLevel");
0085
0086 string tmpPath = ps.getUntrackedParameter<string>("pathToWeightFiles", "RecoEcal/EgammaClusterProducers/data/");
0087
0088 presh_pi0_algo = new EndcapPiZeroDiscriminatorAlgo(preshStripECut, preshNst, tmpPath);
0089
0090 LogTrace("EcalClusters") << "PreshowerClusterShapeProducer:presh_pi0_algo class instantiated ";
0091
0092 nEvt_ = 0;
0093 }
0094
0095 PreshowerClusterShapeProducer::~PreshowerClusterShapeProducer() { delete presh_pi0_algo; }
0096
0097 void PreshowerClusterShapeProducer::produce(Event& evt, const EventSetup& es) {
0098 ostringstream ostr;
0099
0100 LogTrace("EcalClusters") << "\n ....... Event " << evt.id() << " with Number = " << nEvt_ + 1
0101 << " is analyzing ....... ";
0102
0103 Handle<EcalRecHitCollection> pRecHits;
0104 Handle<SuperClusterCollection> pSuperClusters;
0105
0106
0107 ESHandle<CaloGeometry> geoHandle = es.getHandle(caloGeometryToken_);
0108 const CaloSubdetectorGeometry* geometry = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0109 const CaloSubdetectorGeometry*& geometry_p = geometry;
0110
0111
0112 auto ps_cl_for_pi0_disc_x = std::make_unique<reco::PreshowerClusterShapeCollection>();
0113 auto ps_cl_for_pi0_disc_y = std::make_unique<reco::PreshowerClusterShapeCollection>();
0114
0115 std::unique_ptr<CaloSubdetectorTopology> topology_p;
0116 if (geometry)
0117 topology_p = std::make_unique<EcalPreshowerTopology>();
0118
0119
0120 evt.getByToken(preshHitToken_, pRecHits);
0121
0122 const EcalRecHitCollection* rechits = pRecHits.product();
0123
0124 LogTrace("EcalClusters") << "PreshowerClusterShapeProducer: ### Total # of preshower RecHits: " << rechits->size();
0125
0126
0127
0128
0129 map<DetId, EcalRecHit> rechits_map;
0130 EcalRecHitCollection::const_iterator it;
0131 for (it = rechits->begin(); it != rechits->end(); it++) {
0132 rechits_map.insert(make_pair(it->id(), *it));
0133 }
0134
0135 LogTrace("EcalClusters") << "PreshowerClusterShapeProducer: ### Preshower RecHits_map of size " << rechits_map.size()
0136 << " was created!";
0137
0138 reco::PreshowerClusterShapeCollection ps_cl_x, ps_cl_y;
0139
0140
0141 int SC_index = 0;
0142
0143
0144
0145
0146
0147 evt.getByToken(endcapSClusterToken_, pSuperClusters);
0148 const reco::SuperClusterCollection* SClusts = pSuperClusters.product();
0149 LogTrace("EcalClusters") << "### Total # Endcap Superclusters: " << SClusts->size();
0150
0151 SuperClusterCollection::const_iterator it_s;
0152 for (it_s = SClusts->begin(); it_s != SClusts->end(); it_s++) {
0153 SuperClusterRef it_super(reco::SuperClusterRef(pSuperClusters, SC_index));
0154
0155 float SC_eta = it_super->eta();
0156
0157 LogTrace("EcalClusters") << "PreshowerClusterShapeProducer: superCl_E = " << it_super->energy()
0158 << " superCl_Et = " << it_super->energy() * sin(2 * atan(exp(-it_super->eta())))
0159 << " superCl_Eta = " << SC_eta << " superCl_Phi = " << it_super->phi();
0160
0161 if (fabs(SC_eta) >= 1.65 && fabs(SC_eta) <= 2.5) {
0162 if (geometry) {
0163 const GlobalPoint pointSC(it_super->x(), it_super->y(), it_super->z());
0164 LogTrace("EcalClusters") << "SC centroind = " << pointSC;
0165
0166
0167
0168 DetId tmp_stripX = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(pointSC, 1);
0169 DetId tmp_stripY = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(pointSC, 2);
0170 ESDetId stripX = (tmp_stripX == DetId(0)) ? ESDetId(0) : ESDetId(tmp_stripX);
0171 ESDetId stripY = (tmp_stripY == DetId(0)) ? ESDetId(0) : ESDetId(tmp_stripY);
0172
0173 vector<float> vout_stripE1 = presh_pi0_algo->findPreshVector(stripX, &rechits_map, topology_p.get());
0174 vector<float> vout_stripE2 = presh_pi0_algo->findPreshVector(stripY, &rechits_map, topology_p.get());
0175
0176 LogTrace("EcalClusters") << "PreshowerClusterShapeProducer : ES Energy vector associated to the given SC = ";
0177 for (int k1 = 0; k1 < 11; k1++) {
0178 LogTrace("EcalClusters") << vout_stripE1[k1] << " ";
0179 }
0180
0181 for (int k1 = 0; k1 < 11; k1++) {
0182 LogTrace("EcalClusters") << vout_stripE2[k1] << " ";
0183 }
0184
0185 reco::PreshowerClusterShape ps1 = reco::PreshowerClusterShape(vout_stripE1, 1);
0186 ps1.setSCRef(it_super);
0187 ps_cl_x.push_back(ps1);
0188
0189 reco::PreshowerClusterShape ps2 = reco::PreshowerClusterShape(vout_stripE2, 2);
0190 ps2.setSCRef(it_super);
0191 ps_cl_y.push_back(ps2);
0192 }
0193 SC_index++;
0194 }
0195 }
0196
0197 ps_cl_for_pi0_disc_x->assign(ps_cl_x.begin(), ps_cl_x.end());
0198 ps_cl_for_pi0_disc_y->assign(ps_cl_y.begin(), ps_cl_y.end());
0199
0200 evt.put(std::move(ps_cl_for_pi0_disc_x), PreshowerClusterShapeCollectionX_);
0201 evt.put(std::move(ps_cl_for_pi0_disc_y), PreshowerClusterShapeCollectionY_);
0202 LogTrace("EcalClusters") << "PreshowerClusterShapeCollection added to the event";
0203
0204 nEvt_++;
0205
0206 LogDebug("PiZeroDiscriminatorDebug") << ostr.str();
0207 }