EgammaClusterProducers/src/PreshowerClusterShapeProducer.cc

0001 // authors A. Kyriakis, D. Maletic
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_;  // internal counter of events
0044
0045   //clustering parameters:
0046
0047   edm::EDGetTokenT<EcalRecHitCollection> preshHitToken_;                // name of module/plugin/producer
0048                                                                         // producing hits
0049   edm::EDGetTokenT<reco::SuperClusterCollection> endcapSClusterToken_;  // likewise for producer
0050                                                                         // of endcap superclusters
0051   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0052
0053   std::string PreshowerClusterShapeCollectionX_;
0054   std::string PreshowerClusterShapeCollectionY_;
0055
0056   EndcapPiZeroDiscriminatorAlgo* presh_pi0_algo;  // algorithm doing the real work
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   // use configuration file to setup input/output collection names
0069   // Parameters to identify the hit collections
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;  // use this stream for all messages in produce
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   // get the ECAL -> Preshower geometry and topology:
0107   ESHandle<CaloGeometry> geoHandle = es.getHandle(caloGeometryToken_);
0108   const CaloSubdetectorGeometry* geometry = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0109   const CaloSubdetectorGeometry*& geometry_p = geometry;
0110
0111   // create a unique_ptr to a PreshowerClusterShapeCollection
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   // fetch the Preshower product (RecHits)
0120   evt.getByToken(preshHitToken_, pRecHits);
0121   // pointer to the object in the product
0122   const EcalRecHitCollection* rechits = pRecHits.product();
0123
0124   LogTrace("EcalClusters") << "PreshowerClusterShapeProducer: ### Total # of preshower RecHits: " << rechits->size();
0125
0126   //  if ( rechits->size() <= 0 ) return;
0127
0128   // make the map of Preshower rechits:
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   //make cycle over Photon Collection
0141   int SC_index = 0;
0142   //  Handle<PhotonCollection> correctedPhotonHandle;
0143   //  evt.getByLabel(photonCorrCollectionProducer_, correctedPhotonCollection_ , correctedPhotonHandle);
0144   //  const PhotonCollection corrPhoCollection = *(correctedPhotonHandle.product());
0145   //  cout << " Photon Collection size : " << corrPhoCollection.size() << endl;
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) {  //  Use Preshower region only
0162       if (geometry) {
0163         const GlobalPoint pointSC(it_super->x(), it_super->y(), it_super->z());  // get the centroid of the SC
0164         LogTrace("EcalClusters") << "SC centroind = " << pointSC;
0165
0166         // Get the Preshower 2-planes RecHit vectors associated with the given SC
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     }  // end of cycle over Endcap SC
0195   }
0196   // put collection of PreshowerClusterShape in the Event:
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 }