Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:14:45

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   //Retrieve tracker topology from geometry
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   //  iEvent.getByLabel(theClustersLabel, clusters);
0143   iEvent.getByToken(theClustersToken_, clusters);
0144 
0145   edm::Handle<edm::DetSetVector<SiStripProcessedRawDigi>> rawProcessedDigis;
0146   //  iEvent.getByLabel("siStripProcessedRawDigis", "", rawProcessedDigis);
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 }