18 theLsb(conf.getParameter<edm::
ParameterSet>(
"TripletsPSet")),
19 writeTriplets_(conf.getParameter<bool>(
"writeTriplets")),
20 seedOnMiddle_(conf.existsAs<bool>(
"seedOnMiddle") ? conf.getParameter<bool>(
"seedOnMiddle") :
false),
21 rescaleError_(conf.existsAs<double>(
"rescaleError") ? conf.getParameter<double>(
"rescaleError") : 1.0),
22 tripletsVerbosity_(conf.getParameter<edm::
ParameterSet>(
"TripletsPSet").getUntrackedParameter<uint32_t>(
"debugLevel",0)),
23 seedVerbosity_(conf.getUntrackedParameter<uint32_t>(
"seedDebugLevel",0)),
24 helixVerbosity_(conf.getUntrackedParameter<uint32_t>(
"helixDebugLevel",0)),
25 check_(conf.getParameter<edm::
ParameterSet>(
"ClusterCheckPSet")),
26 maxTriplets_(conf.getParameter<int32_t>(
"maxTriplets")),
27 maxSeeds_(conf.getParameter<int32_t>(
"maxSeeds"))
31 float originradius = regionConf.
getParameter<
double>(
"originRadius");
32 float halflength = regionConf.
getParameter<
double>(
"originHalfLength");
33 float originz = regionConf.
getParameter<
double>(
"originZPosition");
44 produces<TrajectorySeedCollection>();
45 if (
writeTriplets_) produces<edm::OwnVector<TrackingRecHit> >(
"cosmicTriplets");
92 edm::LogError(
"TooManyClusters") <<
"Found too many clusters (" << clustsOrZero <<
"), bailing out.\n";
98 bool seedsOk =
seeds(*output,es);
106 outtriplets->push_back(hit1->
clone());
107 outtriplets->push_back(hit2->
clone());
108 outtriplets->push_back(hit3->
clone());
117 ev.
put(outtriplets,
"cosmicTriplets");
150 float iy1 = ihit1->globalPosition().y();
151 float oy1 = ohit1->globalPosition().y();
152 float iy2 = ihit2->globalPosition().y();
153 float oy2 = ohit2->globalPosition().y();
157 return (iy1 != iy2 ? (iy1 > iy2) : (oy1 > oy2));
159 }
else if (oy2 - iy2 > 0) {
163 return (iy1 != iy2 ? (iy1 < iy2) : (oy1 < oy2));
169 using namespace ctfseeding;
174 SeedingLayerSets::const_iterator iLss;
178 for (iLss = lss.begin(); iLss != lss.end(); iLss++){
181 throw cms::Exception(
"CtfSpecialSeedGenerator") <<
"You are using " << ls.size() <<
" layers in set instead of 3 ";
185 std::vector<SeedingHit> innerHits =
region_.
hits(e, es, &ls[0]);
186 std::vector<SeedingHit> middleHits =
region_.
hits(e, es, &ls[1]);
187 std::vector<SeedingHit> outerHits =
region_.
hits(e, es, &ls[2]);
188 std::vector<SeedingHit>::const_iterator iOuterHit,iMiddleHit,iInnerHit;
192 <<
" (" << (iLss - lss.begin()) <<
"): # = "
193 << innerHits.size() <<
"/" << middleHits.size() <<
"/" << outerHits.size() << std::endl;
198 std::vector<TTRH> innerTTRHs, middleTTRHs, outerTTRHs;
201 std::vector<bool> innerOk( innerHits.size(),
true);
202 std::vector<bool> middleOk(middleHits.size(),
true);
203 std::vector<bool> outerOk( outerHits.size(),
true);
208 for (iOuterHit = outerHits.begin(), idx = 0; iOuterHit != outerHits.end(); ++idx, ++iOuterHit){
209 outerTTRHs.push_back(ls[2].hitBuilder()->
build((**iOuterHit).hit()));
212 for (iMiddleHit = middleHits.begin(), idx = 0; iMiddleHit != middleHits.end(); ++idx, ++iMiddleHit){
213 middleTTRHs.push_back(ls[1].hitBuilder()->
build((**iMiddleHit).hit()));
216 for (iInnerHit = innerHits.begin(), idx = 0; iInnerHit != innerHits.end(); ++idx, ++iInnerHit){
217 innerTTRHs.push_back(ls[0].hitBuilder()->
build((**iInnerHit).hit()));
226 for (iOuterHit = outerHits.begin(); iOuterHit != outerHits.end(); iOuterHit++){
227 idx = iOuterHit - outerHits.begin();
228 TTRH & outerTTRH = outerTTRHs[idx];
229 GlobalPoint outerpos = outerTTRH->globalPosition();
230 bool outerok = outerOk[idx];
237 for (iMiddleHit = middleHits.begin(); iMiddleHit != middleHits.end(); iMiddleHit++){
238 idx = iMiddleHit - middleHits.begin();
239 TTRH & middleTTRH = middleTTRHs[idx];
240 GlobalPoint middlepos = middleTTRH->globalPosition();
241 bool middleok = middleOk[idx];
248 for (iInnerHit = innerHits.begin(); iInnerHit != innerHits.end(); iInnerHit++){
249 idx = iInnerHit - innerHits.begin();
250 TTRH & innerTTRH = innerTTRHs[idx];
251 GlobalPoint innerpos = innerTTRH->globalPosition();
252 bool innerok = innerOk[idx];
261 || outerpos.
y() < innerpos.
y()
264 || outerpos.
y() > innerpos.
y()
269 if (
goodTriplet(innerpos,middlepos,outerpos,minRho)) {
275 edm::LogError(
"TooManyTriplets") <<
"Found too many triplets, bailing out.\n";
280 << innerpos <<
" + " << middlepos <<
" + " << outerpos << std::endl;
284 << innerpos <<
" + " << middlepos <<
" + " << outerpos << std::endl;
295 <<
" (" << (iLss - lss.begin()) <<
"): # = "
296 << innerHits.size() <<
"/" << middleHits.size() <<
"/" << outerHits.size()
308 int subdet =
detid.subdetId();
312 if (
typeid(*hit) ==
typeid(SiStripMatchedRecHit2D)) {
313 const SiStripMatchedRecHit2D *mhit =
static_cast<const SiStripMatchedRecHit2D *
>(hit);
319 }
else if (
typeid(*hit) ==
typeid(SiStripRecHit2D)) {
320 return checkCharge(static_cast<const SiStripRecHit2D &>(*hit), subdet);
329 const SiStripCluster *clust = (hit.cluster().isNonnull() ? hit.cluster().get() : hit.cluster_regional().get());
333 <<
", detid = " << hit.geographicalId().rawId() <<
", firstStrip = " << clust->
firstStrip() << std::endl;
336 <<
", detid = " << hit.geographicalId().rawId() <<
", firstStrip = " << clust->
firstStrip() << std::endl;
343 std::vector<TTRH>::const_iterator it = hits.begin(),
start = it,
end = hits.end();
344 std::vector<bool>::iterator
ok = oks.begin(), okStart =
ok;
345 while (
start < end) {
346 DetId lastid = (*start)->geographicalId();
347 for (it =
start + 1; (it <
end) && ((*it)->geographicalId() == lastid); ++it) {
352 std::cerr <<
"SimpleCosmicBONSeeder: Marking noisy module " << lastid.
rawId() <<
", it has " << (it-
start) <<
" rechits"
358 std::cerr <<
"SimpleCosmicBONSeeder: Not marking noisy module " << lastid.
rawId() <<
", it has " << (it-
start) <<
" rechits"
367 float dyOM = outer.
y() - middle.
y(), dyIM = inner.
y() - middle.
y();
368 if ((dyOM * dyIM > 0) && (fabs(dyOM)>10) && (fabs(dyIM)>10)) {
372 float dzOM = outer.
z() - middle.
z(), dzIM = inner.
z() - middle.
z();
373 if ((dzOM * dzIM > 0) && (fabs(dzOM)>50) && (fabs(dzIM)>50)) {
379 if (theCircle.
rho() < minRho) {
387 std::pair<GlobalVector,int>
391 std::cout <<
"DEBUG PZ =====" << std::endl;
392 FastHelix helix(inner,middle,outer,iSetup);
394 std::cout <<
"FastHelix P = " << gv <<
"\n";
401 double rho = theCircle.
rho();
405 double pt = 0.01 * rho * (0.3*tesla.
z());
408 double dx1 = outer.
x()-theCircle.
x0();
409 double dy1 = outer.
y()-theCircle.
y0();
410 double py = pt*dx1/
rho, px = -pt*dy1/
rho;
411 if(px*(middle.
x() - outer.
x()) + py*(middle.
y() - outer.
y()) < 0.) {
412 px *= -1.; py *= -1.;
416 double dz = inner.
z() - outer.
z();
417 double sinphi = ( dx1*(inner.
y()-theCircle.
y0()) - dy1*(inner.
x()-theCircle.
x0())) / (rho *
rho);
418 double dphi =
std::abs(std::asin(sinphi));
419 double pz = pt * dz / (dphi *
rho);
421 int myq = ((theCircle.
x0()*py - theCircle.
y0()*px) / tesla.
z()) > 0. ? +1 : -1;
423 std::pair<GlobalVector,int> mypq(
GlobalVector(px,py,pz),myq);
426 std::cout <<
"Gio: pt = " << pt << std::endl;
427 std::cout <<
"Gio: dz = " << dz <<
", sinphi = " << sinphi <<
", dphi = " << dphi <<
", dz/drphi = " << (dz/dphi/
rho) << std::endl;
430 std::cout <<
"Gio's fit P = " << mypq.first <<
"\n";
431 std::cout <<
"Gio's fit Q = " << myq <<
"\n";
445 toGlobal((*(trip.
inner())).localPosition());
448 toGlobal((*(trip.
middle())).localPosition());
451 toGlobal((*(trip.
outer())).localPosition());
454 std::cout <<
"Processing triplet " << it <<
": " << inner <<
" + " << middle <<
" + " << outer << std::endl;
456 if ( (outer.y()-inner.
y())*outer.y() < 0 ) {
458 std::swap(const_cast<TransientTrackingRecHit::ConstRecHitPointer &>(trip.
inner()),
459 const_cast<TransientTrackingRecHit::ConstRecHitPointer &>(trip.
outer()) );
464 std::cout <<
"The seed was going away from CMS! swapped in <-> out" << std::endl;
465 std::cout <<
"Processing swapped triplet " << it <<
": " << inner <<
" + " << middle <<
" + " << outer << std::endl;
470 std::pair<GlobalVector,int> pq =
pqFromHelixFit(inner,middle,outer,iSetup);
472 float ch = pq.second;
473 float Mom =
sqrt( gv.
x()*gv.
x() + gv.
y()*gv.
y() + gv.
z()*gv.
z() );
475 if(Mom > 10000 ||
isnan(Mom)) {
477 std::cout <<
"Processing triplet " << it <<
": fail for momentum." << std::endl;
483 std::cout <<
"Processing triplet " << it <<
": fail for pt = " << gv.
perp() <<
" < ptMin = " <<
region_.
ptMin() << std::endl;
488 if((outer.y()-inner.
y())>0){
490 std::cout <<
"Processing triplet " << it <<
": downgoing." << std::endl;
493 gv = -1*gv; ch = -1.*ch;
496 std::cout <<
"Processing triplet " << it <<
": upgoing." << std::endl;
500 if (( gv.
z() * (outer.z()-inner.
z()) > 0 ) && ( fabs(outer.z()-inner.
z()) > 5) && (fabs(gv.
z()) > .01)) {
501 std::cout <<
"ORRORE: outer.z()-inner.z() = " << (outer.z()-inner.
z()) <<
", gv.z() = " << gv.
z() << std::endl;
512 std::cout <<
"Processing triplet " << it <<
". start from " << std::endl;
513 std::cout <<
" X = " << outer <<
", P = " << gv << std::endl;
519 TSOS propagated, updated;
521 for (
size_t ih = 0; ih < 3; ++ih) {
524 std::cout <<
"Stopping at middle hit, as requested." << std::endl;
528 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
"." << std::endl;
530 propagated = propagator->
propagate(CosmicSeed,
tracker->idToDet((*seedHits[ih]).geographicalId())->surface());
532 propagated = propagator->
propagate(updated,
tracker->idToDet((*seedHits[ih]).geographicalId())->surface());
534 if (!propagated.isValid()) {
536 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
": failed propagation." << std::endl;
540 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
": propagated state = " << propagated;
546 if (!updated.isValid()) {
548 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
": failed update." << std::endl;
552 std::cout <<
"Processing triplet " << it <<
", hit " << ih <<
": updated state = " << updated;
555 if (!fail && updated.isValid() && (updated.globalMomentum().perp() <
region_.
ptMin())) {
557 std::cout <<
"Processing triplet " << it <<
558 ": failed for final pt " << updated.globalMomentum().perp() <<
" < " <<
region_.
ptMin() << std::endl;
561 if (!fail && updated.isValid() && (updated.globalMomentum().mag() <
pMin_)) {
563 std::cout <<
"Processing triplet " << it <<
564 ": failed for final p " << updated.globalMomentum().perp() <<
" < " <<
pMin_ << std::endl;
570 std::cout <<
"Processing triplet " << it <<
", rescale error by " <<
rescaleError_ <<
": state BEFORE rescaling " << updated;
571 std::cout <<
" Cartesian error (X,P) before rescaling= \n" << updated.cartesianError().matrix() << std::endl;
576 std::cout <<
"Processed triplet " << it <<
": success (saved as #"<<output.size()<<
") : "
577 << inner <<
" + " << middle <<
" + " << outer << std::endl;
578 std::cout <<
" pt = " << updated.globalMomentum().perp() <<
579 " eta = " << updated.globalMomentum().eta() <<
580 " phi = " << updated.globalMomentum().phi() <<
581 " ch = " << updated.charge() << std::endl;
585 std::cout <<
" X = " << updated.globalPosition() <<
", P = " << updated.globalMomentum() << std::endl;
587 std::cout <<
" Cartesian error (X,P) = \n" << updated.cartesianError().matrix() << std::endl;
596 edm::LogError(
"TooManySeeds") <<
"Found too many seeds, bailing out.\n";
T getParameter(std::string const &) const
GlobalTrackingRegion region_
FTS stateAtVertex() const
const OuterRecHit & outer() const
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_
const GlobalTrajectoryParameters & parameters() const
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
OrderedHitTriplets hitTriplets
Global3DPoint GlobalPoint
const InnerRecHit & inner() const
virtual void produce(edm::Event &e, const edm::EventSetup &c)
const MiddleRecHit & middle() const
uint16_t firstStrip() const
void checkNoisyModules(const std::vector< TransientTrackingRecHit::RecHitPointer > &hits, std::vector< bool > &oks) const
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
bool checkCharge(const TrackingRecHit *hit) const
SimpleCosmicBONSeeder(const edm::ParameterSet &conf)
uint32_t rawId() const
get the raw id
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
virtual TrackingRegion::Hits hits(const edm::Event &ev, const edm::EventSetup &es, const ctfseeding::SeedingLayer *layer) const
get hits from layer compatible with region constraints
SeedingLayerSetsBuilder theLsb
std::vector< TrajectorySeed > TrajectorySeedCollection
const T & max(const T &a, const T &b)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
virtual float ptMin() const
minimal pt of interest
PropagatorWithMaterial * thePropagatorOp
GlobalVector momentum() const
edm::ESHandle< MagneticField > magfield
size_t tooManyClusters(const edm::Event &e) const
edm::ESHandle< TransientTrackingRecHitBuilder > TTTRHBuilder
uint32_t tripletsVerbosity_
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
PropagatorWithMaterial * thePropagatorAl
const AlgebraicSymMatrix66 & matrix() const
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Surface &) const
virtual TrackingRecHit * clone() const =0
TrajectoryStateOnSurface TSOS
void rescaleError(double factor)
bool seeds(TrajectorySeedCollection &output, const edm::EventSetup &iSetup)
bool goodTriplet(const GlobalPoint &inner, const GlobalPoint &middle, const GlobalPoint &outer, const double &minRho) const
bool operator()(const OrderedHitTriplet &trip1, const OrderedHitTriplet &trip2) const
TrajectoryStateOnSurface update(const TrajectoryStateOnSurface &, const TransientTrackingRecHit &) const
ctfseeding::SeedingLayerSets layers(const edm::EventSetup &es) const
std::vector< std::string > layerTripletNames_
TransientTrackingRecHit::ConstRecHitPointer SeedingHit
std::vector< int32_t > chargeThresholds_
std::vector< int32_t > maxHitsPerModule_
DetId geographicalId() const
TrackCharge charge() const
virtual unsigned int size() const
const std::vector< uint8_t > & amplitudes() const
std::vector< SeedingLayer > SeedingLayers
Global3DVector GlobalVector
std::vector< std::vector< SeedingLayer > > SeedingLayerSets