MeasurementDet/plugins/MeasurementTrackerEventProducer.cc

0001 #include "MeasurementTrackerEventProducer.h"
0002
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0006
0007 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0008 #include "CondFormats/SiPixelObjects/interface/CablingPathToDetUnit.h"
0009 #include "CondFormats/SiPixelObjects/interface/PixelROC.h"
0010 #include "CondFormats/SiPixelObjects/interface/LocalPixel.h"
0011
0012 #include <algorithm>
0013
0014 MeasurementTrackerEventProducer::MeasurementTrackerEventProducer(const edm::ParameterSet& iConfig)
0015     : measurementTrackerToken_(
0016           esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("measurementTracker")))),
0017       switchOffPixelsIfEmpty_(iConfig.getParameter<bool>("switchOffPixelsIfEmpty")) {
0018   std::vector<edm::InputTag> inactivePixelDetectorTags(
0019       iConfig.getParameter<std::vector<edm::InputTag>>("inactivePixelDetectorLabels"));
0020   for (auto& t : inactivePixelDetectorTags)
0021     theInactivePixelDetectorLabels.push_back(consumes<DetIdCollection>(t));
0022
0023   std::vector<edm::InputTag> badPixelFEDChannelCollectionTags =
0024       iConfig.getParameter<std::vector<edm::InputTag>>("badPixelFEDChannelCollectionLabels");
0025   if (!badPixelFEDChannelCollectionTags.empty()) {
0026     for (auto& t : badPixelFEDChannelCollectionTags)
0027       theBadPixelFEDChannelsLabels.push_back(consumes<PixelFEDChannelCollection>(t));
0028     pixelCablingMapToken_ = esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("pixelCablingMapLabel")));
0029   }
0030
0031   std::vector<edm::InputTag> inactiveStripDetectorTags(
0032       iConfig.getParameter<std::vector<edm::InputTag>>("inactiveStripDetectorLabels"));
0033   for (auto& t : inactiveStripDetectorTags)
0034     theInactiveStripDetectorLabels.push_back(consumes<DetIdCollection>(t));
0035
0036   //the measurement tracking is set to skip clusters, the other option is set from outside
0037   edm::InputTag skip = iConfig.getParameter<edm::InputTag>("skipClusters");
0038   selfUpdateSkipClusters_ = !(skip == edm::InputTag(""));
0039   LogDebug("MeasurementTracker") << "skipping clusters: " << selfUpdateSkipClusters_;
0040   isPhase2_ = false;
0041   useVectorHits_ = false;
0042
0043   if (!iConfig.getParameter<std::string>("stripClusterProducer").empty()) {
0044     theStripClusterLabel = consumes<edmNew::DetSetVector<SiStripCluster>>(
0045         edm::InputTag(iConfig.getParameter<std::string>("stripClusterProducer")));
0046     if (selfUpdateSkipClusters_)
0047       theStripClusterMask = consumes<edm::ContainerMask<edmNew::DetSetVector<SiStripCluster>>>(
0048           iConfig.getParameter<edm::InputTag>("skipClusters"));
0049   }
0050   if (!iConfig.getParameter<std::string>("pixelClusterProducer").empty()) {
0051     thePixelClusterLabel = consumes<edmNew::DetSetVector<SiPixelCluster>>(
0052         edm::InputTag(iConfig.getParameter<std::string>("pixelClusterProducer")));
0053     if (selfUpdateSkipClusters_)
0054       thePixelClusterMask = consumes<edm::ContainerMask<edmNew::DetSetVector<SiPixelCluster>>>(
0055           iConfig.getParameter<edm::InputTag>("skipClusters"));
0056   }
0057   if (!iConfig.getParameter<std::string>("Phase2TrackerCluster1DProducer").empty()) {
0058     thePh2OTClusterLabel = consumes<edmNew::DetSetVector<Phase2TrackerCluster1D>>(
0059         edm::InputTag(iConfig.getParameter<std::string>("Phase2TrackerCluster1DProducer")));
0060     isPhase2_ = true;
0061   }
0062   if (!(iConfig.getParameter<edm::InputTag>("vectorHits") == edm::InputTag("") ||
0063         iConfig.getParameter<edm::InputTag>("vectorHitsRej") == edm::InputTag(""))) {
0064     thePh2OTVectorHitsLabel = consumes<VectorHitCollection>(iConfig.getParameter<edm::InputTag>("vectorHits"));
0065     thePh2OTVectorHitsRejLabel = consumes<VectorHitCollection>(iConfig.getParameter<edm::InputTag>("vectorHitsRej"));
0066     isPhase2_ = true;
0067     useVectorHits_ = true;
0068   }
0069
0070   produces<MeasurementTrackerEvent>();
0071 }
0072
0073 void MeasurementTrackerEventProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0074   edm::ParameterSetDescription desc;
0075
0076   desc.add<std::string>("measurementTracker", "");
0077   desc.add<edm::InputTag>("skipClusters", edm::InputTag());
0078   desc.add<std::string>("pixelClusterProducer", "siPixelClusters");
0079   desc.add<std::string>("stripClusterProducer", "siStripClusters");
0080   desc.add<std::string>("Phase2TrackerCluster1DProducer", "");
0081   desc.add<edm::InputTag>("vectorHits", edm::InputTag(""));
0082   desc.add<edm::InputTag>("vectorHitsRej", edm::InputTag(""));
0083
0084   desc.add<std::vector<edm::InputTag>>("inactivePixelDetectorLabels",
0085                                        std::vector<edm::InputTag>{{edm::InputTag("siPixelDigis")}})
0086       ->setComment("One or more DetIdCollections of modules to mask on the fly for a given event");
0087   desc.add<std::vector<edm::InputTag>>("badPixelFEDChannelCollectionLabels", std::vector<edm::InputTag>())
0088       ->setComment("One or more PixelFEDChannelCollections of modules+ROCs to mask on the fly for a given event");
0089   desc.add<std::string>("pixelCablingMapLabel", "");
0090
0091   desc.add<std::vector<edm::InputTag>>("inactiveStripDetectorLabels",
0092                                        std::vector<edm::InputTag>{{edm::InputTag("siStripDigis")}})
0093       ->setComment("One or more DetIdCollections of modules to mask on the fly for a given event");
0094
0095   desc.add<bool>("switchOffPixelsIfEmpty", true)->setComment("let's keep it like this, for cosmics");
0096
0097   descriptions.add("measurementTrackerEventDefault", desc);
0098 }
0099
0100 void MeasurementTrackerEventProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0101   auto const& measurementTracker = iSetup.getData(measurementTrackerToken_);
0102
0103   // create new data structures from templates
0104   auto stripData = std::make_unique<StMeasurementDetSet>(measurementTracker.stripDetConditions());
0105   auto pixelData = std::make_unique<PxMeasurementDetSet>(measurementTracker.pixelDetConditions());
0106   auto phase2OTData = std::make_unique<Phase2OTMeasurementDetSet>(measurementTracker.phase2DetConditions());
0107
0108   std::vector<bool> stripClustersToSkip;
0109   std::vector<bool> pixelClustersToSkip;
0110   std::vector<bool> phase2ClustersToSkip;
0111   // fill them
0112   updateStrips(iEvent, *stripData, stripClustersToSkip);
0113   updatePixels(iEvent,
0114                *pixelData,
0115                pixelClustersToSkip,
0116                dynamic_cast<const TrackerGeometry&>(*(measurementTracker.geomTracker())),
0117                iSetup);
0118   updatePhase2OT(iEvent, *phase2OTData);
0119   updateStacks(iEvent, *phase2OTData);
0120
0121   // put into MTE
0122   // put into event
0123   //
0124
0125   const VectorHitCollection* phase2OTVectorHits = useVectorHits_ ? &iEvent.get(thePh2OTVectorHitsLabel) : nullptr;
0126   const VectorHitCollection* phase2OTVectorHitsRej = useVectorHits_ ? &iEvent.get(thePh2OTVectorHitsRejLabel) : nullptr;
0127   iEvent.put(std::make_unique<MeasurementTrackerEvent>(measurementTracker,
0128                                                        stripData.release(),
0129                                                        pixelData.release(),
0130                                                        phase2OTData.release(),
0131                                                        phase2OTVectorHits,
0132                                                        phase2OTVectorHitsRej,
0133                                                        stripClustersToSkip,
0134                                                        pixelClustersToSkip,
0135                                                        phase2ClustersToSkip));
0136 }
0137
0138 void MeasurementTrackerEventProducer::updatePixels(const edm::Event& event,
0139                                                    PxMeasurementDetSet& thePxDets,
0140                                                    std::vector<bool>& pixelClustersToSkip,
0141                                                    const TrackerGeometry& trackerGeom,
0142                                                    const edm::EventSetup& iSetup) const {
0143   // start by clearinng everything
0144   thePxDets.setEmpty();
0145
0146   std::vector<uint32_t> rawInactiveDetIds;
0147   if (!theInactivePixelDetectorLabels.empty()) {
0148     edm::Handle<DetIdCollection> detIds;
0149     for (const edm::EDGetTokenT<DetIdCollection>& tk : theInactivePixelDetectorLabels) {
0150       if (event.getByToken(tk, detIds)) {
0151         rawInactiveDetIds.insert(rawInactiveDetIds.end(), detIds->begin(), detIds->end());
0152       } else {
0153         static std::atomic<bool> iFailedAlready{false};
0154         bool expected = false;
0155         if (iFailedAlready.compare_exchange_strong(expected, true, std::memory_order_acq_rel)) {
0156           edm::LogError("MissingProduct")
0157               << "I fail to get the list of inactive pixel modules, because of 4.2/4.4 event content change.";
0158         }
0159       }
0160     }
0161     if (!rawInactiveDetIds.empty())
0162       std::sort(rawInactiveDetIds.begin(), rawInactiveDetIds.end());
0163     // mark as inactive if in rawInactiveDetIds
0164     int i = 0, endDet = thePxDets.size();
0165     unsigned int idp = 0;
0166     for (auto id : rawInactiveDetIds) {
0167       if (id == idp)
0168         continue;  // skip multiple id
0169       idp = id;
0170       i = thePxDets.find(id, i);
0171       assert(i != endDet && id == thePxDets.id(i));
0172       thePxDets.setActiveThisEvent(i, false);
0173     }
0174   }
0175
0176   if (!theBadPixelFEDChannelsLabels.empty()) {
0177     auto const& cablingMap = iSetup.getData(pixelCablingMapToken_);
0178
0179     edm::Handle<PixelFEDChannelCollection> pixelFEDChannelCollectionHandle;
0180     for (const edm::EDGetTokenT<PixelFEDChannelCollection>& tk : theBadPixelFEDChannelsLabels) {
0181       if (!event.getByToken(tk, pixelFEDChannelCollectionHandle))
0182         continue;
0183       int i = 0;
0184       for (const auto& disabledChannels : *pixelFEDChannelCollectionHandle) {
0185         PxMeasurementDetSet::BadFEDChannelPositions positions;
0186         for (const auto& ch : disabledChannels) {
0187           const sipixelobjects::PixelROC *roc_first = nullptr, *roc_last = nullptr;
0188           sipixelobjects::CablingPathToDetUnit path = {ch.fed, ch.link, 0};
0189           // PixelFEDChannelCollection addresses the ROCs by their 'idInDetUnit' (from 0 to 15), ROCs also know their on 'idInDetUnit',
0190           // however the cabling map uses a numbering [1,numberOfROCs], see sipixelobjects::PixelFEDLink::roc(unsigned int id), not necessarily sorted in the same direction.
0191           // PixelFEDChannelCollection MUST be filled such that ch.roc_first (ch.roc_last) correspond to the lowest (highest) 'idInDetUnit' in the channel
0192           for (path.roc = 1; path.roc <= (ch.roc_last - ch.roc_first) + 1; path.roc++) {
0193             const sipixelobjects::PixelROC* roc = cablingMap.findItem(path);
0194             if (roc == nullptr)
0195               continue;
0196             assert(roc->rawId() == disabledChannels.detId());
0197             if (roc->idInDetUnit() == ch.roc_first)
0198               roc_first = roc;
0199             if (roc->idInDetUnit() == ch.roc_last)
0200               roc_last = roc;
0201           }
0202           if (roc_first == nullptr || roc_last == nullptr) {
0203             edm::LogError("PixelFEDChannelCollection")
0204                 << "Do not find either roc_first or roc_last in the cabling map.";
0205             continue;
0206           }
0207           const PixelGeomDetUnit* theGeomDet =
0208               dynamic_cast<const PixelGeomDetUnit*>(trackerGeom.idToDet(roc_first->rawId()));
0209           PixelTopology const* topology = &(theGeomDet->specificTopology());
0210           sipixelobjects::LocalPixel::RocRowCol local = {
0211               topology->rowsperroc() / 2, topology->colsperroc() / 2};  //corresponding to center of ROC row, col
0212           sipixelobjects::GlobalPixel global = roc_first->toGlobal(sipixelobjects::LocalPixel(local));
0213           LocalPoint lp1 = topology->localPosition(MeasurementPoint(global.row, global.col));
0214           global = roc_last->toGlobal(sipixelobjects::LocalPixel(local));
0215           LocalPoint lp2 = topology->localPosition(MeasurementPoint(global.row, global.col));
0216           LocalPoint ll(std::min(lp1.x(), lp2.x()), std::min(lp1.y(), lp2.y()), std::min(lp1.z(), lp2.z()));
0217           LocalPoint ur(std::max(lp1.x(), lp2.x()), std::max(lp1.y(), lp2.y()), std::max(lp1.z(), lp2.z()));
0218           positions.push_back(std::make_pair(ll, ur));
0219         }  // loop on channels
0220         if (!positions.empty()) {
0221           i = thePxDets.find(disabledChannels.detId(), i);
0222           assert(i != thePxDets.size() && thePxDets.id(i) == disabledChannels.detId());
0223           thePxDets.addBadFEDChannelPositions(i, positions);
0224         }
0225       }  // loop on DetId-s
0226     }    // loop on labels
0227   }      // if collection labels are populated
0228
0229   // Pixel Clusters
0230   if (thePixelClusterLabel.isUninitialized()) {  //clusters have not been produced
0231     if (switchOffPixelsIfEmpty_) {
0232       thePxDets.setActiveThisEvent(false);
0233     }
0234   } else {
0235     edm::Handle<edmNew::DetSetVector<SiPixelCluster>>& pixelClusters = thePxDets.handle();
0236     if (event.getByToken(thePixelClusterLabel, pixelClusters)) {
0237       const edmNew::DetSetVector<SiPixelCluster>* pixelCollection = pixelClusters.product();
0238
0239       if (switchOffPixelsIfEmpty_ && pixelCollection->empty()) {
0240         thePxDets.setActiveThisEvent(false);
0241       } else {
0242         //std::cout <<"updatePixels "<<pixelCollection->dataSize()<<std::endl;
0243         pixelClustersToSkip.resize(pixelCollection->dataSize());
0244         std::fill(pixelClustersToSkip.begin(), pixelClustersToSkip.end(), false);
0245
0246         if (selfUpdateSkipClusters_) {
0247           edm::Handle<edm::ContainerMask<edmNew::DetSetVector<SiPixelCluster>>> pixelClusterMask;
0248           //and get the collection of pixel ref to skip
0249           event.getByToken(thePixelClusterMask, pixelClusterMask);
0250           LogDebug("MeasurementTracker") << "getting pxl refs to skip";
0251           if (pixelClusterMask.failedToGet())
0252             edm::LogError("MeasurementTracker") << "not getting the pixel clusters to skip";
0253           if (pixelClusterMask->refProd().id() != pixelClusters.id()) {
0254             edm::LogError("ProductIdMismatch")
0255                 << "The pixel masking does not point to the proper collection of clusters: "
0256                 << pixelClusterMask->refProd().id() << "!=" << pixelClusters.id();
0257           }
0258           pixelClusterMask->copyMaskTo(pixelClustersToSkip);
0259         }
0260
0261         // FIXME: should check if lower_bound is better
0262         int i = 0, endDet = thePxDets.size();
0263         for (edmNew::DetSetVector<SiPixelCluster>::const_iterator it = pixelCollection->begin(),
0264                                                                   ed = pixelCollection->end();
0265              it != ed;
0266              ++it) {
0267           edmNew::DetSet<SiPixelCluster> set(*it);
0268           unsigned int id = set.id();
0269           while (id != thePxDets.id(i)) {
0270             ++i;
0271             if (endDet == i)
0272               throw "we have a problem!!!!";
0273           }
0274           // push cluster range in det
0275           if (thePxDets.isActive(i)) {
0276             thePxDets.update(i, set);
0277           }
0278         }
0279       }
0280     } else {
0281       edm::EDConsumerBase::Labels labels;
0282       labelsForToken(thePixelClusterLabel, labels);
0283       edm::LogWarning("MeasurementTrackerEventProducer")
0284           << "input pixel clusters collection " << labels.module << " is not valid";
0285     }
0286   }
0287 }
0288
0289 void MeasurementTrackerEventProducer::updateStrips(const edm::Event& event,
0290                                                    StMeasurementDetSet& theStDets,
0291                                                    std::vector<bool>& stripClustersToSkip) const {
0292   typedef edmNew::DetSet<SiStripCluster> StripDetSet;
0293
0294   std::vector<uint32_t> rawInactiveDetIds;
0295   getInactiveStrips(event, rawInactiveDetIds);
0296
0297   // Strip Clusters
0298   //first clear all of them
0299   theStDets.setEmpty();
0300
0301   if (theStripClusterLabel.isUninitialized())
0302     return;  //clusters have not been produced
0303
0304   const int endDet = theStDets.size();
0305
0306   // mark as inactive if in rawInactiveDetIds
0307   int i = 0;
0308   unsigned int idp = 0;
0309   for (auto id : rawInactiveDetIds) {
0310     if (id == idp)
0311       continue;  // skip multiple id
0312     idp = id;
0313     i = theStDets.find(id, i);
0314     assert(i != endDet && id == theStDets.id(i));
0315     theStDets.setActiveThisEvent(i, false);
0316   }
0317
0318   //=========  actually load cluster =============
0319   {
0320     edm::Handle<edmNew::DetSetVector<SiStripCluster>> clusterHandle;
0321     if (event.getByToken(theStripClusterLabel, clusterHandle)) {
0322       const edmNew::DetSetVector<SiStripCluster>* clusterCollection = clusterHandle.product();
0323
0324       if (selfUpdateSkipClusters_) {
0325         edm::Handle<edm::ContainerMask<edmNew::DetSetVector<SiStripCluster>>> stripClusterMask;
0326         //and get the collection of pixel ref to skip
0327         LogDebug("MeasurementTracker") << "getting strp refs to skip";
0328         event.getByToken(theStripClusterMask, stripClusterMask);
0329         if (stripClusterMask.failedToGet())
0330           edm::LogError("MeasurementTracker") << "not getting the strip clusters to skip";
0331         if (stripClusterMask->refProd().id() != clusterHandle.id()) {
0332           edm::LogError("ProductIdMismatch")
0333               << "The strip masking does not point to the proper collection of clusters: "
0334               << stripClusterMask->refProd().id() << "!=" << clusterHandle.id();
0335         }
0336         stripClusterMask->copyMaskTo(stripClustersToSkip);
0337       }
0338
0339       theStDets.handle() = clusterHandle;
0340       int i = 0;
0341       // cluster and det and in order (both) and unique so let's use set intersection
0342       for (auto j = 0U; j < (*clusterCollection).size(); ++j) {
0343         unsigned int id = (*clusterCollection).id(j);
0344         while (id != theStDets.id(i)) {  // eventually change to lower_bound
0345           ++i;
0346           if (endDet == i)
0347             throw "we have a problem in strips!!!!";
0348         }
0349
0350         // push cluster range in det
0351         if (theStDets.isActive(i))
0352           theStDets.update(i, j);
0353       }
0354     } else {
0355       edm::EDConsumerBase::Labels labels;
0356       labelsForToken(theStripClusterLabel, labels);
0357       edm::LogWarning("MeasurementTrackerEventProducer")
0358           << "input strip cluster collection " << labels.module << " is not valid";
0359     }
0360   }
0361 }
0362
0363 //FIXME: just a temporary solution for phase2!
0364 void MeasurementTrackerEventProducer::updatePhase2OT(const edm::Event& event,
0365                                                      Phase2OTMeasurementDetSet& thePh2OTDets) const {
0366   thePh2OTDets.setEmpty();
0367
0368   // Phase2OT Clusters
0369   if (isPhase2_) {
0370     if (thePh2OTClusterLabel.isUninitialized()) {  //clusters have not been produced
0371       thePh2OTDets.setActiveThisEvent(false);
0372     } else {
0373       edm::Handle<edmNew::DetSetVector<Phase2TrackerCluster1D>>& phase2OTClusters = thePh2OTDets.handle();
0374       if (event.getByToken(thePh2OTClusterLabel, phase2OTClusters)) {
0375         const edmNew::DetSetVector<Phase2TrackerCluster1D>* phase2OTCollection = phase2OTClusters.product();
0376
0377         int i = 0, endDet = thePh2OTDets.size();
0378         for (edmNew::DetSetVector<Phase2TrackerCluster1D>::const_iterator it = phase2OTCollection->begin(),
0379                                                                           ed = phase2OTCollection->end();
0380              it != ed;
0381              ++it) {
0382           edmNew::DetSet<Phase2TrackerCluster1D> set(*it);
0383           unsigned int id = set.id();
0384           while (id != thePh2OTDets.id(i)) {
0385             ++i;
0386             if (endDet == i)
0387               throw "we have a problem!!!!";
0388           }
0389           // push cluster range in det
0390           if (thePh2OTDets.isActive(i)) {
0391             thePh2OTDets.update(i, set);
0392           }
0393         }
0394       } else {
0395         edm::EDConsumerBase::Labels labels;
0396         labelsForToken(thePh2OTClusterLabel, labels);
0397         edm::LogWarning("MeasurementTrackerEventProducer")
0398             << "input Phase2TrackerCluster1D collection " << labels.module << " is not valid";
0399       }
0400     }
0401   }
0402   return;
0403 }
0404
0405 void MeasurementTrackerEventProducer::getInactiveStrips(const edm::Event& event,
0406                                                         std::vector<uint32_t>& rawInactiveDetIds) const {
0407   if (!theInactiveStripDetectorLabels.empty()) {
0408     edm::Handle<DetIdCollection> detIds;
0409     for (const edm::EDGetTokenT<DetIdCollection>& tk : theInactiveStripDetectorLabels) {
0410       if (event.getByToken(tk, detIds)) {
0411         rawInactiveDetIds.insert(rawInactiveDetIds.end(), detIds->begin(), detIds->end());
0412       }
0413     }
0414     if (!rawInactiveDetIds.empty())
0415       std::sort(rawInactiveDetIds.begin(), rawInactiveDetIds.end());
0416   }
0417 }
0418
0419 DEFINE_FWK_MODULE(MeasurementTrackerEventProducer);