diff --git a/Configuration/StandardSequences/python/DigiToRaw_Repack_cff.py b/Configuration/StandardSequences/python/DigiToRaw_Repack_cff.py index 24af3647931b1..56bf9825f1c8f 100644 --- a/Configuration/StandardSequences/python/DigiToRaw_Repack_cff.py +++ b/Configuration/StandardSequences/python/DigiToRaw_Repack_cff.py @@ -72,7 +72,7 @@ cms.InputTag('siStripZeroSuppressionHLT','ScopeMode')), ) -from RecoLocalTracker.SiStripClusterizer.SiStripClusters2ApproxClusters_cff import hltSiStripClusters2ApproxClusters +from RecoLocalTracker.SiStripClusterizer.SiStripClusters2ApproxClusters_cff import * from EventFilter.Utilities.EvFFEDExcluder_cfi import EvFFEDExcluder as _EvFFEDExcluder rawPrimeDataRepacker = _EvFFEDExcluder.clone( @@ -80,5 +80,9 @@ fedsToExclude = [foo for foo in range(50, 490)] ) -DigiToApproxClusterRawTask = cms.Task(siStripDigisHLT,siStripZeroSuppressionHLT,siStripClustersHLT,hltSiStripClusters2ApproxClusters,rawPrimeDataRepacker) +hltScalersRawToDigi = cms.EDProducer( "ScalersRawToDigi", + scalersInputTag = cms.InputTag( "rawDataRepacker" ) +) + +DigiToApproxClusterRawTask = cms.Task(siStripDigisHLT,siStripZeroSuppressionHLT,hltScalersRawToDigi,hltBeamSpotProducer,siStripClustersHLT,hltSiStripClusters2ApproxClusters,rawPrimeDataRepacker) DigiToApproxClusterRaw = cms.Sequence(DigiToApproxClusterRawTask) diff --git a/DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h b/DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h index 8ff3317666376..baf87f791d9f6 100644 --- a/DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h +++ b/DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h @@ -8,27 +8,38 @@ class SiStripApproximateCluster { public: SiStripApproximateCluster() {} - explicit SiStripApproximateCluster(cms_uint16_t barycenter, - cms_uint8_t width, - cms_uint8_t avgCharge, - bool isSaturated) { + explicit SiStripApproximateCluster( + cms_uint16_t barycenter, cms_uint8_t width, cms_uint8_t avgCharge, bool filter, bool isSaturated) { barycenter_ = barycenter; width_ = width; avgCharge_ = avgCharge; + filter_ = filter; isSaturated_ = isSaturated; } - explicit SiStripApproximateCluster(const SiStripCluster& cluster, unsigned int maxNSat); + explicit SiStripApproximateCluster(const SiStripCluster& cluster, + unsigned int maxNSat, + float hitPredPos, + bool peakFilter); cms_uint16_t barycenter() const { return barycenter_; } cms_uint8_t width() const { return width_; } cms_uint8_t avgCharge() const { return avgCharge_; } + bool filter() const { return filter_; } bool isSaturated() const { return isSaturated_; } + bool peakFilter() const { return peakFilter_; } private: cms_uint16_t barycenter_ = 0; cms_uint8_t width_ = 0; cms_uint8_t avgCharge_ = 0; + bool filter_ = false; bool isSaturated_ = false; + bool peakFilter_ = false; + static constexpr double trimMaxADC_ = 30.; + static constexpr double trimMaxFracTotal_ = .15; + static constexpr double trimMaxFracNeigh_ = .25; + static constexpr double maxTrimmedSizeDiffNeg_ = .7; + static constexpr double maxTrimmedSizeDiffPos_ = 1.; }; #endif // DataFormats_SiStripCluster_SiStripApproximateCluster_h diff --git a/DataFormats/SiStripCluster/interface/SiStripCluster.h b/DataFormats/SiStripCluster/interface/SiStripCluster.h index 383e3960f629e..d95f379251036 100644 --- a/DataFormats/SiStripCluster/interface/SiStripCluster.h +++ b/DataFormats/SiStripCluster/interface/SiStripCluster.h @@ -82,6 +82,10 @@ class SiStripCluster { */ int charge() const; + bool filter() const; + + bool isFromApprox() const; + /** Test (set) the merged status of the cluster * */ @@ -99,6 +103,7 @@ class SiStripCluster { //these are used if amplitude information is not available (using approximate cluster constructor) float barycenter_ = 0; int charge_ = 0; + bool filter_ = false; // ggiurgiu@fnal.gov, 01/05/12 // Add cluster errors to be used by rechits from split clusters. diff --git a/DataFormats/SiStripCluster/src/SiStripApproximateCluster.cc b/DataFormats/SiStripCluster/src/SiStripApproximateCluster.cc index 10a7c29c2f60c..01b30963a083a 100644 --- a/DataFormats/SiStripCluster/src/SiStripApproximateCluster.cc +++ b/DataFormats/SiStripCluster/src/SiStripApproximateCluster.cc @@ -3,11 +3,16 @@ #include #include -SiStripApproximateCluster::SiStripApproximateCluster(const SiStripCluster& cluster, unsigned int maxNSat) { +SiStripApproximateCluster::SiStripApproximateCluster(const SiStripCluster& cluster, + unsigned int maxNSat, + float hitPredPos, + bool peakFilter) { barycenter_ = std::round(cluster.barycenter() * 10); width_ = cluster.size(); avgCharge_ = cluster.charge() / cluster.size(); + filter_ = false; isSaturated_ = false; + peakFilter_ = peakFilter; //mimicing the algorithm used in StripSubClusterShapeTrajectoryFilter... //Looks for 3 adjacent saturated strips (ADC>=254) @@ -25,6 +30,28 @@ SiStripApproximateCluster::SiStripApproximateCluster(const SiStripCluster& clust maxSat = std::max(maxSat, thisSat); } if (maxSat >= maxNSat) { + filter_ = true; isSaturated_ = true; } + + unsigned int hitStripsTrim = ampls.size(); + int sum = std::accumulate(ampls.begin(), ampls.end(), 0); + uint8_t trimCut = std::min(trimMaxADC_, std::floor(trimMaxFracTotal_ * sum)); + auto begin = ampls.begin(); + auto last = ampls.end() - 1; + while (hitStripsTrim > 1 && (*begin < std::max(trimCut, trimMaxFracNeigh_ * (*(begin + 1))))) { + hitStripsTrim--; + ++begin; + } + while (hitStripsTrim > 1 && (*last < std::max(trimCut, trimMaxFracNeigh_ * (*(last - 1))))) { + hitStripsTrim--; + --last; + } + if (hitStripsTrim < std::floor(std::abs(hitPredPos) - maxTrimmedSizeDiffNeg_)) { + filter_ = false; + } else if (hitStripsTrim <= std::ceil(std::abs(hitPredPos) + maxTrimmedSizeDiffPos_)) { + filter_ = true; + } else { + filter_ = peakFilter_; + } } diff --git a/DataFormats/SiStripCluster/src/SiStripCluster.cc b/DataFormats/SiStripCluster/src/SiStripCluster.cc index 3d02bf727648e..9deb73b4c4a22 100644 --- a/DataFormats/SiStripCluster/src/SiStripCluster.cc +++ b/DataFormats/SiStripCluster/src/SiStripCluster.cc @@ -26,6 +26,7 @@ SiStripCluster::SiStripCluster(const SiStripApproximateCluster cluster, const ui barycenter_ = cluster.barycenter() / 10.0; charge_ = cluster.width() * cluster.avgCharge(); amplitudes_.resize(cluster.width(), cluster.avgCharge()); + filter_ = cluster.filter(); float halfwidth_ = 0.5f * float(cluster.width()); @@ -60,3 +61,10 @@ float SiStripCluster::barycenter() const { // Need to mask off the high bit of firstStrip_, which contains the merged status. return float((firstStrip_ & stripIndexMask)) + float(sumx) / float(suma) + 0.5f; } +bool SiStripCluster::filter() const { + if (barycenter_ > 0) + return filter_; + return false; +} + +bool SiStripCluster::isFromApprox() const { return (barycenter_ > 0); } diff --git a/DataFormats/SiStripCluster/src/classes_def.xml b/DataFormats/SiStripCluster/src/classes_def.xml index 8c8c0a49d911a..3efcd26a23881 100755 --- a/DataFormats/SiStripCluster/src/classes_def.xml +++ b/DataFormats/SiStripCluster/src/classes_def.xml @@ -1,6 +1,7 @@ - + + @@ -24,7 +25,8 @@ - + + diff --git a/DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h b/DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h index 7f5d86d420854..a8c9a755f6cbb 100644 --- a/DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h +++ b/DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h @@ -44,6 +44,7 @@ namespace sistrip { static const uint16_t STRIPS_PER_FEDCH = STRIPS_PER_APV * APVS_PER_FEDCH; static const uint16_t STRIPS_PER_FEUNIT = STRIPS_PER_FEDCH * FEDCH_PER_FEUNIT; // 3072 static const uint16_t STRIPS_PER_FED = STRIPS_PER_FEUNIT * FEUNITS_PER_FED; // 24576 + static constexpr float MeVperADCStrip = 9.5665E-4; // -------------------- FED buffers -------------------- diff --git a/HLTrigger/Configuration/python/customizeHLTforCMSSW.py b/HLTrigger/Configuration/python/customizeHLTforCMSSW.py index cb80c17f2bb22..1cb644d600e12 100644 --- a/HLTrigger/Configuration/python/customizeHLTforCMSSW.py +++ b/HLTrigger/Configuration/python/customizeHLTforCMSSW.py @@ -226,6 +226,16 @@ def customizeHLTfor41495(process): return process +def customizeHLTfor41815(process): + # use hlt online BeamSpot for SiStripClusters2ApproxClusters + for producer in producers_by_type(process, 'SiStripClusters2ApproxClusters'): + producer.beamSpot = cms.InputTag('hltOnlineBeamSpot') + + if hasattr(process, 'HLT_HIRandom_v4'): + getattr(process,'HLT_HIRandom_v4').insert(2,process.HLTBeamSpot) + + return process + # CMSSW version specific customizations def customizeHLTforCMSSW(process, menuType="GRun"): @@ -236,5 +246,6 @@ def customizeHLTforCMSSW(process, menuType="GRun"): process = customizeHLTfor41058(process) process = customizeHLTfor41495(process) + process = customizeHLTfor41815(process) return process diff --git a/RecoLocalTracker/SiStripClusterizer/plugins/BuildFile.xml b/RecoLocalTracker/SiStripClusterizer/plugins/BuildFile.xml index 27ae390132063..8933e253d7be6 100644 --- a/RecoLocalTracker/SiStripClusterizer/plugins/BuildFile.xml +++ b/RecoLocalTracker/SiStripClusterizer/plugins/BuildFile.xml @@ -6,5 +6,8 @@ + + + diff --git a/RecoLocalTracker/SiStripClusterizer/plugins/SiStripApprox2ApproxClusters.cc b/RecoLocalTracker/SiStripClusterizer/plugins/SiStripApprox2ApproxClusters.cc index 910e8c8eb54a2..2c6d8f838cf78 100644 --- a/RecoLocalTracker/SiStripClusterizer/plugins/SiStripApprox2ApproxClusters.cc +++ b/RecoLocalTracker/SiStripClusterizer/plugins/SiStripApprox2ApproxClusters.cc @@ -59,6 +59,7 @@ void SiStripApprox2ApproxClusters::produce(edm::Event& event, edm::EventSetup co float barycenter = cluster.barycenter(); uint8_t width = cluster.width(); float avgCharge = cluster.avgCharge(); + bool filter = cluster.filter(); bool isSaturated = cluster.isSaturated(); switch (approxVersion) { @@ -86,7 +87,7 @@ void SiStripApprox2ApproxClusters::produce(edm::Event& event, edm::EventSetup co break; } - ff.push_back(SiStripApproximateCluster(barycenter, width, avgCharge, isSaturated)); + ff.push_back(SiStripApproximateCluster(barycenter, width, avgCharge, filter, isSaturated)); } } diff --git a/RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusters2ApproxClusters.cc b/RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusters2ApproxClusters.cc index 2f4d32f672aea..c842c384c4971 100644 --- a/RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusters2ApproxClusters.cc +++ b/RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusters2ApproxClusters.cc @@ -1,16 +1,33 @@ - - -#include "FWCore/Framework/interface/MakerMacros.h" +#include "CalibFormats/SiStripObjects/interface/SiStripDetInfo.h" +#include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h" +#include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h" +#include "CondFormats/SiStripObjects/interface/SiStripNoises.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h" +#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackBase.h" +#include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h" #include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/Framework/interface/stream/EDProducer.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/Utilities/interface/InputTag.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ParameterSet/interface/ParameterSetDescription.h" -#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h" -#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" -#include "DataFormats/Common/interface/DetSetVectorNew.h" -#include "DataFormats/Common/interface/DetSetVector.h" +#include "FWCore/Utilities/interface/ESInputTag.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h" +#include "RecoTracker/PixelLowPtUtilities/interface/SlidingPeakFinder.h" +#include "RecoTracker/Record/interface/CkfComponentsRecord.h" #include #include @@ -27,6 +44,25 @@ class SiStripClusters2ApproxClusters : public edm::stream::EDProducer<> { edm::EDGetTokenT > clusterToken; unsigned int maxNSat; + static constexpr double subclusterWindow_ = .7; + static constexpr double seedCutMIPs_ = .35; + static constexpr double seedCutSN_ = 7.; + static constexpr double subclusterCutMIPs_ = .45; + static constexpr double subclusterCutSN_ = 12.; + + edm::InputTag beamSpot_; + edm::EDGetTokenT beamSpotToken_; + + edm::ESGetToken tkGeomToken_; + + edm::FileInPath fileInPath_; + SiStripDetInfo detInfo_; + + std::string csfLabel_; + edm::ESGetToken csfToken_; + + edm::ESGetToken stripNoiseToken_; + edm::ESHandle theNoise_; }; SiStripClusters2ApproxClusters::SiStripClusters2ApproxClusters(const edm::ParameterSet& conf) { @@ -34,18 +70,79 @@ SiStripClusters2ApproxClusters::SiStripClusters2ApproxClusters(const edm::Parame maxNSat = conf.getParameter("maxSaturatedStrips"); clusterToken = consumes >(inputClusters); + + beamSpot_ = conf.getParameter("beamSpot"); + beamSpotToken_ = consumes(beamSpot_); + + tkGeomToken_ = esConsumes(); + + fileInPath_ = edm::FileInPath(SiStripDetInfoFileReader::kDefaultFile); + detInfo_ = SiStripDetInfoFileReader::read(fileInPath_.fullPath()); + + csfLabel_ = conf.getParameter("clusterShapeHitFilterLabel"); + csfToken_ = esConsumes(edm::ESInputTag("", csfLabel_)); + + stripNoiseToken_ = esConsumes(); + produces >(); } -void SiStripClusters2ApproxClusters::produce(edm::Event& event, edm::EventSetup const&) { +void SiStripClusters2ApproxClusters::produce(edm::Event& event, edm::EventSetup const& iSetup) { auto result = std::make_unique >(); const auto& clusterCollection = event.get(clusterToken); + auto const beamSpotHandle = event.getHandle(beamSpotToken_); + auto const& bs = beamSpotHandle.isValid() ? *beamSpotHandle : reco::BeamSpot(); + if (not beamSpotHandle.isValid()) { + edm::LogError("SiStripClusters2ApproxClusters") + << "didn't find a valid beamspot with label \"" << beamSpot_.encode() << "\" -> using (0,0,0)"; + } + + const auto& tkGeom = &iSetup.getData(tkGeomToken_); + const auto& theFilter = &iSetup.getData(csfToken_); + const auto& theNoise_ = &iSetup.getData(stripNoiseToken_); + for (const auto& detClusters : clusterCollection) { edmNew::DetSetVector::FastFiller ff{*result, detClusters.id()}; - for (const auto& cluster : detClusters) - ff.push_back(SiStripApproximateCluster(cluster, maxNSat)); + unsigned int detId = detClusters.id(); + const GeomDet* det = tkGeom->idToDet(detId); + double nApvs = detInfo_.getNumberOfApvsAndStripLength(detId).first; + double stripLength = detInfo_.getNumberOfApvsAndStripLength(detId).second; + double barycenter_ypos = 0.5 * stripLength; + + const StripGeomDetUnit* stripDet = dynamic_cast(det); + float mip = 3.9 / (sistrip::MeVperADCStrip / stripDet->surface().bounds().thickness()); + + for (const auto& cluster : detClusters) { + const LocalPoint& lp = LocalPoint(((cluster.barycenter() * 10 / (sistrip::STRIPS_PER_APV * nApvs)) - + ((stripDet->surface().bounds().width()) * 0.5f)), + barycenter_ypos - (0.5f * stripLength), + 0.); + const GlobalPoint& gpos = det->surface().toGlobal(lp); + GlobalPoint beamspot(bs.position().x(), bs.position().y(), bs.position().z()); + const GlobalVector& gdir = gpos - beamspot; + const LocalVector& ldir = det->toLocal(gdir); + + int hitStrips; + float hitPredPos; + theFilter->getSizes(detId, cluster, lp, ldir, hitStrips, hitPredPos); + + bool peakFilter = false; + SlidingPeakFinder pf(std::max(2, std::ceil(std::abs(hitPredPos) + subclusterWindow_))); + float mipnorm = mip / std::abs(ldir.z()); + PeakFinderTest test(mipnorm, + detId, + cluster.firstStrip(), + theNoise_, + seedCutMIPs_, + seedCutSN_, + subclusterCutMIPs_, + subclusterCutSN_); + peakFilter = pf.apply(cluster.amplitudes(), test); + + ff.push_back(SiStripApproximateCluster(cluster, maxNSat, hitPredPos, peakFilter)); + } } event.put(std::move(result)); @@ -55,6 +152,8 @@ void SiStripClusters2ApproxClusters::fillDescriptions(edm::ConfigurationDescript edm::ParameterSetDescription desc; desc.add("inputClusters", edm::InputTag("siStripClusters")); desc.add("maxSaturatedStrips", 3); + desc.add("clusterShapeHitFilterLabel", "ClusterShapeHitFilter"); // add CSF label + desc.add("beamSpot", edm::InputTag("offlineBeamSpot")); // add BeamSpot tag descriptions.add("SiStripClusters2ApproxClusters", desc); } diff --git a/RecoLocalTracker/SiStripClusterizer/python/SiStripClusters2ApproxClusters_cff.py b/RecoLocalTracker/SiStripClusterizer/python/SiStripClusters2ApproxClusters_cff.py index 2454d2f64292f..ee674d4c60a9f 100644 --- a/RecoLocalTracker/SiStripClusterizer/python/SiStripClusters2ApproxClusters_cff.py +++ b/RecoLocalTracker/SiStripClusterizer/python/SiStripClusters2ApproxClusters_cff.py @@ -1,5 +1,14 @@ from RecoLocalTracker.SiStripClusterizer.SiStripClusters2ApproxClusters_cfi import * from Configuration.ProcessModifiers.approxSiStripClusters_cff import approxSiStripClusters + +from RecoTracker.PixelLowPtUtilities.ClusterShapeHitFilterESProducer_cfi import ClusterShapeHitFilterESProducer as _ClusterShapeHitFilterESProducer +hltClusterShapeHitFilterESProducer = _ClusterShapeHitFilterESProducer.clone(ComponentName = 'hltClusterShapeHitFilterESProducer') + +from RecoVertex.BeamSpotProducer.BeamSpotOnline_cfi import onlineBeamSpotProducer +hltBeamSpotProducer = onlineBeamSpotProducer.clone(src = 'hltScalersRawToDigi') hltSiStripClusters2ApproxClusters = SiStripClusters2ApproxClusters.clone() -approxSiStripClusters.toModify(hltSiStripClusters2ApproxClusters, inputClusters = "siStripClustersHLT") +approxSiStripClusters.toModify(hltSiStripClusters2ApproxClusters, + beamSpot = "hltBeamSpotProducer", + inputClusters = "siStripClustersHLT", + clusterShapeHitFilterLabel = "hltClusterShapeHitFilterESProducer") diff --git a/RecoTracker/IterativeTracking/python/PixelLessStep_cff.py b/RecoTracker/IterativeTracking/python/PixelLessStep_cff.py index 5d6a4684d50f6..967a02c1a9ff8 100644 --- a/RecoTracker/IterativeTracking/python/PixelLessStep_cff.py +++ b/RecoTracker/IterativeTracking/python/PixelLessStep_cff.py @@ -201,15 +201,12 @@ FilterStripHits = cms.bool(True), ClusterShapeHitFilterName = cms.string('pixelLessStepClusterShapeHitFilter'), ClusterShapeCacheSrc = cms.InputTag('siPixelClusterShapeCache') # not really needed here since FilterPixelHits=False - ) + ), + _StripSubClusterShapeSeedFilter.clone() ) ) ) -from RecoTracker.PixelLowPtUtilities.StripSubClusterShapeSeedFilter_cfi import StripSubClusterShapeSeedFilter as _StripSubClusterShapeSeedFilter -from Configuration.ProcessModifiers.approxSiStripClusters_cff import approxSiStripClusters -(~approxSiStripClusters).toModify(pixelLessStepSeeds.SeedComparitorPSet.comparitors, func = lambda list: list.append(_StripSubClusterShapeSeedFilter.clone()) ) - trackingLowPU.toModify(pixelLessStepHitDoublets, produceSeedingHitSets=True, produceIntermediateHitDoublets=False) trackingLowPU.toModify(pixelLessStepSeeds, seedingHitSets = 'pixelLessStepHitDoublets', diff --git a/RecoTracker/IterativeTracking/python/TobTecStep_cff.py b/RecoTracker/IterativeTracking/python/TobTecStep_cff.py index 1df35f49af81c..97d7296501729 100644 --- a/RecoTracker/IterativeTracking/python/TobTecStep_cff.py +++ b/RecoTracker/IterativeTracking/python/TobTecStep_cff.py @@ -93,6 +93,7 @@ extraPhiKDBox = 0.01, ) from RecoTracker.TkSeedGenerator.seedCreatorFromRegionConsecutiveHitsEDProducer_cff import seedCreatorFromRegionConsecutiveHitsEDProducer as _seedCreatorFromRegionConsecutiveHitsTripletOnlyEDProducer +from RecoTracker.PixelLowPtUtilities.StripSubClusterShapeSeedFilter_cfi import StripSubClusterShapeSeedFilter as _StripSubClusterShapeSeedFilter _tobTecStepSeedComparitorPSet = dict( ComponentName = 'CombinedSeedComparitor', mode = cms.string('and'), @@ -104,7 +105,8 @@ FilterStripHits = cms.bool(True), ClusterShapeHitFilterName = cms.string('tobTecStepClusterShapeHitFilter'), ClusterShapeCacheSrc = cms.InputTag('siPixelClusterShapeCache') # not really needed here since FilterPixelHits=False - ) + ), + _StripSubClusterShapeSeedFilter.clone() ) ) @@ -113,10 +115,6 @@ SeedComparitorPSet = _tobTecStepSeedComparitorPSet, ) -from RecoTracker.PixelLowPtUtilities.StripSubClusterShapeSeedFilter_cfi import StripSubClusterShapeSeedFilter as _StripSubClusterShapeSeedFilter -from Configuration.ProcessModifiers.approxSiStripClusters_cff import approxSiStripClusters -(~approxSiStripClusters).toModify(tobTecStepSeedsTripl.SeedComparitorPSet.comparitors, func = lambda list: list.append(_StripSubClusterShapeSeedFilter.clone()) ) - #fastsim import FastSimulation.Tracking.TrajectorySeedProducer_cfi from FastSimulation.Tracking.SeedingMigration import _hitSetProducerToFactoryPSet diff --git a/RecoTracker/PixelLowPtUtilities/interface/SlidingPeakFinder.h b/RecoTracker/PixelLowPtUtilities/interface/SlidingPeakFinder.h new file mode 100644 index 0000000000000..40e3189b58f4b --- /dev/null +++ b/RecoTracker/PixelLowPtUtilities/interface/SlidingPeakFinder.h @@ -0,0 +1,102 @@ +#ifndef RecoTracker_PixelLowPtUtilities_SlidingPeakFinder_h +#define RecoTracker_PixelLowPtUtilities_SlidingPeakFinder_h + +#include +#include +#include +#include +#include "CondFormats/SiStripObjects/interface/SiStripNoises.h" +#include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h" + +class SlidingPeakFinder { +public: + SlidingPeakFinder(unsigned int size) : size_(size), half_((size + 1) / 2) {} + + template + bool apply(const uint8_t *x, + const uint8_t *begin, + const uint8_t *end, + const Test &test, + bool verbose = false, + int firststrip = 0) { + const uint8_t *ileft = (x != begin) ? std::min_element(x - 1, x + half_) : begin - 1; + const uint8_t *iright = ((x + size_) < end) ? std::min_element(x + half_, std::min(x + size_ + 1, end)) : end; + uint8_t left = (ileft < begin ? 0 : *ileft); + uint8_t right = (iright >= end ? 0 : *iright); + uint8_t center = *std::max_element(x, std::min(x + size_, end)); + uint8_t maxmin = std::max(left, right); + if (maxmin < center) { + bool ret = test(center, maxmin); + if (ret) { + ret = test(ileft, iright, begin, end); + } + return ret; + } else { + return false; + } + } + + template + bool apply(const V &ls, const Test &test, bool verbose = false, int firststrip = 0) { + const uint8_t *begin = &*ampls.begin(); + const uint8_t *end = &*ampls.end(); + for (const uint8_t *x = begin; x < end - (half_ - 1); ++x) { + if (apply(x, begin, end, test, verbose, firststrip)) { + return true; + } + } + return false; + } + +private: + unsigned int size_, half_; +}; + +struct PeakFinderTest { + PeakFinderTest(float mip, + uint32_t detid, + uint32_t firstStrip, + const SiStripNoises *theNoise, + float seedCutMIPs, + float seedCutSN, + float subclusterCutMIPs, + float subclusterCutSN) + : mip_(mip), + detid_(detid), + firstStrip_(firstStrip), + noiseObj_(theNoise), + noises_(theNoise->getRange(detid)), + subclusterCutMIPs_(subclusterCutMIPs), + sumCut_(subclusterCutMIPs_ * mip_), + subclusterCutSN2_(subclusterCutSN * subclusterCutSN) { + cut_ = std::min(seedCutMIPs * mip, seedCutSN * noiseObj_->getNoise(firstStrip + 1, noises_)); + } + + bool operator()(uint8_t max, uint8_t min) const { return max - min > cut_; } + bool operator()(const uint8_t *left, const uint8_t *right, const uint8_t *begin, const uint8_t *end) const { + int yleft = (left < begin ? 0 : *left); + int yright = (right >= end ? 0 : *right); + float sum = 0.0; + int maxval = 0; + float noise = 0; + for (const uint8_t *x = left + 1; x < right; ++x) { + int baseline = (yleft * int(right - x) + yright * int(x - left)) / int(right - left); + sum += int(*x) - baseline; + noise += std::pow(noiseObj_->getNoise(firstStrip_ + int(x - begin), noises_), 2); + maxval = std::max(maxval, int(*x) - baseline); + } + if (sum > sumCut_ && sum * sum > noise * subclusterCutSN2_) + return true; + return false; + } + +private: + float mip_; + unsigned int detid_; + int firstStrip_; + const SiStripNoises *noiseObj_; + SiStripNoises::Range noises_; + uint8_t cut_; + float subclusterCutMIPs_, sumCut_, subclusterCutSN2_; +}; +#endif diff --git a/RecoTracker/PixelLowPtUtilities/src/StripSubClusterShapeTrajectoryFilter.cc b/RecoTracker/PixelLowPtUtilities/src/StripSubClusterShapeTrajectoryFilter.cc index de333a956d6bf..362d15005167a 100644 --- a/RecoTracker/PixelLowPtUtilities/src/StripSubClusterShapeTrajectoryFilter.cc +++ b/RecoTracker/PixelLowPtUtilities/src/StripSubClusterShapeTrajectoryFilter.cc @@ -10,6 +10,7 @@ #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h" #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" +#include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h" #include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -184,13 +185,13 @@ StripSubClusterShapeFilterBase::StripSubClusterShapeFilterBase(const edm::Parame StripSubClusterShapeFilterBase::~StripSubClusterShapeFilterBase() { #if 0 - std::cout << "StripSubClusterShapeFilterBase " << label_ <<": called " << called_ << std::endl; - std::cout << "StripSubClusterShapeFilterBase " << label_ <<": saturated " << saturated_ << std::endl; - std::cout << "StripSubClusterShapeFilterBase " << label_ <<": test " << test_ << std::endl; - std::cout << "StripSubClusterShapeFilterBase " << label_ <<": failTooNarrow " << failTooNarrow_ << std::endl; - std::cout << "StripSubClusterShapeFilterBase " << label_ <<": passTrim " << passTrim_ << std::endl; - std::cout << "StripSubClusterShapeFilterBase " << label_ <<": passSC " << passSC_ << std::endl; - std::cout << "StripSubClusterShapeFilterBase " << label_ <<": failTooLarge " << failTooLarge_ << std::endl; + std::cout << "StripSubClusterShapeFilterBase " << label_ << ": called " << called_ << std::endl; + std::cout << "StripSubClusterShapeFilterBase " << label_ << ": saturated " << saturated_ << std::endl; + std::cout << "StripSubClusterShapeFilterBase " << label_ << ": test " << test_ << std::endl; + std::cout << "StripSubClusterShapeFilterBase " << label_ << ": failTooNarrow " << failTooNarrow_ << std::endl; + std::cout << "StripSubClusterShapeFilterBase " << label_ << ": passTrim " << passTrim_ << std::endl; + std::cout << "StripSubClusterShapeFilterBase " << label_ << ": passSC " << passSC_ << std::endl; + std::cout << "StripSubClusterShapeFilterBase " << label_ << ": failTooLarge " << failTooLarge_ << std::endl; #endif } @@ -265,74 +266,77 @@ bool StripSubClusterShapeFilterBase::testLastHit(const TrackingRecHit *hit, return true; } - // compute number of consecutive saturated strips - // (i.e. with adc count >= 254, see SiStripCluster class for documentation) - unsigned int thisSat = (ampls[0] >= 254), maxSat = thisSat; - for (unsigned int i = 1, n = ampls.size(); i < n; ++i) { - if (ampls[i] >= 254) { - thisSat++; - } else if (thisSat > 0) { + if (!cluster.isFromApprox()) { + // compute number of consecutive saturated strips + // (i.e. with adc count >= 254, see SiStripCluster class for documentation) + unsigned int thisSat = (ampls[0] >= 254), maxSat = thisSat; + for (unsigned int i = 1, n = ampls.size(); i < n; ++i) { + if (ampls[i] >= 254) { + thisSat++; + } else if (thisSat > 0) { + maxSat = std::max(maxSat, thisSat); + thisSat = 0; + } + } + if (thisSat > 0) { maxSat = std::max(maxSat, thisSat); - thisSat = 0; } - } - if (thisSat > 0) { - maxSat = std::max(maxSat, thisSat); - } - if (maxSat >= maxNSat_) { - INC_COUNTER(saturated_) - return true; - } + if (maxSat >= maxNSat_) { + INC_COUNTER(saturated_) + return true; + } - // trimming - INC_COUNTER(test_) - unsigned int hitStripsTrim = ampls.size(); - int sum = std::accumulate(ampls.begin(), ampls.end(), 0); - uint8_t trimCut = std::min(trimMaxADC_, std::floor(trimMaxFracTotal_ * sum)); - auto begin = ampls.begin(); - auto last = ampls.end() - 1; - while (hitStripsTrim > 1 && (*begin < std::max(trimCut, trimMaxFracNeigh_ * (*(begin + 1))))) { - hitStripsTrim--; - ++begin; - } - while (hitStripsTrim > 1 && (*last < std::max(trimCut, trimMaxFracNeigh_ * (*(last - 1))))) { - hitStripsTrim--; - --last; - } + // trimming + INC_COUNTER(test_) + unsigned int hitStripsTrim = ampls.size(); + int sum = std::accumulate(ampls.begin(), ampls.end(), 0); + uint8_t trimCut = std::min(trimMaxADC_, std::floor(trimMaxFracTotal_ * sum)); + auto begin = ampls.begin(); + auto last = ampls.end() - 1; + while (hitStripsTrim > 1 && (*begin < std::max(trimCut, trimMaxFracNeigh_ * (*(begin + 1))))) { + hitStripsTrim--; + ++begin; + } + while (hitStripsTrim > 1 && (*last < std::max(trimCut, trimMaxFracNeigh_ * (*(last - 1))))) { + hitStripsTrim--; + --last; + } - if (hitStripsTrim < std::floor(std::abs(hitPredPos) - maxTrimmedSizeDiffNeg_)) { - INC_COUNTER(failTooNarrow_) - return false; - } else if (hitStripsTrim <= std::ceil(std::abs(hitPredPos) + maxTrimmedSizeDiffPos_)) { - INC_COUNTER(passTrim_) - return true; - } + if (hitStripsTrim < std::floor(std::abs(hitPredPos) - maxTrimmedSizeDiffNeg_)) { + INC_COUNTER(failTooNarrow_) + return false; + } else if (hitStripsTrim <= std::ceil(std::abs(hitPredPos) + maxTrimmedSizeDiffPos_)) { + INC_COUNTER(passTrim_) + return true; + } - const StripGeomDetUnit *stripDetUnit = dynamic_cast(det); - if (det == nullptr) { - edm::LogError("Strip not a StripGeomDetUnit?") << " on " << detId.rawId() << "\n"; - return true; - } + const StripGeomDetUnit *stripDetUnit = dynamic_cast(det); + if (det == nullptr) { + edm::LogError("Strip not a StripGeomDetUnit?") << " on " << detId.rawId() << "\n"; + return true; + } - float MeVperADCStrip = 9.5665E-4; // conversion constant from ADC counts to MeV for the SiStrip detector - float mip = - 3.9 / (MeVperADCStrip / stripDetUnit->surface().bounds().thickness()); // 3.9 MeV/cm = ionization in silicon - float mipnorm = mip / std::abs(ldir.z()); - ::SlidingPeakFinder pf(std::max(2, std::ceil(std::abs(hitPredPos) + subclusterWindow_))); - ::PeakFinderTest test(mipnorm, - detId(), - cluster.firstStrip(), - &*theNoise, - seedCutMIPs_, - seedCutSN_, - subclusterCutMIPs_, - subclusterCutSN_); - if (pf.apply(cluster.amplitudes(), test)) { - INC_COUNTER(passSC_) - return true; + float mip = 3.9 / (sistrip::MeVperADCStrip / + stripDetUnit->surface().bounds().thickness()); // 3.9 MeV/cm = ionization in silicon + float mipnorm = mip / std::abs(ldir.z()); + ::SlidingPeakFinder pf(std::max(2, std::ceil(std::abs(hitPredPos) + subclusterWindow_))); + ::PeakFinderTest test(mipnorm, + detId(), + cluster.firstStrip(), + &*theNoise, + seedCutMIPs_, + seedCutSN_, + subclusterCutMIPs_, + subclusterCutSN_); + if (pf.apply(cluster.amplitudes(), test)) { + INC_COUNTER(passSC_) + return true; + } else { + INC_COUNTER(failTooLarge_) + return false; + } } else { - INC_COUNTER(failTooLarge_) - return false; + return cluster.filter(); } } return true;