363 edm::LogInfo(
"RPCRecHitValid") <<
"Cannot find simHit collection\n";
370 edm::LogInfo(
"RPCRecHitValid") <<
"Cannot find recHit collection\n";
377 edm::LogInfo(
"RPCRecHitValid") <<
"Cannot find TrackingParticle collection\n";
384 edm::LogInfo(
"RPCRecHitValid") <<
"Cannot find TrackingParticle to SimHit association map\n";
391 edm::LogInfo(
"RPCRecHitValid") <<
"Cannot find muon collection\n";
395 typedef edm::PSimHitContainer::const_iterator SimHitIter;
397 typedef std::vector<TrackPSimHitRef> SimHitRefs;
400 SimHitRefs muonSimHits, pthrSimHits;
402 for (
int i = 0,
n = simParticleHandle->size();
i <
n; ++
i) {
404 if (simParticle->pt() < 1.0
or simParticle->p() < 2.5)
408 SimHitRefs simHitsFromParticle;
409 auto range = std::equal_range(simHitsTPAssoc->begin(),
410 simHitsTPAssoc->end(),
413 for (
auto simParticleToHit =
range.first; simParticleToHit !=
range.second; ++simParticleToHit) {
414 auto simHit = simParticleToHit->second;
419 simHitsFromParticle.push_back(simParticleToHit->second);
421 const int nRPCHit = simHitsFromParticle.size();
422 const bool hasRPCHit = nRPCHit > 0;
424 if (
abs(simParticle->pdgId()) == 13) {
425 muonSimHits.insert(muonSimHits.end(), simHitsFromParticle.begin(), simHitsFromParticle.end());
428 int nRPCHitBarrel = 0;
429 int nRPCHitEndcap = 0;
430 for (
const auto &
simHit : simHitsFromParticle) {
431 const RPCDetId rpcDetId = static_cast<const RPCDetId>(
simHit->detUnitId());
432 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rpcDetId));
436 if (rpcDetId.
region() == 0)
444 if (nRPCHitBarrel and nRPCHitEndcap) {
449 }
else if (nRPCHitBarrel) {
454 }
else if (nRPCHitEndcap) {
465 pthrSimHits.insert(pthrSimHits.end(), simHitsFromParticle.begin(), simHitsFromParticle.end());
469 switch (simParticle->pdgId()) {
508 int nRefHitBarrel = 0, nRefHitEndcap = 0;
509 for (
const auto &
simHit : muonSimHits) {
510 const RPCDetId detId = static_cast<const RPCDetId>(
simHit->detUnitId());
511 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
513 const int region = roll->id().region();
514 const int ring = roll->id().ring();
516 const int station = roll->id().station();
538 for (
const auto &
simHit : pthrSimHits) {
539 const RPCDetId detId = static_cast<const RPCDetId>(
simHit->detUnitId());
540 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId()));
542 const int region = roll->id().region();
543 const int ring = roll->id().ring();
545 const int station = roll->id().station();
568 int sumClusterSizeBarrel = 0, sumClusterSizeEndcap = 0;
569 int nRecHitBarrel = 0, nRecHitEndcap = 0;
570 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
571 const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
572 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId()));
576 const int region = roll->id().region();
577 const int ring = roll->id().ring();
579 const int station = roll->id().station();
583 const double time = recHitIter->timeError() >= 0 ? recHitIter->time() : recHitIter->BunchX() * 25;
589 sumClusterSizeBarrel += recHitIter->clusterSize();
598 sumClusterSizeEndcap += recHitIter->clusterSize();
606 if (roll->isIRPC()) {
612 const double nRecHit = nRecHitBarrel + nRecHitEndcap;
616 const int sumClusterSize = sumClusterSizeBarrel + sumClusterSizeEndcap;
619 if (nRecHitBarrel > 0) {
622 if (nRecHitEndcap > 0) {
628 typedef std::map<TrackPSimHitRef, RecHitIter> SimToRecHitMap;
629 SimToRecHitMap simToRecHitMap;
631 for (
const auto &
simHit : muonSimHits) {
632 const RPCDetId simDetId = static_cast<const RPCDetId>(
simHit->detUnitId());
636 const double simX =
simHit->localPosition().x();
638 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
639 const RPCDetId recDetId = static_cast<const RPCDetId>(recHitIter->rpcId());
640 const RPCRoll *recRoll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(recDetId));
644 if (simDetId != recDetId)
647 const double recX = recHitIter->localPosition().x();
648 const double newDx = fabs(recX - simX);
651 SimToRecHitMap::const_iterator prevSimToReco = simToRecHitMap.find(
simHit);
652 if (prevSimToReco == simToRecHitMap.end()) {
653 simToRecHitMap.insert(std::make_pair(
simHit, recHitIter));
655 const double oldDx = fabs(prevSimToReco->second->localPosition().x() - simX);
658 simToRecHitMap[
simHit] = recHitIter;
666 int nMatchHitBarrel = 0, nMatchHitEndcap = 0;
667 for (
const auto &
match : simToRecHitMap) {
669 RecHitIter recHitIter =
match.second;
671 const RPCDetId detId = static_cast<const RPCDetId>(
simHit->detUnitId());
672 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
674 const int region = roll->id().region();
675 const int ring = roll->id().ring();
677 const int station = roll->id().station();
681 const double simX =
simHit->localPosition().x();
682 const double recX = recHitIter->localPosition().x();
683 const double errX =
sqrt(recHitIter->localPositionError().xx());
684 const double dX = recX - simX;
685 const double pull =
errX == 0 ? -999 : dX /
errX;
723 for (reco::MuonCollection::const_iterator
muon = muonHandle->begin();
muon != muonHandle->end(); ++
muon) {
724 if (!
muon->isGlobalMuon())
727 int nRPCHitBarrel = 0;
728 int nRPCHitEndcap = 0;
732 if (!(*recHit)->isValid())
734 const DetId detId = (*recHit)->geographicalId();
737 const RPCDetId rpcDetId = static_cast<const RPCDetId>(detId);
739 if (rpcDetId.
region() == 0)
745 const int nRPCHit = nRPCHitBarrel + nRPCHitEndcap;
747 if (nRPCHitBarrel and nRPCHitEndcap) {
752 }
else if (nRPCHitBarrel) {
757 }
else if (nRPCHitEndcap) {
770 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
771 const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
772 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
774 const int region = roll->id().region();
775 const int ring = roll->id().ring();
777 const int station = roll->id().station();
782 for (
const auto &
match : simToRecHitMap) {
783 if (recHitIter ==
match.second) {
790 int nPunchMatched = 0;
792 for (
const auto &
simHit : pthrSimHits) {
793 const int absSimHitPType =
abs(
simHit->particleType());
794 if (absSimHitPType == 13)
797 const RPCDetId simDetId = static_cast<const RPCDetId>(
simHit->detUnitId());
798 if (simDetId == detId)
802 if (nPunchMatched > 0) {
816 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
817 const RPCDetId recDetId = static_cast<const RPCDetId>(recHitIter->rpcId());
818 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(recDetId));
820 const int region = roll->id().region();
827 const double recX = recHitIter->localPosition().x();
828 const double recErrX =
sqrt(recHitIter->localPositionError().xx());
831 for (SimHitIter simHitIter = simHitHandle->begin(); simHitIter != simHitHandle->end(); ++simHitIter) {
832 const RPCDetId simDetId = static_cast<const RPCDetId>(simHitIter->detUnitId());
833 const RPCRoll *simRoll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(simDetId));
837 if (simDetId != recDetId)
840 const double simX = simHitIter->localPosition().x();
841 const double dX = fabs(recX - simX);
843 if (dX / recErrX < 5) {