MeasurementDet/plugins/TkStackMeasurementDet.cc

0001 #include "TkStackMeasurementDet.h"
0002
0003 #include "TrackingTools/MeasurementDet/interface/MeasurementDetException.h"
0004
0005 using namespace std;
0006
0007 TkStackMeasurementDet::TkStackMeasurementDet(const StackGeomDet* gdet, const PixelClusterParameterEstimator* cpe)
0008     : MeasurementDet(gdet), thePixelCPE(cpe), theLowerDet(nullptr), theUpperDet(nullptr) {}
0009
0010 void TkStackMeasurementDet::init(const MeasurementDet* lowerDet, const MeasurementDet* upperDet) {
0011   theLowerDet = dynamic_cast<const TkPhase2OTMeasurementDet*>(lowerDet);
0012   theUpperDet = dynamic_cast<const TkPhase2OTMeasurementDet*>(upperDet);
0013
0014   if ((theLowerDet == nullptr) || (theUpperDet == nullptr)) {
0015     throw MeasurementDetException(
0016         "TkStackMeasurementDet ERROR: Trying to glue a det which is not a TkPhase2OTMeasurementDet");
0017   }
0018 }
0019
0020 TkStackMeasurementDet::RecHitContainer TkStackMeasurementDet::recHits(const TrajectoryStateOnSurface& ts,
0021                                                                       const MeasurementTrackerEvent& data) const {
0022   RecHitContainer result;
0023
0024   if (data.phase2OTVectorHits().empty())
0025     return result;
0026   LogTrace("MeasurementTracker") << " is not empty";
0027   if (!isActive(data))
0028     return result;
0029   LogTrace("MeasurementTracker") << " and is active";
0030
0031   //find clusters to skip
0032   const detset& lowerDetSet = data.phase2OTData().detSet(lowerDet()->index());
0033   const detset& upperDetSet = data.phase2OTData().detSet(upperDet()->index());
0034   std::vector<bool> skipClustersUpper(data.phase2OTClustersToSkip().empty() ? 0 : upperDetSet.size(), false);
0035   std::vector<bool> skipClustersLower(data.phase2OTClustersToSkip().empty() ? 0 : lowerDetSet.size(), false);
0036
0037   const Phase2TrackerCluster1D* begin = nullptr;
0038   if (!data.phase2OTData().handle()->data().empty()) {
0039     begin = &(data.phase2OTData().handle()->data().front());
0040   }
0041   if (!data.phase2OTClustersToSkip().empty()) {
0042     if (!lowerDetSet.empty()) {
0043       for (const_iterator cil = lowerDetSet.begin(); cil != lowerDetSet.end(); ++cil) {
0044         if (cil < begin) {
0045           edm::LogError("IndexMisMatch") << "TkStackMeasurementDet cannot create hit because of index mismatch.";
0046           return result;
0047         }
0048         unsigned int indexl = cil - begin;
0049         if (data.phase2OTClustersToSkip()[indexl]) {
0050           int iLocalL = std::distance(lowerDetSet.begin(), cil);
0051           skipClustersLower[iLocalL] = true;
0052         }
0053       }
0054     }
0055     if (!upperDetSet.empty()) {
0056       for (const_iterator ciu = upperDetSet.begin(); ciu != upperDetSet.end(); ++ciu) {
0057         if (ciu < begin) {
0058           edm::LogError("IndexMisMatch") << "TkStackMeasurementDet cannot create hit because of index mismatch.";
0059           return result;
0060         }
0061         unsigned int indexu = ciu - begin;
0062         if (data.phase2OTClustersToSkip()[indexu]) {
0063           int iLocalU = std::distance(upperDetSet.begin(), ciu);
0064           skipClustersUpper[iLocalU] = true;
0065         }
0066       }
0067     }
0068   }
0069   DetId detIdStack = specificGeomDet().geographicalId();
0070
0071   auto iterator = data.phase2OTVectorHits().find(detIdStack);
0072   if (iterator == data.phase2OTVectorHits().end())
0073     return result;
0074   for (const auto& vecHit : data.phase2OTVectorHits()[detIdStack]) {
0075     if (!data.phase2OTClustersToSkip().empty()) {
0076       if (skipClustersLower[vecHit.lowerCluster().key() - lowerDetSet.offset()])
0077         continue;
0078       if (skipClustersUpper[vecHit.upperCluster().key() - upperDetSet.offset()])
0079         continue;
0080     }
0081     result.push_back(std::make_shared<VectorHit>(vecHit));
0082   }
0083
0084   iterator = data.phase2OTVectorHitsRej().find(detIdStack);
0085   if (iterator == data.phase2OTVectorHitsRej().end())
0086     return result;
0087   for (const auto& vecHit : data.phase2OTVectorHitsRej()[detIdStack]) {
0088     if (!data.phase2OTClustersToSkip().empty()) {
0089       if (skipClustersLower[vecHit.lowerCluster().key() - lowerDetSet.offset()])
0090         continue;
0091       if (skipClustersUpper[vecHit.upperCluster().key() - upperDetSet.offset()])
0092         continue;
0093     }
0094     result.push_back(std::make_shared<VectorHit>(vecHit));
0095   }
0096
0097   return result;
0098 }
0099
0100 bool TkStackMeasurementDet::measurements(const TrajectoryStateOnSurface& stateOnThisDet,
0101                                          const MeasurementEstimator& est,
0102                                          const MeasurementTrackerEvent& data,
0103                                          TempMeasurements& result) const {
0104   LogDebug("MeasurementTracker") << "TkStackMeasurementDet::measurements";
0105
0106   if (!isActive(data)) {
0107     result.add(theInactiveHit, 0.F);
0108     return true;
0109   }
0110
0111   LogTrace("MeasurementTracker") << " is active";
0112
0113   auto oldSize = result.size();
0114   MeasurementDet::RecHitContainer&& allHits = recHits(stateOnThisDet, data);
0115
0116   for (auto&& hit : allHits) {
0117     std::pair<bool, double> diffEst = est.estimate(stateOnThisDet, *hit);
0118     if (diffEst.first) {
0119       LogDebug("MeasurementTracker") << "New vh added with chi2: " << diffEst.second;
0120       result.add(std::move(hit), diffEst.second);
0121     }
0122   }
0123
0124   if (result.size() > oldSize)
0125     return true;
0126
0127   // create a TrajectoryMeasurement with an invalid RecHit and zero estimate
0128   result.add(theMissingHit, 0.F);
0129   LogDebug("MeasurementTracker") << "adding missing hit";
0130   return false;
0131 }