File indexing completed on 2024-04-06 11:59:46
0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/ESHandle.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "DataFormats/SiStripDigi/interface/SiStripProcessedRawDigi.h"
0005 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0006 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0007 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0008 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0009 #include "FWCore/Framework/interface/stream/EDProducer.h"
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012 #include "FWCore/Utilities/interface/EDGetToken.h"
0013 #include "DataFormats/Common/interface/DetSetVector.h"
0014 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0015 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0016 #include "RecoLocalTracker/SiStripClusterizer/interface/SiStripClusterInfo.h"
0017
0018 class ShallowClustersProducer : public edm::stream::EDProducer<> {
0019 public:
0020 explicit ShallowClustersProducer(const edm::ParameterSet&);
0021
0022 private:
0023 edm::InputTag theClustersLabel;
0024 std::string Prefix;
0025 void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0026
0027 struct moduleVars {
0028 moduleVars(uint32_t, const TrackerTopology*);
0029 int subdetid, side, layerwheel, stringringrod, petal, stereo;
0030 uint32_t module;
0031 };
0032
0033 struct NearDigis {
0034 NearDigis(const SiStripClusterInfo&);
0035 NearDigis(const SiStripClusterInfo&, const edm::DetSetVector<SiStripProcessedRawDigi>&);
0036 float max, left, right, first, last, Lleft, Rright;
0037 float etaX() const { return ((left + right) / max) / 2.; }
0038 float eta() const { return right > left ? max / (max + right) : left / (left + max); }
0039 float etaasymm() const { return right > left ? (right - max) / (right + max) : (max - left) / (max + left); }
0040 float outsideasymm() const { return (last - first) / (last + first); }
0041 };
0042
0043 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0044 edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster>> theClustersToken_;
0045 edm::EDGetTokenT<edm::DetSetVector<SiStripProcessedRawDigi>> theDigisToken_;
0046 SiStripClusterInfo siStripClusterInfo_;
0047 };
0048
0049 #include "FWCore/Framework/interface/MakerMacros.h"
0050 DEFINE_FWK_MODULE(ShallowClustersProducer);
0051
0052 ShallowClustersProducer::ShallowClustersProducer(const edm::ParameterSet& iConfig)
0053 : Prefix(iConfig.getParameter<std::string>("Prefix")),
0054 topoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0055 siStripClusterInfo_(consumesCollector()) {
0056 produces<std::vector<unsigned>>(Prefix + "number");
0057 produces<std::vector<unsigned>>(Prefix + "width");
0058 produces<std::vector<float>>(Prefix + "variance");
0059 produces<std::vector<float>>(Prefix + "barystrip");
0060 produces<std::vector<float>>(Prefix + "middlestrip");
0061 produces<std::vector<unsigned>>(Prefix + "charge");
0062 produces<std::vector<float>>(Prefix + "noise");
0063 produces<std::vector<float>>(Prefix + "ston");
0064 produces<std::vector<unsigned>>(Prefix + "seedstrip");
0065 produces<std::vector<unsigned>>(Prefix + "seedindex");
0066 produces<std::vector<unsigned>>(Prefix + "seedcharge");
0067 produces<std::vector<float>>(Prefix + "seednoise");
0068 produces<std::vector<float>>(Prefix + "seedgain");
0069 produces<std::vector<unsigned>>(Prefix + "qualityisbad");
0070
0071 produces<std::vector<float>>(Prefix + "rawchargeC");
0072 produces<std::vector<float>>(Prefix + "rawchargeL");
0073 produces<std::vector<float>>(Prefix + "rawchargeR");
0074 produces<std::vector<float>>(Prefix + "rawchargeLL");
0075 produces<std::vector<float>>(Prefix + "rawchargeRR");
0076 produces<std::vector<float>>(Prefix + "eta");
0077 produces<std::vector<float>>(Prefix + "foldedeta");
0078 produces<std::vector<float>>(Prefix + "etaX");
0079 produces<std::vector<float>>(Prefix + "etaasymm");
0080 produces<std::vector<float>>(Prefix + "outsideasymm");
0081 produces<std::vector<float>>(Prefix + "neweta");
0082 produces<std::vector<float>>(Prefix + "newetaerr");
0083
0084 produces<std::vector<unsigned>>(Prefix + "detid");
0085 produces<std::vector<int>>(Prefix + "subdetid");
0086 produces<std::vector<int>>(Prefix + "module");
0087 produces<std::vector<int>>(Prefix + "side");
0088 produces<std::vector<int>>(Prefix + "layerwheel");
0089 produces<std::vector<int>>(Prefix + "stringringrod");
0090 produces<std::vector<int>>(Prefix + "petal");
0091 produces<std::vector<int>>(Prefix + "stereo");
0092
0093 theClustersToken_ = consumes<edmNew::DetSetVector<SiStripCluster>>(iConfig.getParameter<edm::InputTag>("Clusters"));
0094 theDigisToken_ = consumes<edm::DetSetVector<SiStripProcessedRawDigi>>(edm::InputTag("siStripProcessedRawDigis", ""));
0095 }
0096
0097 void ShallowClustersProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0098
0099 edm::ESHandle<TrackerTopology> tTopoHandle = iSetup.getHandle(topoToken_);
0100 const TrackerTopology* const tTopo = tTopoHandle.product();
0101
0102 siStripClusterInfo_.initEvent(iSetup);
0103
0104 auto number = std::make_unique<std::vector<unsigned>>(7, 0);
0105 auto width = std::make_unique<std::vector<unsigned>>();
0106 auto variance = std::make_unique<std::vector<float>>();
0107 auto barystrip = std::make_unique<std::vector<float>>();
0108 auto middlestrip = std::make_unique<std::vector<float>>();
0109 auto charge = std::make_unique<std::vector<unsigned>>();
0110 auto noise = std::make_unique<std::vector<float>>();
0111 auto ston = std::make_unique<std::vector<float>>();
0112 auto seedstrip = std::make_unique<std::vector<unsigned>>();
0113 auto seedindex = std::make_unique<std::vector<unsigned>>();
0114 auto seedcharge = std::make_unique<std::vector<unsigned>>();
0115 auto seednoise = std::make_unique<std::vector<float>>();
0116 auto seedgain = std::make_unique<std::vector<float>>();
0117 auto qualityisbad = std::make_unique<std::vector<unsigned>>();
0118
0119 auto rawchargeC = std::make_unique<std::vector<float>>();
0120 auto rawchargeL = std::make_unique<std::vector<float>>();
0121 auto rawchargeR = std::make_unique<std::vector<float>>();
0122 auto rawchargeLL = std::make_unique<std::vector<float>>();
0123 auto rawchargeRR = std::make_unique<std::vector<float>>();
0124 auto etaX = std::make_unique<std::vector<float>>();
0125 auto eta = std::make_unique<std::vector<float>>();
0126 auto foldedeta = std::make_unique<std::vector<float>>();
0127 auto etaasymm = std::make_unique<std::vector<float>>();
0128 auto outsideasymm = std::make_unique<std::vector<float>>();
0129 auto neweta = std::make_unique<std::vector<float>>();
0130 auto newetaerr = std::make_unique<std::vector<float>>();
0131
0132 auto detid = std::make_unique<std::vector<unsigned>>();
0133 auto subdetid = std::make_unique<std::vector<int>>();
0134 auto side = std::make_unique<std::vector<int>>();
0135 auto module = std::make_unique<std::vector<int>>();
0136 auto layerwheel = std::make_unique<std::vector<int>>();
0137 auto stringringrod = std::make_unique<std::vector<int>>();
0138 auto petal = std::make_unique<std::vector<int>>();
0139 auto stereo = std::make_unique<std::vector<int>>();
0140
0141 edm::Handle<edmNew::DetSetVector<SiStripCluster>> clusters;
0142
0143 iEvent.getByToken(theClustersToken_, clusters);
0144
0145 edm::Handle<edm::DetSetVector<SiStripProcessedRawDigi>> rawProcessedDigis;
0146
0147 iEvent.getByToken(theDigisToken_, rawProcessedDigis);
0148
0149 edmNew::DetSetVector<SiStripCluster>::const_iterator itClusters = clusters->begin();
0150 for (; itClusters != clusters->end(); ++itClusters) {
0151 uint32_t id = itClusters->id();
0152 const moduleVars moduleV(id, tTopo);
0153 for (edmNew::DetSet<SiStripCluster>::const_iterator cluster = itClusters->begin(); cluster != itClusters->end();
0154 ++cluster) {
0155 siStripClusterInfo_.setCluster(*cluster, id);
0156 const SiStripClusterInfo& info = siStripClusterInfo_;
0157 const NearDigis digis = rawProcessedDigis.isValid() ? NearDigis(info, *rawProcessedDigis) : NearDigis(info);
0158
0159 (number->at(0))++;
0160 (number->at(moduleV.subdetid))++;
0161 width->push_back(cluster->amplitudes().size());
0162 barystrip->push_back(cluster->barycenter());
0163 variance->push_back(info.variance());
0164 middlestrip->push_back(info.firstStrip() + info.width() / 2.0);
0165 charge->push_back(info.charge());
0166 noise->push_back(info.noiseRescaledByGain());
0167 ston->push_back(info.signalOverNoise());
0168 seedstrip->push_back(info.maxStrip());
0169 seedindex->push_back(info.maxIndex());
0170 seedcharge->push_back(info.maxCharge());
0171 seednoise->push_back(info.stripNoisesRescaledByGain().at(info.maxIndex()));
0172 seedgain->push_back(info.stripGains().at(info.maxIndex()));
0173 qualityisbad->push_back(info.IsAnythingBad());
0174
0175 rawchargeC->push_back(digis.max);
0176 rawchargeL->push_back(digis.left);
0177 rawchargeR->push_back(digis.right);
0178 rawchargeLL->push_back(digis.Lleft);
0179 rawchargeRR->push_back(digis.Rright);
0180 etaX->push_back(digis.etaX());
0181 eta->push_back(digis.eta());
0182 etaasymm->push_back(digis.etaasymm());
0183 outsideasymm->push_back(digis.outsideasymm());
0184 neweta->push_back((digis.last - digis.first) / info.charge());
0185 newetaerr->push_back((sqrt(digis.last + digis.first)) / pow(info.charge(), 1.5));
0186
0187 detid->push_back(id);
0188 subdetid->push_back(moduleV.subdetid);
0189 side->push_back(moduleV.side);
0190 module->push_back(moduleV.module);
0191 layerwheel->push_back(moduleV.layerwheel);
0192 stringringrod->push_back(moduleV.stringringrod);
0193 petal->push_back(moduleV.petal);
0194 stereo->push_back(moduleV.stereo);
0195 }
0196 }
0197
0198 iEvent.put(std::move(number), Prefix + "number");
0199 iEvent.put(std::move(width), Prefix + "width");
0200 iEvent.put(std::move(variance), Prefix + "variance");
0201 iEvent.put(std::move(barystrip), Prefix + "barystrip");
0202 iEvent.put(std::move(middlestrip), Prefix + "middlestrip");
0203 iEvent.put(std::move(charge), Prefix + "charge");
0204 iEvent.put(std::move(noise), Prefix + "noise");
0205 iEvent.put(std::move(ston), Prefix + "ston");
0206 iEvent.put(std::move(seedstrip), Prefix + "seedstrip");
0207 iEvent.put(std::move(seedindex), Prefix + "seedindex");
0208 iEvent.put(std::move(seedcharge), Prefix + "seedcharge");
0209 iEvent.put(std::move(seednoise), Prefix + "seednoise");
0210 iEvent.put(std::move(seedgain), Prefix + "seedgain");
0211 iEvent.put(std::move(qualityisbad), Prefix + "qualityisbad");
0212
0213 iEvent.put(std::move(rawchargeC), Prefix + "rawchargeC");
0214 iEvent.put(std::move(rawchargeL), Prefix + "rawchargeL");
0215 iEvent.put(std::move(rawchargeR), Prefix + "rawchargeR");
0216 iEvent.put(std::move(rawchargeLL), Prefix + "rawchargeLL");
0217 iEvent.put(std::move(rawchargeRR), Prefix + "rawchargeRR");
0218 iEvent.put(std::move(etaX), Prefix + "etaX");
0219 iEvent.put(std::move(eta), Prefix + "eta");
0220 iEvent.put(std::move(foldedeta), Prefix + "foldedeta");
0221 iEvent.put(std::move(etaasymm), Prefix + "etaasymm");
0222 iEvent.put(std::move(outsideasymm), Prefix + "outsideasymm");
0223 iEvent.put(std::move(neweta), Prefix + "neweta");
0224 iEvent.put(std::move(newetaerr), Prefix + "newetaerr");
0225
0226 iEvent.put(std::move(detid), Prefix + "detid");
0227 iEvent.put(std::move(subdetid), Prefix + "subdetid");
0228 iEvent.put(std::move(module), Prefix + "module");
0229 iEvent.put(std::move(side), Prefix + "side");
0230 iEvent.put(std::move(layerwheel), Prefix + "layerwheel");
0231 iEvent.put(std::move(stringringrod), Prefix + "stringringrod");
0232 iEvent.put(std::move(petal), Prefix + "petal");
0233 iEvent.put(std::move(stereo), Prefix + "stereo");
0234 }
0235
0236 ShallowClustersProducer::NearDigis::NearDigis(const SiStripClusterInfo& info) {
0237 max = info.maxCharge();
0238 left = info.maxIndex() > uint16_t(0) ? info.stripCharges()[info.maxIndex() - 1] : 0;
0239 Lleft = info.maxIndex() > uint16_t(1) ? info.stripCharges()[info.maxIndex() - 2] : 0;
0240 right = unsigned(info.maxIndex() + 1) < info.stripCharges().size() ? info.stripCharges()[info.maxIndex() + 1] : 0;
0241 Rright = unsigned(info.maxIndex() + 2) < info.stripCharges().size() ? info.stripCharges()[info.maxIndex() + 2] : 0;
0242 first = info.stripCharges()[0];
0243 last = info.stripCharges()[info.width() - 1];
0244 }
0245
0246 ShallowClustersProducer::NearDigis::NearDigis(const SiStripClusterInfo& info,
0247 const edm::DetSetVector<SiStripProcessedRawDigi>& rawProcessedDigis) {
0248 edm::DetSetVector<SiStripProcessedRawDigi>::const_iterator digiframe = rawProcessedDigis.find(info.detId());
0249 if (digiframe != rawProcessedDigis.end()) {
0250 max = digiframe->data.at(info.maxStrip()).adc();
0251 left = info.maxStrip() > uint16_t(0) ? digiframe->data.at(info.maxStrip() - 1).adc() : 0;
0252 Lleft = info.maxStrip() > uint16_t(1) ? digiframe->data.at(info.maxStrip() - 2).adc() : 0;
0253 right = unsigned(info.maxStrip() + 1) < digiframe->data.size() ? digiframe->data.at(info.maxStrip() + 1).adc() : 0;
0254 Rright = unsigned(info.maxStrip() + 2) < digiframe->data.size() ? digiframe->data.at(info.maxStrip() + 2).adc() : 0;
0255 first = digiframe->data.at(info.firstStrip()).adc();
0256 last = digiframe->data.at(info.firstStrip() + info.width() - 1).adc();
0257 } else {
0258 *this = NearDigis(info);
0259 }
0260 }
0261
0262 ShallowClustersProducer::moduleVars::moduleVars(uint32_t detid, const TrackerTopology* tTopo) {
0263 SiStripDetId subdet(detid);
0264 subdetid = subdet.subDetector();
0265 if (SiStripDetId::TIB == subdetid) {
0266 module = tTopo->tibModule(detid);
0267 side = tTopo->tibIsZMinusSide(detid) ? -1 : 1;
0268 layerwheel = tTopo->tibLayer(detid);
0269 stringringrod = tTopo->tibString(detid);
0270 stereo = tTopo->tibIsStereo(detid) ? 1 : 0;
0271 } else if (SiStripDetId::TID == subdetid) {
0272 module = tTopo->tidModule(detid);
0273 side = tTopo->tidIsZMinusSide(detid) ? -1 : 1;
0274 layerwheel = tTopo->tidWheel(detid);
0275 stringringrod = tTopo->tidRing(detid);
0276 stereo = tTopo->tidIsStereo(detid) ? 1 : 0;
0277 } else if (SiStripDetId::TOB == subdetid) {
0278 module = tTopo->tobModule(detid);
0279 side = tTopo->tobIsZMinusSide(detid) ? -1 : 1;
0280 layerwheel = tTopo->tobLayer(detid);
0281 stringringrod = tTopo->tobRod(detid);
0282 stereo = tTopo->tobIsStereo(detid) ? 1 : 0;
0283 } else if (SiStripDetId::TEC == subdetid) {
0284 module = tTopo->tecModule(detid);
0285 side = tTopo->tecIsZMinusSide(detid) ? -1 : 1;
0286 layerwheel = tTopo->tecWheel(detid);
0287 stringringrod = tTopo->tecRing(detid);
0288 petal = tTopo->tecPetalNumber(detid);
0289 stereo = tTopo->tecIsStereo(detid) ? 1 : 0;
0290 } else {
0291 module = 0;
0292 side = 0;
0293 layerwheel = -1;
0294 stringringrod = -1;
0295 petal = -1;
0296 }
0297 }
