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<DetIdVector>(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 DetIdVectors 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 (not 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           assert(ch.roc_last >= ch.roc_first);
0193           for (path.roc = 1; path.roc <= (ch.roc_last - ch.roc_first) + 1; ++path.roc) {
0194             const sipixelobjects::PixelROC* roc = cablingMap.findItem(path);
0195             if (roc == nullptr)
0196               continue;
0197             assert(roc->rawId() == disabledChannels.detId());
0198             if (roc->idInDetUnit() == ch.roc_first)
0199               roc_first = roc;
0200             if (roc->idInDetUnit() == ch.roc_last)
0201               roc_last = roc;
0202           }
0203           if (roc_first == nullptr || roc_last == nullptr) {
0204             edm::LogError("PixelFEDChannelCollection")
0205                 << "Do not find either roc_first or roc_last in the cabling map.";
0206             continue;
0207           }
0208           const PixelGeomDetUnit* theGeomDet =
0209               dynamic_cast<const PixelGeomDetUnit*>(trackerGeom.idToDet(roc_first->rawId()));
0210           PixelTopology const* topology = &(theGeomDet->specificTopology());
0211           sipixelobjects::LocalPixel::RocRowCol local = {
0212               topology->rowsperroc() / 2, topology->colsperroc() / 2};  //corresponding to center of ROC row, col
0213           sipixelobjects::GlobalPixel global = roc_first->toGlobal(sipixelobjects::LocalPixel(local));
0214           LocalPoint lp1 = topology->localPosition(MeasurementPoint(global.row, global.col));
0215           global = roc_last->toGlobal(sipixelobjects::LocalPixel(local));
0216           LocalPoint lp2 = topology->localPosition(MeasurementPoint(global.row, global.col));
0217           LocalPoint ll(std::min(lp1.x(), lp2.x()), std::min(lp1.y(), lp2.y()), std::min(lp1.z(), lp2.z()));
0218           LocalPoint ur(std::max(lp1.x(), lp2.x()), std::max(lp1.y(), lp2.y()), std::max(lp1.z(), lp2.z()));
0219           positions.push_back(std::make_pair(ll, ur));
0220         }  // loop on channels
0221         if (not positions.empty()) {
0222           i = thePxDets.find(disabledChannels.detId(), i);
0223           assert(i != thePxDets.size() && thePxDets.id(i) == disabledChannels.detId());
0224           thePxDets.addBadFEDChannelPositions(i, positions);
0225         }
0226       }  // loop on DetId-s
0227     }  // loop on labels
0228   }  // if collection labels are populated
0229
0230   // Pixel Clusters
0231   if (thePixelClusterLabel.isUninitialized()) {  //clusters have not been produced
0232     if (switchOffPixelsIfEmpty_) {
0233       thePxDets.setActiveThisEvent(false);
0234     }
0235   } else {
0236     edm::Handle<edmNew::DetSetVector<SiPixelCluster>>& pixelClusters = thePxDets.handle();
0237     if (event.getByToken(thePixelClusterLabel, pixelClusters)) {
0238       const edmNew::DetSetVector<SiPixelCluster>* pixelCollection = pixelClusters.product();
0239
0240       if (switchOffPixelsIfEmpty_ && pixelCollection->empty()) {
0241         thePxDets.setActiveThisEvent(false);
0242       } else {
0243         //std::cout <<"updatePixels "<<pixelCollection->dataSize()<<std::endl;
0244         pixelClustersToSkip.resize(pixelCollection->dataSize());
0245         std::fill(pixelClustersToSkip.begin(), pixelClustersToSkip.end(), false);
0246
0247         if (selfUpdateSkipClusters_) {
0248           edm::Handle<edm::ContainerMask<edmNew::DetSetVector<SiPixelCluster>>> pixelClusterMask;
0249           //and get the collection of pixel ref to skip
0250           event.getByToken(thePixelClusterMask, pixelClusterMask);
0251           LogDebug("MeasurementTracker") << "getting pxl refs to skip";
0252           if (pixelClusterMask.failedToGet())
0253             edm::LogError("MeasurementTracker") << "not getting the pixel clusters to skip";
0254           if (pixelClusterMask->refProd().id() != pixelClusters.id()) {
0255             edm::LogError("ProductIdMismatch")
0256                 << "The pixel masking does not point to the proper collection of clusters: "
0257                 << pixelClusterMask->refProd().id() << "!=" << pixelClusters.id();
0258           }
0259           pixelClusterMask->copyMaskTo(pixelClustersToSkip);
0260         }
0261
0262         // FIXME: should check if lower_bound is better
0263         int i = 0, endDet = thePxDets.size();
0264         for (edmNew::DetSetVector<SiPixelCluster>::const_iterator it = pixelCollection->begin(),
0265                                                                   ed = pixelCollection->end();
0266              it != ed;
0267              ++it) {
0268           edmNew::DetSet<SiPixelCluster> set(*it);
0269           unsigned int id = set.id();
0270           while (id != thePxDets.id(i)) {
0271             ++i;
0272             if (endDet == i)
0273               throw "we have a problem!!!!";
0274           }
0275           // push cluster range in det
0276           if (thePxDets.isActive(i)) {
0277             thePxDets.update(i, set);
0278           }
0279         }
0280       }
0281     } else {
0282       edm::EDConsumerBase::Labels labels;
0283       labelsForToken(thePixelClusterLabel, labels);
0284       edm::LogWarning("MeasurementTrackerEventProducer")
0285           << "input pixel clusters collection " << labels.module << " is not valid";
0286     }
0287   }
0288 }
0289
0290 void MeasurementTrackerEventProducer::updateStrips(const edm::Event& event,
0291                                                    StMeasurementDetSet& theStDets,
0292                                                    std::vector<bool>& stripClustersToSkip) const {
0293   typedef edmNew::DetSet<SiStripCluster> StripDetSet;
0294
0295   std::vector<uint32_t> rawInactiveDetIds;
0296   getInactiveStrips(event, rawInactiveDetIds);
0297
0298   // Strip Clusters
0299   //first clear all of them
0300   theStDets.setEmpty();
0301
0302   if (theStripClusterLabel.isUninitialized())
0303     return;  //clusters have not been produced
0304
0305   const int endDet = theStDets.size();
0306
0307   // mark as inactive if in rawInactiveDetIds
0308   int i = 0;
0309   unsigned int idp = 0;
0310   for (auto id : rawInactiveDetIds) {
0311     if (id == idp)
0312       continue;  // skip multiple id
0313     idp = id;
0314     i = theStDets.find(id, i);
0315     assert(i != endDet && id == theStDets.id(i));
0316     theStDets.setActiveThisEvent(i, false);
0317   }
0318
0319   //=========  actually load cluster =============
0320   {
0321     edm::Handle<edmNew::DetSetVector<SiStripCluster>> clusterHandle;
0322     if (event.getByToken(theStripClusterLabel, clusterHandle)) {
0323       const edmNew::DetSetVector<SiStripCluster>* clusterCollection = clusterHandle.product();
0324
0325       if (selfUpdateSkipClusters_) {
0326         edm::Handle<edm::ContainerMask<edmNew::DetSetVector<SiStripCluster>>> stripClusterMask;
0327         //and get the collection of pixel ref to skip
0328         LogDebug("MeasurementTracker") << "getting strp refs to skip";
0329         event.getByToken(theStripClusterMask, stripClusterMask);
0330         if (stripClusterMask.failedToGet())
0331           edm::LogError("MeasurementTracker") << "not getting the strip clusters to skip";
0332         if (stripClusterMask->refProd().id() != clusterHandle.id()) {
0333           edm::LogError("ProductIdMismatch")
0334               << "The strip masking does not point to the proper collection of clusters: "
0335               << stripClusterMask->refProd().id() << "!=" << clusterHandle.id();
0336         }
0337         stripClusterMask->copyMaskTo(stripClustersToSkip);
0338       }
0339
0340       theStDets.handle() = clusterHandle;
0341       int i = 0;
0342       // cluster and det and in order (both) and unique so let's use set intersection
0343       for (auto j = 0U; j < (*clusterCollection).size(); ++j) {
0344         unsigned int id = (*clusterCollection).id(j);
0345         while (id != theStDets.id(i)) {  // eventually change to lower_bound
0346           ++i;
0347           if (endDet == i)
0348             throw "we have a problem in strips!!!!";
0349         }
0350
0351         // push cluster range in det
0352         if (theStDets.isActive(i))
0353           theStDets.update(i, j);
0354       }
0355     } else {
0356       edm::EDConsumerBase::Labels labels;
0357       labelsForToken(theStripClusterLabel, labels);
0358       edm::LogWarning("MeasurementTrackerEventProducer")
0359           << "input strip cluster collection " << labels.module << " is not valid";
0360     }
0361   }
0362 }
0363
0364 //FIXME: just a temporary solution for phase2!
0365 void MeasurementTrackerEventProducer::updatePhase2OT(const edm::Event& event,
0366                                                      Phase2OTMeasurementDetSet& thePh2OTDets) const {
0367   thePh2OTDets.setEmpty();
0368
0369   // Phase2OT Clusters
0370   if (isPhase2_) {
0371     if (thePh2OTClusterLabel.isUninitialized()) {  //clusters have not been produced
0372       thePh2OTDets.setActiveThisEvent(false);
0373     } else {
0374       edm::Handle<edmNew::DetSetVector<Phase2TrackerCluster1D>>& phase2OTClusters = thePh2OTDets.handle();
0375       if (event.getByToken(thePh2OTClusterLabel, phase2OTClusters)) {
0376         const edmNew::DetSetVector<Phase2TrackerCluster1D>* phase2OTCollection = phase2OTClusters.product();
0377
0378         int i = 0, endDet = thePh2OTDets.size();
0379         for (edmNew::DetSetVector<Phase2TrackerCluster1D>::const_iterator it = phase2OTCollection->begin(),
0380                                                                           ed = phase2OTCollection->end();
0381              it != ed;
0382              ++it) {
0383           edmNew::DetSet<Phase2TrackerCluster1D> set(*it);
0384           unsigned int id = set.id();
0385           while (id != thePh2OTDets.id(i)) {
0386             ++i;
0387             if (endDet == i)
0388               throw "we have a problem!!!!";
0389           }
0390           // push cluster range in det
0391           if (thePh2OTDets.isActive(i)) {
0392             thePh2OTDets.update(i, set);
0393           }
0394         }
0395       } else {
0396         edm::EDConsumerBase::Labels labels;
0397         labelsForToken(thePh2OTClusterLabel, labels);
0398         edm::LogWarning("MeasurementTrackerEventProducer")
0399             << "input Phase2TrackerCluster1D collection " << labels.module << " is not valid";
0400       }
0401     }
0402   }
0403   return;
0404 }
0405
0406 void MeasurementTrackerEventProducer::getInactiveStrips(const edm::Event& event,
0407                                                         std::vector<uint32_t>& rawInactiveDetIds) const {
0408   if (!theInactiveStripDetectorLabels.empty()) {
0409     edm::Handle<DetIdVector> detIds;
0410     for (const edm::EDGetTokenT<DetIdVector>& tk : theInactiveStripDetectorLabels) {
0411       if (event.getByToken(tk, detIds)) {
0412         rawInactiveDetIds.insert(rawInactiveDetIds.end(), detIds->begin(), detIds->end());
0413       }
0414     }
0415     if (!rawInactiveDetIds.empty())
0416       std::sort(rawInactiveDetIds.begin(), rawInactiveDetIds.end());
0417   }
0418 }
0419
0420 DEFINE_FWK_MODULE(MeasurementTrackerEventProducer);