37 #ifdef StripSubClusterShapeFilterBase_COUNTERS
38 #define INC_COUNTER(X) X++;
40 #define INC_COUNTER(X)
43 class SlidingPeakFinder {
45 SlidingPeakFinder(
unsigned int size) :
46 size_(size), half_((size+1)/2) {}
48 template<
typename Test>
49 bool apply(
const uint8_t *
x,
const uint8_t *
begin,
const uint8_t *
end,
const Test &
test,
bool verbose=
false,
int firststrip=0) {
50 const uint8_t * ileft = (x !=
begin) ? std::min_element(x-1,x+half_) : begin - 1;
52 uint8_t left = (ileft < begin ? 0 : *ileft );
53 uint8_t right = (iright >= end ? 0 : *iright);
55 uint8_t maxmin =
std::max(left,right);
56 if (maxmin < center) {
57 bool ret =
test(center, maxmin);
59 ret =
test(ileft,iright,begin,end);
67 template<
typename V,
typename Test>
68 bool apply(
const V & ampls,
const Test & test,
bool verbose=
false,
int firststrip=0) {
69 const uint8_t *begin = &*ampls.begin();
70 const uint8_t *end = &*ampls.end();
71 for (
const uint8_t *x = begin; x < end - (half_-1); ++
x) {
80 unsigned int size_, half_;
84 struct PeakFinderTest {
85 PeakFinderTest(
float mip, uint32_t
detid, uint32_t firstStrip,
const SiStripNoises *theNoise,
86 float seedCutMIPs,
float seedCutSN,
87 float subclusterCutMIPs,
float subclusterCutSN) :
88 mip_(mip), detid_(detid), firstStrip_(firstStrip), noiseObj_(theNoise), noises_(theNoise->
getRange(detid)),
89 subclusterCutMIPs_(subclusterCutMIPs), sumCut_(subclusterCutMIPs_*mip_), subclusterCutSN2_(subclusterCutSN*subclusterCutSN)
91 cut_ = std::min<float>(seedCutMIPs*mip, seedCutSN*noiseObj_->getNoise(firstStrip+1, noises_));
94 bool operator()(uint8_t
max, uint8_t
min)
const {
95 return max-min > cut_;
97 bool operator()(
const uint8_t *left,
const uint8_t *right,
const uint8_t *
begin,
const uint8_t *
end)
const {
98 int yleft = (left < begin ? 0 : *left);
99 int yright = (right >= end ? 0 : *right);
101 int maxval = 0;
float noise = 0;
102 for (
const uint8_t *
x = left+1;
x < right; ++
x) {
103 int baseline = (yleft * int(right-
x) + yright * int(
x-left)) /
int(right-left);
104 sum += int(*
x) - baseline;
105 noise +=
std::pow(noiseObj_->getNoise(firstStrip_ +
int(
x-begin), noises_), 2);
106 maxval =
std::max(maxval,
int(*
x) - baseline);
108 if (sum > sumCut_ && sum*sum > noise*subclusterCutSN2_)
return true;
114 unsigned int detid_;
int firstStrip_;
118 float subclusterCutMIPs_, sumCut_, subclusterCutSN2_;
129 trimMaxFracTotal_(iCfg.
getParameter<
double>(
"trimMaxFracTotal")),
130 trimMaxFracNeigh_(iCfg.
getParameter<
double>(
"trimMaxFracNeigh")),
131 maxTrimmedSizeDiffPos_(iCfg.
getParameter<
double>(
"maxTrimmedSizeDiffPos")),
132 maxTrimmedSizeDiffNeg_(iCfg.
getParameter<
double>(
"maxTrimmedSizeDiffNeg")),
133 subclusterWindow_(iCfg.
getParameter<
double>(
"subclusterWindow")),
136 subclusterCutMIPs_(iCfg.
getParameter<
double>(
"subclusterCutMIPs")),
137 subclusterCutSN_(iCfg.
getParameter<
double>(
"subclusterCutSN"))
138 #ifdef StripSubClusterShapeFilterBase_COUNTERS
139 ,called_(0),saturated_(0),test_(0),passTrim_(0),failTooLarge_(0),passSC_(0),failTooNarrow_(0)
142 if (iCfg.
exists(
"layerMask")) {
144 const char *ndets[4] = {
"TIB",
"TID",
"TOB",
"TEC" };
145 const int idets[4] = { 3, 4, 5, 6 };
146 for (
unsigned int i = 0;
i < 4; ++
i) {
150 layerMask_[idets[
i]][0] = 2;
152 for (uint32_t lay : iLM.
getParameter<std::vector<uint32_t>>(ndets[i])) {
153 layerMask_[idets[
i]][lay] = 1;
158 for (
auto & arr : layerMask_) {
168 std::cout <<
"StripSubClusterShapeFilterBase " <<
label_ <<
": called " << called_ << std::endl;
169 std::cout <<
"StripSubClusterShapeFilterBase " <<
label_ <<
": saturated " << saturated_ << std::endl;
170 std::cout <<
"StripSubClusterShapeFilterBase " <<
label_ <<
": test " << test_ << std::endl;
171 std::cout <<
"StripSubClusterShapeFilterBase " <<
label_ <<
": failTooNarrow " << failTooNarrow_ << std::endl;
172 std::cout <<
"StripSubClusterShapeFilterBase " <<
label_ <<
": passTrim " << passTrim_ << std::endl;
173 std::cout <<
"StripSubClusterShapeFilterBase " <<
label_ <<
": passSC " << passSC_ << std::endl;
174 std::cout <<
"StripSubClusterShapeFilterBase " <<
label_ <<
": failTooLarge " << failTooLarge_ << std::endl;
193 return testLastHit(&mono, gpos, gdir,
true) && testLastHit(&stereo, gpos, gdir,
true);
197 return testLastHit(&orig, gpos, gdir,
true);
198 }
else if ((stripHit = dynamic_cast<const TrackerSingleRecHit *>(hit)) != 0) {
201 if (layerMask_[detId.
subdetId()][0] == 0) {
203 }
else if (layerMask_[detId.
subdetId()][0] == 2) {
204 unsigned int ilayer = theTopology->layer(detId);
205 if (layerMask_[detId.
subdetId()][ilayer] == 0) {
210 const GeomDet *det = theTracker->idToDet(detId);
214 lpos -= ldir * lpos.
z()/ldir.
z();
216 int hitStrips;
float hitPredPos;
218 bool usable = theFilter->getSizes(detId, cluster, lpos, ldir, hitStrips, hitPredPos);
219 if (!usable)
return true;
225 if (
std::abs(hitPredPos) < 1.5
f && hitStrips <= 2) {
232 unsigned int thisSat = (ampls[0] >= 254), maxSat = thisSat;
233 for (
unsigned int i = 1,
n = ampls.size();
i <
n; ++
i) {
234 if (ampls[
i] >= 254) {
236 }
else if (thisSat > 0) {
237 maxSat = std::max<int>(maxSat, thisSat);
242 maxSat = std::max<int>(maxSat, thisSat);
244 if (maxSat >= maxNSat_) {
251 unsigned int hitStripsTrim = ampls.size();
252 int sum = std::accumulate(ampls.begin(), ampls.end(), 0);
253 uint8_t trimCut = std::min<uint8_t>(trimMaxADC_, std::floor(trimMaxFracTotal_ * sum));
254 auto begin = ampls.begin();
255 auto last = ampls.end()-1;
256 while (hitStripsTrim > 1 && (*begin < std::max<uint8_t>(trimCut, trimMaxFracNeigh_*(*(begin+1)))) ) { hitStripsTrim--; ++
begin; }
257 while (hitStripsTrim > 1 && (*
last < std::max<uint8_t>(trimCut, trimMaxFracNeigh_*(*(
last -1)))) ) { hitStripsTrim--; --
last; }
259 if (hitStripsTrim < std::floor(
std::abs(hitPredPos)-maxTrimmedSizeDiffNeg_)) {
262 }
else if (hitStripsTrim <= std::ceil(
std::abs(hitPredPos)+maxTrimmedSizeDiffPos_)) {
268 if (det == 0) {
edm::LogError(
"Strip not a StripGeomDetUnit?") <<
" on " << detId.
rawId() <<
"\n";
return true; }
273 ::SlidingPeakFinder pf(std::max<int>(2,std::ceil(
std::abs(hitPredPos)+subclusterWindow_)));
274 ::PeakFinderTest
test(mipnorm, detId(), cluster.
firstStrip(), &*theNoise, seedCutMIPs_, seedCutSN_, subclusterCutMIPs_, subclusterCutSN_);
306 throw cms::Exception(
"toBeContinued(Traj) instead of toBeContinued(TempTraj)");
314 if (hit == 0 || !hit->
isValid())
return true;
325 return testLastHit(*tms.
rbegin());
333 return testLastHit(*tms.rbegin());
341 return testLastHit(*tms.
rbegin());
349 filterAtHelixStage_(iConfig.
getParameter<
bool>(
"FilterAtHelixStage"))
352 if (filterAtHelixStage_)
edm::LogError(
"Configuration") <<
"StripSubClusterShapeSeedFilter: FilterAtHelixStage is not yet working correctly.\n";
360 if (filterAtHelixStage_)
return true;
362 if (hit == 0 || !hit->
isValid())
return true;
364 return testLastHit(hit, tsos,
false);
370 if (!filterAtHelixStage_)
return true;
372 if (!helix.
isValid())
edm::LogWarning(
"InvalidHelix") <<
"PixelClusterShapeSeedComparitor helix is not valid, result is bad";
378 float x0 = vertex.
x(), y0 = vertex.
y();
379 for (
unsigned int i = 0,
n = hits.
size();
i <
n; ++
i) {
380 auto const & hit = *hits[
i];
384 float x1 = pos.
x(), y1 = pos.
y(), dx1 = x1 - xc, dy1 = y1 - yc;
387 float perpx = -dy1, perpy = dx1;
388 if (perpx * (x1-x0) + perpy * (y1 - y0) < 0) {
389 perpy = -perpy; perpx = -perpx;
393 float pnorm = 1.0/
std::sqrt(perpx*perpx + perpy*perpy + 1);
394 GlobalVector gdir(perpx*pnorm, perpy*pnorm, (momvtx.
z() > 0 ? pnorm : -pnorm));
396 if (!testLastHit(&hit, pos, gdir)) {
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
virtual bool qualityFilter(const TempTrajectory &) const override
bool existsAs(std::string const ¶meterName, bool trackiness=true) const
checks if a parameter exists as a given type
virtual bool toBeContinued(TempTrajectory &) const override
ConstRecHitPointer const & recHit() const
const DataContainer & measurements() const
bool exists(std::string const ¶meterName) const
checks if a parameter exists
GlobalPoint globalPosition() const
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
const Bounds & bounds() const
uint16_t firstStrip() const
const FastCircle & circle() const
virtual GlobalPoint globalPosition() const
const Plane & surface() const
The nominal surface of the GeomDet.
uint32_t rawId() const
get the raw id
virtual float thickness() const =0
virtual ~StripSubClusterShapeFilterBase()
DataContainer const & measurements() const
BaseTrackerRecHit const * ConstRecHitPointer
std::vector< TrajectoryMeasurement > DataContainer
GlobalVector momentum() const
Abs< T >::type abs(const T &t)
bool testLastHit(const TrackingRecHit *hit, const TrajectoryStateOnSurface &tsos, bool mustProject=false) const
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
StripSubClusterShapeFilterBase(const edm::ParameterSet &iConfig, edm::ConsumesCollector &iC)
GlobalPoint position() const
const_iterator rbegin() const
void setEventBase(const edm::Event &, const edm::EventSetup &)
SiStripRecHit2D originalHit() const
virtual TrackingRecHit const * hit() const
SiStripRecHit2D stereoHit() const
SiStripCluster const & stripCluster() const
virtual bool compatible(const TrajectoryStateOnSurface &tsos, SeedingHitSet::ConstRecHitPointer hit) const
SiStripRecHit2D monoHit() const
const Range getRange(const uint32_t detID) const
unsigned int size() const
TrajectoryStateOnSurface const & updatedState() const
DetId geographicalId() const
bool testLastHit(const TrajectoryMeasurement &last) const
std::pair< ContainerIterator, ContainerIterator > Range
tuple size
Write out results.
const std::vector< uint8_t > & amplitudes() const
Power< A, B >::type pow(const A &a, const B &b)
GlobalVector globalDirection() const
StripSubClusterShapeSeedFilter(const edm::ParameterSet &iConfig, edm::ConsumesCollector &iC)