File indexing completed on 2024-09-07 04:37:44
0001
0002
0003
0004
0005
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
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
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
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] ";
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
0121
0122 edm::DetSetVector<StripDigiSimLink>::const_iterator isearch = stripdigisimlink_->find(detID);
0123 if (isearch == stripdigisimlink_->end())
0124 continue;
0125
0126
0127 edmNew::DetSetVector<SiStripCluster>::FastFiller outFill(output, det->id());
0128
0129
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;
0135 }
0136
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
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
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 }
0169
0170 size_t NsimTrk = trackID.size();
0171 if (NsimTrk < 2) {
0172 if (NsimTrk == 1)
0173 outFill.push_back(*clust);
0174
0175
0176 continue;
0177 }
0178
0179
0180
0181
0182
0183
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
0218
0219
0220
0221 outFill.push_back(SiStripCluster((uint16_t)TKfirstStrip[TKidx], TKampl[TKidx].begin(), TKampl[TKidx].end()));
0222 }
0223 }
0224 }
0225 }
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
0237 DEFINE_FWK_MODULE(ClusterMCsplitStrips);