CaloTowersCreator/src/CaloTowersMerger.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CaloTowersMerger
0004 // Class:      CaloTowersMerger
0005 //
0006 /**\class CaloTowersMerger CaloTowersMerger.cc RecoLocalCalo/CaloTowersMerger/src/CaloTowersMerger.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Jean-Roch Vlimant,40 3-A28,+41227671209,
0015 //         Created:  Thu Nov  4 16:36:30 CET 2010
0016 //
0017 //
0018
0019 // Anton Anastassov (Northwestern):
0020 // Add code for actual tower merging, define two distinct inputs of
0021 // "regular" and "extra" towers
0022 // This functionality is subject to some restrictions
0023 // (see notes below)
0024
0025 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0026 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/Framework/interface/stream/EDProducer.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031
0032 #include <memory>
0033
0034 //
0035 // class declaration
0036 //
0037
0038 class CaloTowersMerger : public edm::stream::EDProducer<> {
0039 public:
0040   explicit CaloTowersMerger(const edm::ParameterSet&);
0041   ~CaloTowersMerger() override;
0042
0043   CaloTower mergedTower(const CaloTower& t1, const CaloTower& t2);
0044
0045 private:
0046   void produce(edm::Event&, const edm::EventSetup&) override;
0047
0048   // ----------member data ---------------------------
0049
0050   edm::InputTag regularTowerTag, extraTowerTag;
0051   edm::EDGetTokenT<CaloTowerCollection> tok_reg_;
0052   edm::EDGetTokenT<CaloTowerCollection> tok_ext_;
0053 };
0054
0055 //
0056 // constants, enums and typedefs
0057 //
0058
0059 //
0060 // static data member definitions
0061 //
0062
0063 //
0064 // constructors and destructor
0065 //
0066 CaloTowersMerger::CaloTowersMerger(const edm::ParameterSet& iConfig) {
0067   regularTowerTag = iConfig.getParameter<edm::InputTag>("regularTowerTag");
0068   extraTowerTag = iConfig.getParameter<edm::InputTag>("extraTowerTag");
0069
0070   // register for data access
0071   tok_reg_ = consumes<CaloTowerCollection>(regularTowerTag);
0072   tok_ext_ = consumes<CaloTowerCollection>(extraTowerTag);
0073
0074   //register your products
0075   produces<CaloTowerCollection>();
0076 }
0077
0078 CaloTowersMerger::~CaloTowersMerger() {
0079   // do anything here that needs to be done at desctruction time
0080   // (e.g. close files, deallocate resources etc.)
0081 }
0082
0083 //
0084 // member functions
0085 //
0086
0087 // ------------ method called to produce the data  ------------
0088 void CaloTowersMerger::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0089   edm::Handle<CaloTowerCollection> regTower, extraTower;
0090
0091   iEvent.getByToken(tok_reg_, regTower);
0092   iEvent.getByToken(tok_ext_, extraTower);
0093
0094   if (!regTower.isValid() && !extraTower.isValid()) {
0095     edm::LogError("CaloTowersMerger") << "both input tag:" << regularTowerTag << " and " << extraTowerTag
0096                                       << " are invalid. empty merged collection";
0097     iEvent.put(std::make_unique<CaloTowerCollection>());
0098     return;
0099   } else if (!regTower.isValid() || !extraTower.isValid()) {
0100     if (!regTower.isValid() && extraTower.isValid())
0101       regTower = extraTower;
0102     iEvent.put(std::make_unique<CaloTowerCollection>(*regTower));
0103     return;
0104   } else {
0105     //both valid input collections: merging
0106     auto output = std::make_unique<CaloTowerCollection>();
0107     output->reserve(regTower->size() + extraTower->size());
0108
0109     CaloTowerCollection::const_iterator rt_begin = regTower->begin();
0110     CaloTowerCollection::const_iterator rt_end = regTower->end();
0111     CaloTowerCollection::const_iterator rt_it = rt_begin;
0112
0113     //vector of overlapping towers
0114     std::vector<CaloTowerCollection::const_iterator> overlappingTowers;
0115     overlappingTowers.reserve(extraTower->size());
0116
0117     for (; rt_it != rt_end; ++rt_it) {
0118       CaloTowerCollection::const_iterator et_it = extraTower->find(rt_it->id());
0119       if (et_it != extraTower->end()) {
0120         //need to merge the components
0121         //FIXME
0122
0123         /////   CaloTower mergedTower(*t1);
0124         //one needs to merge t1 and t2 into mergedTower
0125         //end FIXME
0126         CaloTower mt = mergedTower(*rt_it, *et_it);
0127
0128         output->push_back(mt);
0129         overlappingTowers.push_back(et_it);
0130
0131       } else {
0132         //just copy the regular tower over
0133         output->push_back(*rt_it);
0134       }
0135     }
0136     CaloTowerCollection::const_iterator et_begin = extraTower->begin();
0137     CaloTowerCollection::const_iterator et_end = extraTower->end();
0138     CaloTowerCollection::const_iterator et_it = et_begin;
0139     for (; et_it != et_end; ++et_it) {
0140       if (std::find(overlappingTowers.begin(), overlappingTowers.end(), et_it) == overlappingTowers.end())
0141         //non overlapping tower
0142         //copy the extra tower over
0143         output->push_back(*et_it);
0144     }
0145     iEvent.put(std::move(output));
0146   }
0147 }
0148
0149 // Make a new tower by merging two towers containing exclusive hits
0150 // This cannot be done fully consistently and independent of the
0151 // creation mechnaism...
0152 // This functionlaity it to be used only for testing the effects
0153 // of rejected bad hits.
0154
0155 CaloTower CaloTowersMerger::mergedTower(const CaloTower& rt, const CaloTower& et) {
0156   double newOuterE = 0;
0157
0158   // HO energies are always saved (even if useHO=false)
0159   // make sure there is no double counting
0160   // crude test if one has HO energy and the other not (possibly due to bad hit cleanup)
0161   // However: for |iEta|>15 E_outer has different meening:
0162   // it holds either the energy in the outer depths in HCAL, etc
0163
0164   if (rt.ietaAbs() < 16 && (fabs(rt.outerEnergy() - et.outerEnergy()) < 0.00001)) {
0165     // there is no differnce in the store HO energies
0166     newOuterE = rt.outerEnergy();
0167   } else {
0168     newOuterE = rt.outerEnergy() + et.outerEnergy();
0169   }
0170
0171   bool rt_hasEcalConstit = false;
0172   bool et_hasEcalConstit = false;
0173
0174   bool rt_hasHcalConstit = false;
0175   bool et_hasHcalConstit = false;
0176
0177   // check if there are HCAL/ECAL constituents in the towers
0178
0179   std::vector<DetId>::const_iterator rc_begin = rt.constituents().begin();
0180   std::vector<DetId>::const_iterator rc_end = rt.constituents().end();
0181   std::vector<DetId>::const_iterator rc_it;
0182
0183   for (rc_it = rc_begin; rc_it != rc_end; ++rc_it) {
0184     if (rc_it->det() == DetId::Hcal)
0185       rt_hasHcalConstit = true;
0186     break;
0187   }
0188   for (rc_it = rc_begin; rc_it != rc_end; ++rc_it) {
0189     if (rc_it->det() == DetId::Ecal)
0190       rt_hasEcalConstit = true;
0191     break;
0192   }
0193
0194   std::vector<DetId>::const_iterator ec_begin = et.constituents().begin();
0195   std::vector<DetId>::const_iterator ec_end = et.constituents().end();
0196   std::vector<DetId>::const_iterator ec_it;
0197
0198   for (ec_it = ec_begin; ec_it != ec_end; ++ec_it) {
0199     if (ec_it->det() == DetId::Hcal)
0200       et_hasHcalConstit = true;
0201     break;
0202   }
0203   for (ec_it = ec_begin; ec_it != ec_end; ++ec_it) {
0204     if (ec_it->det() == DetId::Ecal)
0205       et_hasEcalConstit = true;
0206     break;
0207   }
0208
0209   std::vector<DetId> combinedConstituents = rt.constituents();
0210   for (ec_it = ec_begin; ec_it != ec_end; ++ec_it) {
0211     // if properly resconstructed, the only possible overlap should be for HO hits that
0212     // are always listed as constituents if above thereshold (old JetMET request)
0213     if (std::find(combinedConstituents.begin(), combinedConstituents.end(), *ec_it) == combinedConstituents.end())
0214       combinedConstituents.push_back(*ec_it);
0215   }
0216
0217   GlobalPoint newEmPosition(0.0, 0.0, 0.0);
0218
0219   // The following assumes the current default
0220   // momentum reconstruction method (1) and
0221   // accepted rechits with some threshod >0
0222
0223   if (rt_hasEcalConstit && et_hasEcalConstit) {
0224     if (rt.emEnergy() > 0 && et.emEnergy() > 0) {
0225       double sumEmE = rt.emEnergy() + et.emEnergy();
0226
0227       double x = rt.emEnergy() * rt.emPosition().x() + et.emEnergy() * et.emPosition().x();
0228       double y = rt.emEnergy() * rt.emPosition().y() + et.emEnergy() * et.emPosition().y();
0229       double z = rt.emEnergy() * rt.emPosition().z() + et.emEnergy() * et.emPosition().z();
0230
0231       GlobalPoint weightedEmdPosition(x / sumEmE, y / sumEmE, z / sumEmE);
0232       newEmPosition = weightedEmdPosition;
0233     }
0234
0235   } else if (rt_hasEcalConstit && !et_hasEcalConstit) {
0236     newEmPosition = rt.emPosition();
0237   } else if (!rt_hasEcalConstit && et_hasEcalConstit) {
0238     newEmPosition = et.emPosition();
0239   }
0240
0241   GlobalPoint newHadPosition(0.0, 0.0, 0.0);
0242   // had positions are the same if there is at least one constituent
0243   if (rt_hasHcalConstit) {
0244     newHadPosition = rt.hadPosition();
0245   } else if (et_hasHcalConstit) {
0246     newHadPosition = et.hadPosition();
0247   }
0248
0249   // MAke the new tower and set all values
0250
0251   CaloTower mergedTower(rt.id(),
0252                         rt.emEnergy() + et.emEnergy(),
0253                         rt.hadEnergy() + et.hadEnergy(),
0254                         newOuterE,
0255                         -1,
0256                         -1,
0257                         rt.p4() + et.p4(),
0258                         newEmPosition,
0259                         newHadPosition);
0260
0261   mergedTower.addConstituents(combinedConstituents);
0262
0263   (rt.hottestCellE() > et.hottestCellE()) ? mergedTower.setHottestCellE(rt.hottestCellE())
0264                                           : mergedTower.setHottestCellE(et.hottestCellE());
0265
0266   unsigned int numBadHcalChan = rt.numBadHcalCells() - et.numProblematicHcalCells() - rt.numRecoveredHcalCells();
0267   unsigned int numBadEcalChan = rt.numBadEcalCells() - et.numProblematicEcalCells() - rt.numRecoveredEcalCells();
0268
0269   unsigned int numProbHcalChan = rt.numProblematicHcalCells() + et.numProblematicHcalCells();
0270   unsigned int numProbEcalChan = rt.numProblematicEcalCells() + et.numProblematicEcalCells();
0271
0272   unsigned int numRecHcalChan = rt.numRecoveredHcalCells() + et.numRecoveredHcalCells();
0273   unsigned int numRecEcalChan = rt.numRecoveredEcalCells() + et.numRecoveredEcalCells();
0274
0275   mergedTower.setCaloTowerStatus(
0276       numBadHcalChan, numBadEcalChan, numRecHcalChan, numRecEcalChan, numProbHcalChan, numProbEcalChan);
0277
0278   // use timing from the good tower only (for now, in default reco we use information from good hits only)
0279   // time is saved as integer but returned as float in (ns)
0280   mergedTower.setEcalTime(int(rt.ecalTime() * 100.0 + 0.5));
0281   mergedTower.setHcalTime(int(rt.hcalTime() * 100.0 + 0.5));
0282
0283   return mergedTower;
0284 }
0285
0286 //define this as a plug-in
0287 DEFINE_FWK_MODULE(CaloTowersMerger);