Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:17:22

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 }