SiStripClusterizer/test/ClusterRefinerTagMCmerged.cc

0001 //
0002 // Use MC truth to identify merged clusters, i.e., those associated with more than one
0003 // (in-time) SimTrack.
0004 //
0005 // Author:  Bill Ford (wtford)  6 March 2015
0006 //
0007
0008 // system includes
0009 #include <memory>
0010
0011 // user includes
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Utilities/interface/InputTag.h"
0016 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0017 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019
0020 class ClusterRefinerTagMCmerged : public edm::stream::EDProducer<> {
0021 public:
0022   explicit ClusterRefinerTagMCmerged(const edm::ParameterSet& conf);
0023   virtual void produce(edm::Event&, const edm::EventSetup&);
0024
0025 private:
0026   template <class T>
0027   bool findInput(const edm::EDGetTokenT<T>&, edm::Handle<T>&, const edm::Event&);
0028   virtual void refineCluster(const edm::Handle<edmNew::DetSetVector<SiStripCluster>>& input,
0029                              edmNew::DetSetVector<SiStripCluster>& output);
0030
0031   const edm::InputTag inputTag;
0032   typedef edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster>> token_t;
0033   token_t inputToken;
0034   edm::ParameterSet confClusterRefiner_;
0035   bool useAssociateHit_;
0036
0037   TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0038   std::unique_ptr<TrackerHitAssociator> associator_;
0039 };
0040
0041 ClusterRefinerTagMCmerged::ClusterRefinerTagMCmerged(const edm::ParameterSet& conf)
0042     : inputTag(conf.getParameter<edm::InputTag>("UntaggedClusterProducer")),
0043       confClusterRefiner_(conf.getParameter<edm::ParameterSet>("ClusterRefiner")),
0044       useAssociateHit_(!confClusterRefiner_.getParameter<bool>("associateRecoTracks")),
0045       trackerHitAssociatorConfig_(confClusterRefiner_, consumesCollector()) {
0046   produces<edmNew::DetSetVector<SiStripCluster>>();
0047   inputToken = consumes<edmNew::DetSetVector<SiStripCluster>>(inputTag);
0048 }
0049
0050 void ClusterRefinerTagMCmerged::produce(edm::Event& event, const edm::EventSetup& es) {
0051   auto output = std::make_unique<edmNew::DetSetVector<SiStripCluster>>();
0052   output->reserve(10000, 4 * 10000);
0053
0054   associator_.reset(new TrackerHitAssociator(event, trackerHitAssociatorConfig_));
0055   edm::Handle<edmNew::DetSetVector<SiStripCluster>> input;
0056
0057   if (findInput(inputToken, input, event))
0058     refineCluster(input, *output);
0059   else
0060     edm::LogError("Input Not Found") << "[ClusterRefinerTagMCmerged::produce] ";  // << inputTag;
0061
0062   LogDebug("Output") << output->dataSize() << " clusters from " << output->size() << " modules";
0063   output->shrink_to_fit();
0064   event.put(std::move(output));
0065 }
0066
0067 void ClusterRefinerTagMCmerged::refineCluster(const edm::Handle<edmNew::DetSetVector<SiStripCluster>>& input,
0068                                               edmNew::DetSetVector<SiStripCluster>& output) {
0069   for (edmNew::DetSetVector<SiStripCluster>::const_iterator det = input->begin(); det != input->end(); det++) {
0070     // DetSetVector filler to receive the clusters we produce
0071     edmNew::DetSetVector<SiStripCluster>::TSFastFiller outFill(output, det->id());
0072     uint32_t detId = det->id();
0073     int ntk = 0;
0074     int NtkAll = 0;
0075     for (edmNew::DetSet<SiStripCluster>::iterator clust = det->begin(); clust != det->end(); clust++) {
0076       auto const& amp = clust->amplitudes();
0077       SiStripCluster* newCluster = new SiStripCluster(clust->firstStrip(), amp.begin(), amp.end());
0078       if (associator_ != 0) {
0079         std::vector<SimHitIdpr> simtrackid;
0080         if (useAssociateHit_) {
0081           std::vector<PSimHit> simhit;
0082           associator_->associateCluster(clust, DetId(detId), simtrackid, simhit);
0083           NtkAll = simtrackid.size();
0084           ntk = 0;
0085           if (simtrackid.size() > 1) {
0086             for (auto const& it : simtrackid) {
0087               int NintimeHit = 0;
0088               for (auto const& ih : simhit) {
0089                 // std::cout << "  hit(tk, evt) trk(tk, evt) bunch ("
0090                 //    << ih.trackId() << ", " << ih.eventId().rawId() << ") ("
0091                 //    << it.first << ", " << it.second.rawId() << ") "
0092                 //    << ih.eventId().bunchCrossing()
0093                 //    << std::endl;
0094                 if (ih.trackId() == it.first && ih.eventId() == it.second && ih.eventId().bunchCrossing() == 0)
0095                   ++NintimeHit;
0096               }
0097               if (NintimeHit > 0)
0098                 ++ntk;
0099             }
0100           }
0101         } else {
0102           associator_->associateSimpleRecHitCluster(clust, DetId(detId), simtrackid);
0103           ntk = NtkAll = simtrackid.size();
0104         }
0105         if (ntk > 1) {
0106           newCluster->setMerged(true);
0107         } else {
0108           newCluster->setMerged(false);
0109         }
0110       }
0111       outFill.push_back(*newCluster);
0112       // std::cout << "t_m_w " << " " << NtkAll << " " << newCluster->isMerged()
0113       //        << " " << clust->amplitudes().size()
0114       //        << std::endl;
0115     }  // traverse clusters
0116   }  // traverse sensors
0117 }
0118
0119 template <class T>
0120 inline bool ClusterRefinerTagMCmerged::findInput(const edm::EDGetTokenT<T>& tag,
0121                                                  edm::Handle<T>& handle,
0122                                                  const edm::Event& e) {
0123   e.getByToken(tag, handle);
0124   return handle.isValid();
0125 }
0126
0127 #include "FWCore/Framework/interface/MakerMacros.h"
0128 DEFINE_FWK_MODULE(ClusterRefinerTagMCmerged);