70 theDTRecHitLabel(par.getParameter<
InputTag>(
"DTRecSegmentLabel")),
71 theCSCRecHitLabel(par.getParameter<
InputTag>(
"CSCRecSegmentLabel")),
72 theCSCSegmentsLabel(par.getParameter<
InputTag>(
"CSCSegmentLabel")),
73 theDT4DRecSegmentLabel(par.getParameter<
InputTag>(
"DT4DRecSegmentLabel"))
91 category_ =
"MuonShowerInformationFiller";
93 for (
int istat = 0; istat < 4; istat++) {
147 for (
int istat = 0; istat < 4; istat++) {
194 DTRecSegment4DCollection::id_iterator chamberIdIt;
195 for (chamberIdIt = dtSegments->id_begin();
196 chamberIdIt != dtSegments->id_end();
199 if (*chamberIdIt != chamberId)
continue;
205 iseg!=range.second;++iseg) {
206 if (iseg->dimension() != 4)
continue;
215 for ( CSCSegmentCollection::id_iterator chamberId = cscSegments->id_begin();
216 chamberId != cscSegments->id_end(); ++chamberId) {
218 if ((*chamberId).chamber() != did.chamber())
continue;
224 iseg!=range.second;++iseg) {
225 if (iseg->dimension() != 3)
continue;
236 if (segments.empty())
return allhitscorrelated;
241 if (segments.size() == 1)
return allhitscorrelated;
243 for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseg = segments.begin() + 1;
244 iseg != segments.end(); ++iseg) {
249 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit1 = hits1.begin();
250 ihit1 != hits1.end(); ++ihit1 ) {
252 bool usedbefore =
false;
255 GlobalPoint gp1dinsegHit = (*ihit1)->globalPosition();
257 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit2 = allhitscorrelated.begin();
258 ihit2 != allhitscorrelated.end(); ++ihit2 ) {
262 GlobalPoint gp1dinsegHit2 = (*ihit2)->globalPosition();
264 if ( (gp1dinsegHit2 - gp1dinsegHit).
mag() < 1.0 ) usedbefore =
true;
267 if ( !usedbefore ) allhitscorrelated.push_back(*ihit1);
271 return allhitscorrelated;
283 if ( muonRecHits.empty() )
return muonRecHits;
289 stable_sort(muonRecHits.begin(), muonRecHits.end(),
AbsLessDPhi(refpoint));
291 for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihit = muonRecHits.begin(); ihit != muonRecHits.end() - 1; ++ihit) {
292 if (fabs(
deltaPhi((*(ihit+1))->globalPosition().
phi(), (*ihit)->globalPosition().phi() )) <
step) {
293 result.push_back(*ihit);
299 LogTrace(
category_) <<
"phi front: " << muonRecHits.front()->globalPosition().phi() << endl;
300 LogTrace(
category_) <<
"phi back: " << muonRecHits.back()->globalPosition().phi() << endl;
313 if ( muonRecHits.empty() )
return muonRecHits;
319 stable_sort(muonRecHits.begin(), muonRecHits.end(),
AbsLessDTheta(refpoint));
321 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = muonRecHits.begin(); ihit != muonRecHits.end() - 1; ++ihit) {
322 if (fabs((*(ihit+1))->globalPosition().
theta() - (*ihit)->globalPosition().theta() ) < step) {
323 result.push_back(*ihit);
339 if ( muonRecHits.empty() )
return muonRecHits;
341 stable_sort(muonRecHits.begin(), muonRecHits.end(),
LessPerp());
343 MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator
344 seedhit = min_element(muonRecHits.begin(), muonRecHits.end(),
LessPerp());
346 MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihigh = seedhit;
347 MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ilow = seedhit;
350 while (ihigh != muonRecHits.end()-1 && ( fabs((*(ihigh+1))->globalPosition().
perp() - (*ihigh)->globalPosition().perp() ) < step) ) {
353 while (ilow != muonRecHits.begin() && ( fabs((*ilow)->globalPosition().perp() - (*(ilow -1))->globalPosition().perp()) < step ) ) {
368 vector<const GeomDet*>
total;
371 LogTrace(
category_) <<
"Consider a track " << track.
p() <<
" eta: " << track.
eta() <<
" phi " << track.
phi() << endl;
380 vector<GlobalPoint> allCrossingPoints;
382 const vector<DetLayer*>& dtlayers =
theService->detLayerGeometry()->allDTLayers();
384 for (vector<DetLayer*>::const_iterator iLayer = dtlayers.begin(); iLayer != dtlayers.end(); ++iLayer) {
390 if ((fabs(xPoint.
y()) < 1000.0) && (fabs(xPoint.
z()) < 1500 ) &&
391 (!(xPoint.
y() == 0 && xPoint.
x() == 0 && xPoint.
z() == 0))) allCrossingPoints.push_back(xPoint);
394 stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(),
LessMag(innerPos) );
396 vector<const GeomDet*> tempDT;
398 for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
401 vector<const GeomDet*>::const_iterator
begin = tempDT.begin();
402 vector<const GeomDet*>::const_iterator
end = tempDT.end();
405 total.push_back(*begin);
409 allCrossingPoints.clear();
411 const vector<DetLayer*>& csclayers =
theService->detLayerGeometry()->allCSCLayers();
412 for (vector<DetLayer*>::const_iterator iLayer = csclayers.begin(); iLayer != csclayers.end(); ++iLayer) {
417 if ((fabs(xPoint.
y()) < 1000.0) && (fabs(xPoint.
z()) < 1500.0)
418 && (!(xPoint.
y() == 0 && xPoint.
x() == 0 && xPoint.
z() == 0))) allCrossingPoints.push_back(xPoint);
420 stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(),
LessMag(innerPos) );
422 vector<const GeomDet*> tempCSC;
423 for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
426 vector<const GeomDet*>::const_iterator
begin = tempCSC.begin();
427 vector<const GeomDet*>::const_iterator
end = tempCSC.end();
430 total.push_back(*begin);
456 float radius = cyl.radius();
460 float a = p1.
x() - slope * p1.
y();
462 float n2 = (1 + slope *
slope);
463 float n1 = 2*a*
slope;
468 if ( n1*n1 - 4*n2*n0 > 0 ) {
469 y1 = (-n1 +
sqrt(n1*n1 - 4*n2*n0) ) / (2 * n2);
470 y2 = (-n1 -
sqrt(n1*n1 - 4*n2*n0) ) / (2 * n2);
473 float x1 = p1.
x() + slope * (y1 - p1.
y());
474 float x2 = p1.
x() + slope * (y2 - p1.
y());
476 float slopeZ = dp.
z()/dp.
y();
478 float z1 = p1.
z() + slopeZ * (y1 - p1.
y());
479 float z2 = p1.
z() + slopeZ * (y2 - p1.
y());
482 if ((p2.x()*x1 > 0) && (y1*p2.y() > 0) && (z1*p2.z() > 0)) {
507 float diskZ = disk.position().z();
508 int endcap = diskZ > 0 ? 1 : (diskZ < 0 ? -1 : 0);
509 diskZ = diskZ + endcap*
dynamic_cast<const SimpleDiskBounds&
>(disk.bounds()).thickness()/2.;
513 float slopeZ = dp.
z()/dp.
y();
514 float y1 = diskZ / slopeZ;
516 float slopeX = dp.
z()/dp.
x();
517 float x1 = diskZ / slopeX;
537 if ( gp.
z() < -680.0 ) { minwheel = -3; maxwheel = -3;}
538 else if ( gp.
z() < -396.0 ) { minwheel = -2; maxwheel = -1;}
539 else if (gp.
z() < -126.8) { minwheel = -2; maxwheel = 0; }
540 else if (gp.
z() < 126.8) { minwheel = -1; maxwheel = 1; }
541 else if (gp.
z() < 396.0) { minwheel = 0; maxwheel = 2; }
542 else if (gp.
z() < 680.0) { minwheel = 1; maxwheel = 2; }
543 else { minwheel = 3; maxwheel = 3; }
546 if ( gp.
perp() > 680.0 && gp.
perp() < 755.0 ) station = 4;
547 else if ( gp.
perp() > 580.0 ) station = 3;
548 else if ( gp.
perp() > 480.0 ) station = 2;
549 else if ( gp.
perp() > 380.0 ) station = 1;
554 float phistep =
M_PI/6;
556 float phigp = (float)gp.
phi();
558 if ( fabs(
deltaPhi(phigp, 0*phistep)) < phistep ) sectors.push_back(1);
559 if ( fabs(
deltaPhi(phigp, phistep)) < phistep ) sectors.push_back(2);
560 if ( fabs(
deltaPhi(phigp, 2*phistep)) < phistep ) sectors.push_back(3);
561 if ( fabs(
deltaPhi(phigp, 3*phistep)) < phistep ) {
562 sectors.push_back(4);
563 if (station == 4) sectors.push_back(13);
565 if ( fabs(
deltaPhi(phigp, 4*phistep)) < phistep ) sectors.push_back(5);
566 if ( fabs(
deltaPhi(phigp, 5*phistep)) < phistep ) sectors.push_back(6);
567 if ( fabs(
deltaPhi(phigp, 6*phistep)) < phistep ) sectors.push_back(7);
568 if ( fabs(
deltaPhi(phigp, 7*phistep)) < phistep ) sectors.push_back(8);
569 if ( fabs(
deltaPhi(phigp, 8*phistep)) < phistep ) sectors.push_back(9);
570 if ( fabs(
deltaPhi(phigp, 9*phistep)) < phistep ) {
571 sectors.push_back(10);
572 if (station == 4) sectors.push_back(14);
574 if ( fabs(
deltaPhi(phigp, 10*phistep)) < phistep ) sectors.push_back(11);
575 if ( fabs(
deltaPhi(phigp, 11*phistep)) < phistep ) sectors.push_back(12);
578 LogTrace(
category_) <<
"number of sectors to consider: " << sectors.size() << endl;
579 LogTrace(
category_) <<
"station: " << station <<
" wheels: " << minwheel <<
" " << maxwheel << endl;
581 vector<const GeomDet*>
result;
582 if (station > 4 || station < 1)
return result;
583 if (minwheel > 2 || maxwheel < -2)
return result;
585 for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector ) {
586 for (
int iwheel = minwheel; iwheel != maxwheel + 1; ++iwheel) {
587 DTChamberId chamberid(iwheel, station, (*isector));
588 result.push_back(
theService->trackingGeometry()->idToDet(chamberid));
592 LogTrace(
category_) <<
"number of GeomDets for this track: " << result.size() << endl;
606 if (gp.
z() > 0) {endcap = 1;}
else {endcap = 2;}
612 if ( fabs(gp.
z()) > 1000. && fabs(gp.
z()) < 1055.0 ) {
615 else if ( fabs(gp.
z()) > 910.0 && fabs(gp.
z()) < 965.0) {
618 else if ( fabs(gp.
z()) > 800.0 && fabs(gp.
z()) < 860.0) {
621 else if ( fabs(gp.
z()) > 570.0 && fabs(gp.
z()) < 730.0) {
627 float phistep1 =
M_PI/18.;
628 float phistep2 =
M_PI/9.;
629 float phigp = (float)gp.
phi();
637 if (gp.
perp() > 100 && gp.
perp() < 270) ring = 1;
638 else if (gp.
perp() > 270 && gp.
perp() < 450) ring = 2;
639 else if (gp.
perp() > 450 && gp.
perp() < 695) ring = 3;
640 else if (gp.
perp() > 100 && gp.
perp() < 270) ring = 4;
643 else if (station == 2) {
645 if (gp.
perp() > 140 && gp.
perp() < 350) ring = 1;
646 else if (gp.
perp() > 350 && gp.
perp() < 700) ring = 2;
649 else if (station == 3) {
651 if (gp.
perp() > 160 && gp.
perp() < 350) ring = 1;
652 else if (gp.
perp() > 350 && gp.
perp() < 700) ring = 2;
655 else if (station == 4) {
657 if (gp.
perp() > 175 && gp.
perp() < 350) ring = 1;
658 else if (gp.
perp() > 350 && gp.
perp() < 700) ring = 2;
662 if (station > 1 && ring == 1) {
665 for (
int i = 0;
i < 18;
i++) {
666 if ( fabs(
deltaPhi(phigp,
i*phistep2)) < phistep2 ) sectors.push_back(
i+1);
672 for (
int i = 0;
i < 36;
i++) {
673 if ( fabs(
deltaPhi(phigp,
i*phistep1)) < phistep1 ) sectors.push_back(
i+1);
681 LogTrace(
category_) <<
"CSC number of sectors to consider: " << sectors.size() << endl;
685 vector<const GeomDet*>
result;
686 if (station > 4 || station < 1)
return result;
687 if (endcap == 0)
return result;
688 if (ring == -1)
return result;
693 for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
694 for (
int ilayer = minlayer; ilayer != maxlayer + 1; ++ ilayer) {
695 CSCDetId cscid(endcap, station, ring, (*isector), ilayer);
696 result.push_back(
theService->trackingGeometry()->idToDet(cscid));
715 vector<MuonRecHitContainer> muonRecHits(4);
718 vector<TransientTrackingRecHit::ConstRecHitContainer> muonCorrelatedHits(4);
725 bool dtOverlapToCheck =
false;
726 bool cscOverlapToCheck =
false;
728 for (vector<const GeomDet*>::const_iterator igd = compatibleLayers.begin(); igd != compatibleLayers.end(); igd++ ) {
731 DetId geoId = (*igd)->geographicalId();
742 int wheel =
detid.wheel();
746 TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
747 TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
749 for (; hits_begin!= hits_end;++hits_begin) {
750 muonCorrelatedHits.at(station-1).push_back(*hits_begin);
754 if (
abs(wheel) == 2 && station != 4 && station != 1) dtOverlapToCheck =
true;
763 vector<const TrackingRecHit*> subrechits = (*rechit).recHits();
764 for (vector<const TrackingRecHit*>::iterator irechit = subrechits.begin(); irechit != subrechits.end(); ++irechit) {
776 int ring = did.ring();
780 TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
781 TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
783 for (; hits_begin!= hits_end;++hits_begin) {
784 muonCorrelatedHits.at(station-1).push_back(*hits_begin);
787 if ((station == 1 && ring == 3) && dtOverlapToCheck) cscOverlapToCheck =
true;
793 if (!cscOverlapToCheck) {
800 if (temp.empty())
continue;
803 if (temp.size() > 1) {
804 center = (temp.front()->globalPosition().perp() + temp.back()->globalPosition().perp())/2.;
806 center = temp.front()->globalPosition().perp();
811 muonRecHits.at(2).insert(muonRecHits.at(2).end(),tmpCSC1.begin(),tmpCSC1.end());
813 muonRecHits.at(1).insert(muonRecHits.at(1).end(),tmpCSC1.begin(),tmpCSC1.end());
824 for (
int stat = 0; stat < 4; stat++) {
845 for (
int stat = 0; stat != 4; stat++ ) {
846 if (!muonRecHits[stat].
empty()) {
847 stable_sort(muonRecHits[stat].
begin(), muonRecHits[stat].
end(),
LessPhi());
850 for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator
iseed = muonRecHits[stat].
begin();
iseed != muonRecHits[stat].end(); ++
iseed) {
851 if (!(*iseed)->isValid())
continue;
853 muonRecHitsPhiTemp.clear();
855 if (muonRecHitsPhiTemp.size() > 1) {
856 float dphi = fabs(
deltaPhi((
float)muonRecHitsPhiTemp.back()->globalPosition().phi(), (float)muonRecHitsPhiTemp.front()->globalPosition().phi()));
857 if (dphi > dphimax) {
859 muonRecHitsPhiBest = muonRecHitsPhiTemp;
865 if (!muonRecHitsPhiBest.empty()) {
866 muonRecHits[stat] = muonRecHitsPhiBest;
868 muonRecHits[stat].front();
869 GlobalPoint refpoint = muonRecHits[stat].front()->globalPosition();
874 if (!muonCorrelatedHits.at(stat).empty()) {
877 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator
iseed = muonCorrelatedHits.at(stat).begin();
iseed != muonCorrelatedHits.at(stat).end(); ++
iseed) {
878 if (!(*iseed)->isValid())
continue;
880 muonRecHitsThetaTemp.clear();
881 muonRecHitsThetaTemp =
findThetaCluster(muonCorrelatedHits.at(stat), refpoint);
883 if (muonRecHitsThetaTemp.size() > 1) {
884 float dtheta = fabs((
float)muonRecHitsThetaTemp.back()->globalPosition().theta() - (float)muonRecHitsThetaTemp.front()->globalPosition().theta());
885 if (dtheta > dthetamax) {
887 muonRecHitsThetaBest = muonRecHitsThetaTemp;
893 if (muonRecHitsThetaBest.size() > 1 && muonRecHitsPhiBest.size() > 1)
894 theStationShowerDeltaR.at(stat) =
sqrt(
pow(muonRecHitsPhiBest.front()->globalPosition().phi()-muonRecHitsPhiBest.back()->globalPosition().phi(),2)+
pow(muonRecHitsThetaBest.front()->globalPosition().theta()-muonRecHitsThetaBest.back()->globalPosition().theta(),2));
899 LogTrace(
category_) <<
"deltaR around a track containing all the station hits, by station "
double p() const
momentum vector magnitude
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
T getParameter(std::string const &) const
std::pair< const_iterator, const_iterator > range
iterator range
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
static const double slope[3]
Geom::Phi< T > phi() const
Global3DPoint GlobalPoint
Geom::Theta< T > theta() const
bool isGlobalMuon() const
GlobalPoint globalPosition() const
double phi() const
azimuthal angle of momentum vector
std::vector< int > nStationHits
number of all the muon RecHits per chamber crossed by a track (1D hits)
bool isStandAloneMuon() const
uint32_t rawId() const
get the raw id
static TransientTrackingRecHit::ConstRecHitContainer breakInSubRecHits(TransientTrackingRecHit::ConstRecHitPointer, int granularity)
takes a muon rechit and returns its sub-rechits given a certain granularity
std::vector< int > nStationCorrelatedHits
number of the muon RecHits used by segments per chamber crossed by a track
double eta() const
pseudorapidity of momentum vector
static const int maxLayerId
highest layer id
Abs< T >::type abs(const T &t)
static const int minLayerId
lowest layer id. 0 indicates a full SL
DetId geographicalId() const
The label of this GeomDet.
static const int maxSuperLayerId
highest superlayer id
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) ...
virtual TrackRef outerTrack() const
reference to Track reconstructed in the muon detector only
std::vector< ConstRecHitPointer > ConstRecHitContainer
static const int minSuperLayerId
loweset super layer id. 0 indicates a full chamber
virtual const BoundDisk & specificSurface() const GCC11_FINAL
T perp() const
Magnitude of transverse component.
std::vector< float > stationShowerSizeT
the transverse size of the hit cluster
Detector det() const
get the detector field from this detid
virtual const BoundCylinder & specificSurface() const GCC11_FINAL
Extension of the interface.
static MuonRecHitPointer specificBuild(const GeomDet *geom, const TrackingRecHit *rh)
std::vector< MuonRecHitPointer > MuonRecHitContainer
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
Power< A, B >::type pow(const A &a, const B &b)
virtual TrackRef globalTrack() const
reference to Track reconstructed in both tracked and muon detector
std::vector< float > stationShowerDeltaR
the radius of the cone containing the all the hits around the track