SubCollectionProducers/test/ClusterMCsplitStrips.cc

0001 //
0002 // Use MC truth to identify merged strip clusters, i.e., those associated with more
0003 // than one SimTrack, and split them into their true components.
0004 //
0005 // Author:  Bill Ford (wtford)  17 March 2015
0006 //
0007
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/stream/EDProducer.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Utilities/interface/InputTag.h"
0015 #include "DataFormats/Common/interface/DetSetVector.h"
0016 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0017 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0018 #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
0019 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0020 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0021
0022 #include <memory>
0023
0024 // Class declaration
0025
0026 class ClusterMCsplitStrips : public edm::stream::EDProducer<> {
0027 public:
0028   explicit ClusterMCsplitStrips(const edm::ParameterSet& conf);
0029   virtual void produce(edm::Event&, const edm::EventSetup&);
0030
0031 private:
0032   template <class T>
0033   bool findInput(const edm::EDGetTokenT<T>&, edm::Handle<T>&, const edm::Event&);
0034   void refineCluster(const edm::Handle<edmNew::DetSetVector<SiStripCluster>>& input,
0035                      edmNew::DetSetVector<SiStripCluster>& output);
0036
0037   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> esTopoToken;
0038   const edm::InputTag inputTag;
0039   typedef edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster>> token_t;
0040   token_t inputToken;
0041   edm::ParameterSet confClusterRefiner_;
0042   edm::EDGetTokenT<edm::DetSetVector<StripDigiSimLink>> stripdigisimlinkToken;
0043   edm::Handle<edm::DetSetVector<StripDigiSimLink>> stripdigisimlink_;
0044   edm::ESHandle<TrackerTopology> tTopoHandle_;
0045
0046   std::vector<std::string> moduleTypeStrings_;
0047   std::vector<SiStripModuleGeometry> moduleTypeCodes_;
0048 };
0049
0050 // Class implementation
0051
0052 ClusterMCsplitStrips::ClusterMCsplitStrips(const edm::ParameterSet& conf)
0053     : esTopoToken(esConsumes()),
0054       inputTag(conf.getParameter<edm::InputTag>("UnsplitClusterProducer")),
0055       confClusterRefiner_(conf.getParameter<edm::ParameterSet>("ClusterRefiner")) {
0056   produces<edmNew::DetSetVector<SiStripCluster>>();
0057   inputToken = consumes<edmNew::DetSetVector<SiStripCluster>>(inputTag);
0058   stripdigisimlinkToken = consumes<edm::DetSetVector<StripDigiSimLink>>(edm::InputTag("simSiStripDigis"));
0059   moduleTypeStrings_ = confClusterRefiner_.getParameter<std::vector<std::string>>("moduleTypes");
0060   for (auto mod = moduleTypeStrings_.begin(); mod != moduleTypeStrings_.end(); ++mod) {
0061     if (*mod == "IB1")
0062       moduleTypeCodes_.push_back(SiStripModuleGeometry::IB1);
0063     if (*mod == "IB2")
0064       moduleTypeCodes_.push_back(SiStripModuleGeometry::IB2);
0065     if (*mod == "OB1")
0066       moduleTypeCodes_.push_back(SiStripModuleGeometry::OB1);
0067     if (*mod == "OB2")
0068       moduleTypeCodes_.push_back(SiStripModuleGeometry::OB2);
0069     if (*mod == "W1A")
0070       moduleTypeCodes_.push_back(SiStripModuleGeometry::W1A);
0071     if (*mod == "W2A")
0072       moduleTypeCodes_.push_back(SiStripModuleGeometry::W2A);
0073     if (*mod == "W3A")
0074       moduleTypeCodes_.push_back(SiStripModuleGeometry::W3A);
0075     if (*mod == "W1B")
0076       moduleTypeCodes_.push_back(SiStripModuleGeometry::W1B);
0077     if (*mod == "W2B")
0078       moduleTypeCodes_.push_back(SiStripModuleGeometry::W2B);
0079     if (*mod == "W3B")
0080       moduleTypeCodes_.push_back(SiStripModuleGeometry::W3B);
0081     if (*mod == "W4")
0082       moduleTypeCodes_.push_back(SiStripModuleGeometry::W4);
0083     if (*mod == "W5")
0084       moduleTypeCodes_.push_back(SiStripModuleGeometry::W5);
0085     if (*mod == "W6")
0086       moduleTypeCodes_.push_back(SiStripModuleGeometry::W6);
0087     if (*mod == "W7")
0088       moduleTypeCodes_.push_back(SiStripModuleGeometry::W7);
0089   }
0090 }
0091
0092 void ClusterMCsplitStrips::produce(edm::Event& event, const edm::EventSetup& evSetup) {
0093   //Retrieve tracker topology from geometry
0094   tTopoHandle_ = evSetup.getHandle(esTopoToken);
0095
0096   auto output = std::make_unique<edmNew::DetSetVector<SiStripCluster>>();
0097   output->reserve(10000, 4 * 10000);
0098
0099   event.getByToken(stripdigisimlinkToken, stripdigisimlink_);
0100
0101   edm::Handle<edmNew::DetSetVector<SiStripCluster>> input;
0102
0103   if (findInput(inputToken, input, event))
0104     refineCluster(input, *output);
0105   else
0106     edm::LogError("Input Not Found") << "[ClusterMCsplitStrips::produce] ";  // << inputTag;
0107
0108   LogDebug("Output") << output->dataSize() << " clusters from " << output->size() << " modules";
0109   output->shrink_to_fit();
0110   event.put(std::move(output));
0111 }
0112
0113 void ClusterMCsplitStrips::refineCluster(const edm::Handle<edmNew::DetSetVector<SiStripCluster>>& input,
0114                                          edmNew::DetSetVector<SiStripCluster>& output) {
0115   const TrackerTopology* const tTopo = tTopoHandle_.product();
0116
0117   for (edmNew::DetSetVector<SiStripCluster>::const_iterator det = input->begin(); det != input->end(); det++) {
0118     uint32_t detID = det->id();
0119     DetId theDet(detID);
0120     // int subDet = theDet.subdetId();
0121
0122     edm::DetSetVector<StripDigiSimLink>::const_iterator isearch = stripdigisimlink_->find(detID);
0123     if (isearch == stripdigisimlink_->end())
0124       continue;  // This sensor has no simlinks;
0125                  // Any clusters must be from noise.
0126     // DetSetVector filler to receive the clusters we produce
0127     edmNew::DetSetVector<SiStripCluster>::FastFiller outFill(output, det->id());
0128
0129     // Consider clusters of selected module types for splitting; else just push original cluster to output.
0130     if (std::find(moduleTypeCodes_.begin(), moduleTypeCodes_.end(), tTopo->moduleGeometry(theDet)) ==
0131         moduleTypeCodes_.end()) {
0132       for (auto clust = det->begin(); clust != det->end(); clust++)
0133         outFill.push_back(*clust);
0134       continue;  // On to the next sensor.
0135     }
0136     // Traverse the clusters for this sensor.
0137     for (edmNew::DetSet<SiStripCluster>::iterator clust = det->begin(); clust != det->end(); clust++) {
0138       int clusiz = clust->amplitudes().size();
0139       int first = clust->firstStrip();
0140       int last = first + clusiz;
0141       edm::DetSet<StripDigiSimLink> link_detset = (*isearch);
0142
0143       //  First pass to count simTracks and set this cluster's range in the simlink vector
0144       edm::DetSet<StripDigiSimLink>::const_iterator firstlink = link_detset.data.begin(),
0145                                                     lastlink = link_detset.data.end();
0146       bool firstlinkInit = false;
0147
0148       std::vector<unsigned int> trackID;
0149       for (edm::DetSet<StripDigiSimLink>::const_iterator linkiter = link_detset.data.begin();
0150            linkiter != link_detset.data.end();
0151            linkiter++) {
0152         // DigiSimLinks are ordered first by channel; there can be > 1 track, and > 1 simHit for each track
0153         int thisChannel = linkiter->channel();
0154         if (thisChannel < first)
0155           continue;
0156         if (thisChannel >= first && !firstlinkInit) {
0157           firstlinkInit = true;
0158           firstlink = linkiter;
0159         }
0160         if (thisChannel >= last)
0161           break;
0162         lastlink = linkiter;
0163         lastlink++;
0164         auto const& thisSimTrk = linkiter->SimTrackId();
0165         if (std::find(trackID.begin(), trackID.end(), thisSimTrk) == trackID.end()) {
0166           trackID.push_back(thisSimTrk);
0167         }
0168       }  // end initial loop over this sensor's simlinks
0169
0170       size_t NsimTrk = trackID.size();
0171       if (NsimTrk < 2) {
0172         if (NsimTrk == 1)
0173           outFill.push_back(*clust);  // Unmerged cluster:  push it to the output.
0174
0175         //  (If NsimTrk = 0, cluster has no matched simlinks; abandon it.)
0176         continue;  // On to the next cluster
0177       }
0178
0179       // std::cout << "subDet " << DetId(detID).subdetId() << ", det " << detID << ": " << NsimTrk << " simTracks:";
0180       // for (unsigned int i=0; i<NsimTrk; i++) std::cout << " " << trackID[i];
0181       // std::cout << std::endl;
0182
0183       //  This cluster matches more than one simTrack, so we proceed to split it.
0184       auto const& amp = clust->amplitudes();
0185       std::vector<int> TKfirstStrip(NsimTrk, -1);
0186       std::vector<std::vector<uint16_t>> TKampl(NsimTrk);
0187       std::vector<int> prevStrip(NsimTrk, -1);
0188
0189       for (edm::DetSet<StripDigiSimLink>::const_iterator linkiter = firstlink; linkiter != lastlink; linkiter++) {
0190         int stripIdx = (int)linkiter->channel() - first;
0191
0192         uint16_t rawAmpl = (uint16_t)(amp[stripIdx]);
0193         uint16_t thisAmpl;
0194         if (rawAmpl < 254) {
0195           thisAmpl = std::min(uint16_t(253), std::max(uint16_t(0), (uint16_t)(rawAmpl * linkiter->fraction() + 0.5)));
0196         } else {
0197           thisAmpl = rawAmpl;
0198         }
0199
0200         unsigned int thisSimTrk = linkiter->SimTrackId();
0201         auto const& TKiter = std::find(trackID.begin(), trackID.end(), thisSimTrk);
0202         unsigned int TKidx = TKiter - trackID.begin();
0203
0204         if (TKfirstStrip[TKidx] == -1)
0205           TKfirstStrip[TKidx] = linkiter->channel();
0206         if (stripIdx != prevStrip[TKidx]) {
0207           prevStrip[TKidx] = stripIdx;
0208           TKampl[TKidx].push_back(thisAmpl);
0209         } else {
0210           if (rawAmpl < 254)
0211             (TKampl[TKidx])[linkiter->channel() - TKfirstStrip[TKidx]] += thisAmpl;
0212         }
0213       }
0214
0215       for (unsigned int TKidx = 0; TKidx < NsimTrk; ++TKidx) {
0216         if (std::accumulate(TKampl[TKidx].begin(), TKampl[TKidx].end(), 0) > 0) {
0217           // std::cout << "SimTrackID, 1st: ampl " << trackID[TKidx] << ", " << TKfirstStrip[TKidx] << ":";
0218           // for (unsigned i=0; i<TKampl[TKidx].size(); ++i) std::cout << " " << TKampl[TKidx][i];
0219           // std::cout << std::endl;
0220
0221           outFill.push_back(SiStripCluster((uint16_t)TKfirstStrip[TKidx], TKampl[TKidx].begin(), TKampl[TKidx].end()));
0222         }
0223       }
0224     }  // end loop over original clusters
0225   }    // end loop over sensors
0226 }
0227
0228 template <class T>
0229 inline bool ClusterMCsplitStrips::findInput(const edm::EDGetTokenT<T>& tag,
0230                                             edm::Handle<T>& handle,
0231                                             const edm::Event& e) {
0232   e.getByToken(tag, handle);
0233   return handle.isValid();
0234 }
0235
0236 //define this as a plug-in
0237 DEFINE_FWK_MODULE(ClusterMCsplitStrips);