EcalTestBeam/plugins/TBHodoActiveVolumeRawInfoProducer.cc

0001 /*
0002  * \file TBHodoActiveVolumeRawInfoProducer.cc
0003  *
0004  *
0005  */
0006
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/stream/EDProducer.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/PluginManager/interface/ModuleDef.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015
0016 #include "Geometry/EcalTestBeam/interface/EcalTBHodoscopeGeometry.h"
0017 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0018 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0019 #include "SimDataFormats/EcalTestBeam/interface/HodoscopeDetId.h"
0020 #include "SimDataFormats/EcalTestBeam/interface/PEcalTBInfo.h"
0021 #include "TBDataFormats/EcalTBObjects/interface/EcalTBHodoscopePlaneRawHits.h"
0022 #include "TBDataFormats/EcalTBObjects/interface/EcalTBHodoscopeRawInfo.h"
0023
0024 #include <iostream>
0025 #include <memory>
0026
0027 class TBHodoActiveVolumeRawInfoProducer : public edm::stream::EDProducer<> {
0028 public:
0029   /// Constructor
0030   explicit TBHodoActiveVolumeRawInfoProducer(const edm::ParameterSet &ps);
0031
0032   /// Destructor
0033   ~TBHodoActiveVolumeRawInfoProducer() override = default;
0034
0035   /// Produce digis out of raw data
0036   void produce(edm::Event &event, const edm::EventSetup &eventSetup) override;
0037
0038 private:
0039   double myThreshold;
0040   edm::EDGetTokenT<edm::PCaloHitContainer> m_EcalToken;
0041
0042   std::unique_ptr<EcalTBHodoscopeGeometry> theTBHodoGeom_;
0043 };
0044
0045 using namespace cms;
0046 using namespace std;
0047
0048 TBHodoActiveVolumeRawInfoProducer::TBHodoActiveVolumeRawInfoProducer(const edm::ParameterSet &ps)
0049     : myThreshold(0.05E-3),
0050       m_EcalToken(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalTBH4BeamHits"))) {
0051   produces<EcalTBHodoscopeRawInfo>();
0052
0053   theTBHodoGeom_ = std::make_unique<EcalTBHodoscopeGeometry>();
0054 }
0055
0056 void TBHodoActiveVolumeRawInfoProducer::produce(edm::Event &event, const edm::EventSetup &eventSetup) {
0057   std::unique_ptr<EcalTBHodoscopeRawInfo> product(new EcalTBHodoscopeRawInfo());
0058
0059   // caloHit container
0060   const edm::Handle<edm::PCaloHitContainer> &pCaloHit = event.getHandle(m_EcalToken);
0061   const edm::PCaloHitContainer *caloHits = nullptr;
0062   if (pCaloHit.isValid()) {
0063     caloHits = pCaloHit.product();
0064     LogDebug("EcalTBHodo") << "total # caloHits: " << caloHits->size();
0065   } else {
0066     edm::LogError("EcalTBHodo") << "Error! can't get the caloHitContainer ";
0067   }
0068   if (!caloHits) {
0069     return;
0070   }
0071
0072   // detid - energy_sum map
0073   std::map<unsigned int, double> energyMap;
0074
0075   for (auto &&aHit : *caloHits) {
0076     double thisHitEne = aHit.energy();
0077
0078     std::map<unsigned int, double>::iterator itmap = energyMap.find(aHit.id());
0079     if (itmap == energyMap.end())
0080       energyMap.insert(pair<unsigned int, double>(aHit.id(), thisHitEne));
0081     else {
0082       (*itmap).second += thisHitEne;
0083     }
0084   }
0085
0086   // planes and fibers
0087   int nPlanes = theTBHodoGeom_->getNPlanes();
0088   int nFibers = theTBHodoGeom_->getNFibres();
0089   product->setPlanes(nPlanes);
0090
0091   bool firedChannels[4][64];
0092   for (int iPlane = 0; iPlane < nPlanes; ++iPlane) {
0093     for (int iFiber = 0; iFiber < nFibers; ++iFiber) {
0094       firedChannels[iPlane][iFiber] = 0.;
0095     }
0096   }
0097   for (std::map<unsigned int, double>::const_iterator itmap = energyMap.begin(); itmap != energyMap.end(); ++itmap) {
0098     if ((*itmap).second > myThreshold) {
0099       HodoscopeDetId myHodoDetId = HodoscopeDetId((*itmap).first);
0100       firedChannels[myHodoDetId.planeId()][myHodoDetId.fibrId()] = true;
0101     }
0102   }
0103   for (int iPlane = 0; iPlane < nPlanes; ++iPlane) {
0104     EcalTBHodoscopePlaneRawHits planeHit(nFibers);
0105
0106     for (int iFiber = 0; iFiber < nFibers; ++iFiber) {
0107       planeHit.setHit(iFiber, firedChannels[iPlane][iFiber]);
0108     }
0109     product->setPlane((unsigned int)iPlane, planeHit);
0110   }
0111
0112   LogDebug("EcalTBHodo") << (*product);
0113
0114   event.put(std::move(product));
0115 }
0116
0117 DEFINE_FWK_MODULE(TBHodoActiveVolumeRawInfoProducer);