20 return layer[0].name() +
"+" + layer[1].name() +
"+" + layer[2].name();
28 writeTriplets_(conf.getParameter<
bool>(
"writeTriplets")),
29 seedOnMiddle_(conf.existsAs<
bool>(
"seedOnMiddle") ? conf.getParameter<
bool>(
"seedOnMiddle") :
false),
30 rescaleError_(conf.existsAs<double>(
"rescaleError") ? conf.getParameter<double>(
"rescaleError") : 1.0),
31 tripletsVerbosity_(conf.getUntrackedParameter<uint32_t>(
"TripletsDebugLevel", 0)),
32 seedVerbosity_(conf.getUntrackedParameter<uint32_t>(
"seedDebugLevel", 0)),
33 helixVerbosity_(conf.getUntrackedParameter<uint32_t>(
"helixDebugLevel", 0)),
34 check_(conf.getParameter<
edm::
ParameterSet>(
"ClusterCheckPSet"), consumesCollector()),
35 maxTriplets_(conf.getParameter<int32_t>(
"maxTriplets")),
36 maxSeeds_(conf.getParameter<int32_t>(
"maxSeeds")) {
39 float originradius = regionConf.
getParameter<
double>(
"originRadius");
40 float halflength = regionConf.
getParameter<
double>(
"originHalfLength");
41 float originz = regionConf.
getParameter<
double>(
"originZPosition");
52 produces<TrajectorySeedCollection>();
54 produces<edm::OwnVector<TrackingRecHit>>(
"cosmicTriplets");
90 auto output = std::make_unique<TrajectorySeedCollection>();
91 auto outtriplets = std::make_unique<edm::OwnVector<TrackingRecHit>>();
97 edm::LogError(
"TooManyClusters") <<
"Found too many clusters (" << clustsOrZero <<
"), bailing out.\n";
111 outtriplets->push_back(hit1->
clone());
112 outtriplets->push_back(hit2->
clone());
113 outtriplets->push_back(hit3->
clone());
152 float iy1 = ihit1->globalPosition().y();
153 float oy1 = ohit1->globalPosition().y();
154 float iy2 = ihit2->globalPosition().y();
155 float oy2 = ohit2->globalPosition().y();
159 return (iy1 != iy2 ? (iy1 > iy2) : (oy1 > oy2));
162 }
else if (oy2 - iy2 > 0) {
166 return (iy1 != iy2 ? (iy1 < iy2) : (oy1 < oy2));
191 std::cout <<
"GenericTripletGenerator iLss = " << seedingLayersToString(ls) <<
" (" << layerIndex
192 <<
"): # = " << innerHits.size() <<
"/" << middleHits.size() <<
"/" << outerHits.size() << std::endl;
197 std::vector<TTRH> innerTTRHs, middleTTRHs, outerTTRHs;
200 std::vector<bool> innerOk(innerHits.size(),
true);
201 std::vector<bool> middleOk(middleHits.size(),
true);
202 std::vector<bool> outerOk(outerHits.size(),
true);
207 for (
auto iOuterHit = outerHits.begin(); iOuterHit != outerHits.end(); ++
idx, ++iOuterHit) {
208 outerTTRHs.push_back(&(**iOuterHit));
210 outerOk[
idx] =
false;
213 for (
auto iMiddleHit = middleHits.begin(); iMiddleHit != middleHits.end(); ++
idx, ++iMiddleHit) {
214 middleTTRHs.push_back(&(**iMiddleHit));
216 middleOk[
idx] =
false;
219 for (
auto iInnerHit = innerHits.begin(); iInnerHit != innerHits.end(); ++
idx, ++iInnerHit) {
220 innerTTRHs.push_back(&(**iInnerHit));
222 innerOk[
idx] =
false;
230 for (
auto iOuterHit = outerHits.begin(); iOuterHit != outerHits.end(); iOuterHit++) {
231 idx = iOuterHit - outerHits.begin();
232 TTRH &outerTTRH = outerTTRHs[
idx];
233 GlobalPoint outerpos = outerTTRH->globalPosition();
234 bool outerok = outerOk[
idx];
241 for (
auto iMiddleHit = middleHits.begin(); iMiddleHit != middleHits.end(); iMiddleHit++) {
242 idx = iMiddleHit - middleHits.begin();
243 TTRH &middleTTRH = middleTTRHs[
idx];
244 GlobalPoint middlepos = middleTTRH->globalPosition();
245 bool middleok = middleOk[
idx];
253 for (
auto iInnerHit = innerHits.begin(); iInnerHit != innerHits.end(); iInnerHit++) {
254 idx = iInnerHit - innerHits.begin();
255 TTRH &innerTTRH = innerTTRHs[
idx];
256 GlobalPoint innerpos = innerTTRH->globalPosition();
257 bool innerok = innerOk[
idx];
267 (innerpos.
y() < 0 || middlepos.
y() < 0 || outerpos.
y() < 0 || outerpos.
y() < innerpos.
y()))
270 (innerpos.
y() > 0 || middlepos.
y() > 0 || outerpos.
y() > 0 || outerpos.
y() > innerpos.
y()))
275 std::cout <<
"Trying seed with: " << innerpos <<
" + " << middlepos <<
" + " << outerpos << std::endl;
276 if (
goodTriplet(innerpos, middlepos, outerpos, minRho)) {
282 edm::LogError(
"TooManyTriplets") <<
"Found too many triplets, bailing out.\n";
287 <<
" + " << outerpos << std::endl;
291 <<
" + " << outerpos << std::endl;
301 std::cout <<
" iLss = " << seedingLayersToString(ls) <<
" (" << layerIndex
302 <<
"): # = " << innerHits.size() <<
"/" << middleHits.size() <<
"/" << outerHits.size() <<
": Found " 326 return checkCharge(static_cast<const SiStripRecHit2D &>(*hit), subdet);
348 std::vector<bool> &oks)
const {
350 std::vector<TTRH>::const_iterator it = hits.begin(),
start = it,
end = hits.end();
351 std::vector<bool>::iterator
ok = oks.begin(), okStart =
ok;
352 while (
start < end) {
353 DetId lastid = (*start)->geographicalId();
354 for (it =
start + 1; (it <
end) && ((*it)->geographicalId() == lastid); ++it) {
359 std::cerr <<
"SimpleCosmicBONSeeder: Marking noisy module " << lastid.
rawId() <<
", it has " << (it -
start)
366 std::cerr <<
"SimpleCosmicBONSeeder: Not marking noisy module " << lastid.
rawId() <<
", it has " << (it -
start)
379 const double &minRho)
const {
380 float dyOM = outer.
y() - middle.
y(), dyIM = inner.
y() - middle.
y();
381 if ((dyOM * dyIM > 0) && (fabs(dyOM) > 10) && (fabs(dyIM) > 10)) {
383 std::cout <<
" fail for non coherent dy" << std::endl;
386 float dzOM = outer.
z() - middle.
z(), dzIM = inner.
z() - middle.
z();
387 if ((dzOM * dzIM > 0) && (fabs(dzOM) > 50) && (fabs(dzIM) > 50)) {
389 std::cout <<
" fail for non coherent dz" << std::endl;
394 if (theCircle.
rho() < minRho) {
396 std::cout <<
" fail for pt cut" << std::endl;
408 std::cout <<
"DEBUG PZ =====" << std::endl;
411 std::cout <<
"FastHelix P = " << gv <<
"\n";
412 std::cout <<
"FastHelix Q = " << helix.stateAtVertex().charge() <<
"\n";
418 double rho = theCircle.
rho();
422 double pt = 0.01 * rho * (0.3 * tesla.
z());
425 double dx1 = outer.
x() - theCircle.
x0();
426 double dy1 = outer.
y() - theCircle.
y0();
427 double py = pt * dx1 /
rho,
px = -pt * dy1 /
rho;
428 if (px * (middle.
x() - outer.
x()) + py * (middle.
y() - outer.
y()) < 0.) {
434 double dz = inner.
z() - outer.
z();
435 double sinphi = (dx1 * (inner.
y() - theCircle.
y0()) - dy1 * (inner.
x() - theCircle.
x0())) / (rho *
rho);
436 double dphi =
std::abs(std::asin(sinphi));
437 double pz = pt * dz / (dphi *
rho);
439 int myq = ((theCircle.
x0() * py - theCircle.
y0() *
px) / tesla.
z()) > 0. ? +1 : -1;
441 std::pair<GlobalVector, int> mypq(
GlobalVector(px, py, pz), myq);
444 std::cout <<
"Gio: pt = " << pt << std::endl;
445 std::cout <<
"Gio: dz = " << dz <<
", sinphi = " << sinphi <<
", dphi = " << dphi
446 <<
", dz/drphi = " << (dz / dphi /
rho) << std::endl;
449 std::cout <<
"Gio's fit P = " << mypq.first <<
"\n";
450 std::cout <<
"Gio's fit Q = " << myq <<
"\n";
472 std::cout <<
"Processing triplet " << it <<
": " << inner <<
" + " << middle <<
" + " << outer << std::endl;
474 if ((outer.
y() - inner.
y()) * outer.
y() < 0) {
479 std::cout <<
"The seed was going away from CMS! swapped in <-> out" << std::endl;
480 std::cout <<
"Processing swapped triplet " << it <<
": " << inner <<
" + " << middle <<
" + " << outer
486 std::pair<GlobalVector, int> pq =
pqFromHelixFit(inner, middle, outer, iSetup);
488 float ch = pq.second;
489 float Mom =
sqrt(gv.
x() * gv.
x() + gv.
y() * gv.
y() + gv.
z() * gv.
z());
493 std::cout <<
"Processing triplet " << it <<
": fail for momentum." << std::endl;
505 if ((outer.
y() - inner.
y()) > 0) {
507 std::cout <<
"Processing triplet " << it <<
": downgoing." << std::endl;
514 std::cout <<
"Processing triplet " << it <<
": upgoing." << std::endl;
518 if ((gv.
z() * (outer.
z() - inner.
z()) > 0) && (fabs(outer.
z() - inner.
z()) > 5) && (fabs(gv.
z()) > .01)) {
519 std::cout <<
"ORRORE: outer.z()-inner.z() = " << (outer.
z() - inner.
z()) <<
", gv.z() = " << gv.
z()
528 std::cout <<
"Processing triplet " << it <<
". start from " << std::endl;
529 std::cout <<
" X = " << outer <<
", P = " << gv << std::endl;
535 TSOS propagated, updated;
537 for (
size_t ih = 0; ih < 3; ++ih) {
540 std::cout <<
"Stopping at middle hit, as requested." << std::endl;
544 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
"." << std::endl;
550 if (!propagated.isValid()) {
552 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
": failed propagation." << std::endl;
557 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
": propagated state = " << propagated;
563 if (!updated.isValid()) {
565 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
": failed update." << std::endl;
570 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
": updated state = " << updated;
573 if (!fail && updated.isValid() && (updated.globalMomentum().perp() <
region_.
ptMin())) {
575 std::cout <<
"Processing triplet " << it <<
": failed for final pt " << updated.globalMomentum().perp() <<
" < " 579 if (!fail && updated.isValid() && (updated.globalMomentum().mag() <
pMin_)) {
581 std::cout <<
"Processing triplet " << it <<
": failed for final p " << updated.globalMomentum().perp() <<
" < " 582 <<
pMin_ << std::endl;
589 <<
": state BEFORE rescaling " << updated;
590 std::cout <<
" Cartesian error (X,P) before rescaling= \n" 591 << updated.cartesianError().matrix() << std::endl;
596 std::cout <<
"Processed triplet " << it <<
": success (saved as #" << output.size() <<
") : " << inner <<
" + " 597 << middle <<
" + " << outer << std::endl;
598 std::cout <<
" pt = " << updated.globalMomentum().perp() <<
" eta = " << updated.globalMomentum().eta()
599 <<
" phi = " << updated.globalMomentum().phi() <<
" ch = " << updated.charge() << std::endl;
603 std::cout <<
" X = " << updated.globalPosition() <<
", P = " << updated.globalMomentum() << std::endl;
605 std::cout <<
" Cartesian error (X,P) = \n" << updated.cartesianError().matrix() << std::endl;
613 edm::LogError(
"TooManySeeds") <<
"Found too many seeds, bailing out.\n";
T getParameter(std::string const &) const
static constexpr auto TEC
GlobalTrackingRegion region_
unsigned short numberOfLayersInSet() const
Get number of layers in each SeedingLayerSets.
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
bool triplets(const edm::Event &e, const edm::EventSetup &c)
void init(const edm::EventSetup &c)
CartesianTrajectoryError cartesianError() const
bool existsAs(std::string const ¶meterName, bool trackiness=true) const
checks if a parameter exists as a given type
bool checkMaxHitsPerModule_
MiddleRecHit middle() const
int nominalValue() const
The nominal field value for this map in kGauss.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
constexpr bool isNotFinite(T x)
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
unsigned short LayerSetIndex
OrderedHitTriplets hitTriplets
Global3DPoint GlobalPoint
constexpr uint32_t rawId() const
get the raw id
uint16_t firstStrip() const
bool checkCharge(const TrackingRecHit *hit) const
SimpleCosmicBONSeeder(const edm::ParameterSet &conf)
InnerRecHit inner() const
edm::ESHandle< TrackerGeometry > tracker
bool matchedRecHitUsesAnd_
std::pair< GlobalVector, int > pqFromHelixFit(const GlobalPoint &inner, const GlobalPoint &middle, const GlobalPoint &outer, const edm::EventSetup &iSetup) const
TrajectoryStateOnSurface update(const TrajectoryStateOnSurface &, const TrackingRecHit &) const override
BaseTrackerRecHit const * ConstRecHitPointer
std::vector< TrajectorySeed > TrajectorySeedCollection
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
PropagatorWithMaterial * thePropagatorOp
edm::ESHandle< MagneticField > magfield
virtual TrackingRecHit * clone() const =0
size_t tooManyClusters(const edm::Event &e) const
edm::ESHandle< TransientTrackingRecHitBuilder > TTTRHBuilder
uint32_t tripletsVerbosity_
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
Abs< T >::type abs(const T &t)
edm::EDGetTokenT< SeedingLayerSetsHits > seedingLayerToken_
ClusterRef cluster() const
T const * get() const
Returns C++ pointer to the item.
static constexpr auto TOB
TrackingRegion::Hits hits(const edm::EventSetup &es, const SeedingLayerSetsHits::SeedingLayer &layer) const override
get hits from layer compatible with region constraints
unsigned int size() const override
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
PropagatorWithMaterial * thePropagatorAl
const AlgebraicSymMatrix66 & matrix() const
TrajectoryStateOnSurface TSOS
virtual TrackingRecHit const * hit() const
void rescaleError(double factor)
static constexpr auto TIB
SiStripRecHit2D stereoHit() const
bool seeds(TrajectorySeedCollection &output, const edm::EventSetup &iSetup)
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
float ptMin() const
minimal pt of interest
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
bool goodTriplet(const GlobalPoint &inner, const GlobalPoint &middle, const GlobalPoint &outer, const double &minRho) const
SiStripRecHit2D monoHit() const
bool operator()(const OrderedHitTriplet &trip1, const OrderedHitTriplet &trip2) const
void produce(edm::Event &e, const edm::EventSetup &c) override
const TrackerGeomDet * idToDet(DetId) const override
SeedingHitSet::ConstRecHitPointer SeedingHit
std::vector< int32_t > chargeThresholds_
void checkNoisyModules(const std::vector< SeedingHitSet::ConstRecHitPointer > &hits, std::vector< bool > &oks) const
std::vector< int32_t > maxHitsPerModule_
DetId geographicalId() const
OuterRecHit outer() const
static constexpr auto TID
unsigned short size() const
Get the number of SeedingLayerSets.
T const * product() const
const std::vector< uint8_t > & amplitudes() const
Global3DVector GlobalVector