#include <ME0SegmentAlgorithm.h>

Inheritance diagram for ME0SegmentAlgorithm:

Public Types
typedef std::vector < HitAndPositionPtrContainer >	ProtoSegments
	Typedefs. More...

Public Types inherited from ME0SegmentAlgorithmBase
typedef std::vector < HitAndPosition >	HitAndPositionContainer

typedef std::vector< const HitAndPosition * >	HitAndPositionPtrContainer

Public Member Functions
	ME0SegmentAlgorithm (const edm::ParameterSet &ps)
	Constructor. More...

std::vector< ME0Segment >	run (const ME0Chamber *chamber, const HitAndPositionContainer &rechits) override

	~ME0SegmentAlgorithm () override
	Destructor. More...

Public Member Functions inherited from ME0SegmentAlgorithmBase
	ME0SegmentAlgorithmBase (const edm::ParameterSet &)
	Constructor. More...

virtual	~ME0SegmentAlgorithmBase ()
	Destructor. More...

Private Member Functions
void	buildSegments (const ME0Chamber *chamber, const HitAndPositionPtrContainer &rechits, std::vector< ME0Segment > &me0segs)

ProtoSegments	chainHits (const ME0Chamber *chamber, const HitAndPositionContainer &rechits)

ProtoSegments	clusterHits (const HitAndPositionContainer &rechits)
	Utility functions. More...

bool	isGoodToMerge (const ME0Chamber *chamber, const HitAndPositionPtrContainer &newChain, const HitAndPositionPtrContainer &oldChain)

Private Attributes
bool	debug

double	dEtaChainBoxMax

double	dPhiChainBoxMax

double	dTimeChainBoxMax

double	dXclusBoxMax

double	dYclusBoxMax

int	maxRecHitsInCluster

unsigned int	minHitsPerSegment

const std::string	myName

bool	preClustering

bool	preClustering_useChaining

std::unique_ptr< MuonSegFit >	sfit_

Static Private Attributes
static constexpr float	running_max = std::numeric_limits<float>::max()

Detailed Description

This algorithm is very basic no attemp to deal with ambiguities , noise etc. The ME0 track segments is built out of the rechit's in a the 6 ME0 Layer denoted as the ME0 Ensabmle .

Authors: Marcello Maggi

Definition at line 24 of file ME0SegmentAlgorithm.h.

Member Typedef Documentation

typedef std::vector<HitAndPositionPtrContainer> ME0SegmentAlgorithm::ProtoSegments

Typedefs.

Definition at line 27 of file ME0SegmentAlgorithm.h.

Constructor & Destructor Documentation

ME0SegmentAlgorithm::ME0SegmentAlgorithm ( const edm::ParameterSet & ps )

explicit

Constructor.

Definition at line 25 of file ME0SegmentAlgorithm.cc.

References debug, dEtaChainBoxMax, dPhiChainBoxMax, dTimeChainBoxMax, dXclusBoxMax, dYclusBoxMax, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), maxRecHitsInCluster, minHitsPerSegment, preClustering, and preClustering_useChaining.

     : ME0SegmentAlgorithmBase(ps), myName("ME0SegmentAlgorithm") {
   debug = ps.getUntrackedParameter<bool>("ME0Debug");
   minHitsPerSegment = ps.getParameter<unsigned int>("minHitsPerSegment");
   preClustering = ps.getParameter<bool>("preClustering");
   dXclusBoxMax = ps.getParameter<double>("dXclusBoxMax");
   dYclusBoxMax = ps.getParameter<double>("dYclusBoxMax");
   preClustering_useChaining = ps.getParameter<bool>("preClusteringUseChaining");
   dPhiChainBoxMax = ps.getParameter<double>("dPhiChainBoxMax");
   dEtaChainBoxMax = ps.getParameter<double>("dEtaChainBoxMax");
   dTimeChainBoxMax = ps.getParameter<double>("dTimeChainBoxMax");
   maxRecHitsInCluster = ps.getParameter<int>("maxRecHitsInCluster");
 
   edm::LogVerbatim("ME0SegmentAlgorithm")
       << "[ME0SegmentAlgorithm::ctor] Parameters to build segments :: "
       << "preClustering = " << preClustering << " preClustering_useChaining = " << preClustering_useChaining
       << " dPhiChainBoxMax = " << dPhiChainBoxMax << " dEtaChainBoxMax = " << dEtaChainBoxMax
       << " dTimeChainBoxMax = " << dTimeChainBoxMax << " minHitsPerSegment = " << minHitsPerSegment
       << " maxRecHitsInCluster = " << maxRecHitsInCluster;
 }

ME0SegmentAlgorithm::~ME0SegmentAlgorithm ( )

override

Destructor.

Definition at line 49 of file ME0SegmentAlgorithm.cc.

49 {}

Member Function Documentation

void ME0SegmentAlgorithm::buildSegments	(	const ME0Chamber *	chamber,
		const HitAndPositionPtrContainer &	rechits,
		std::vector< ME0Segment > &	me0segs
	)

private

Definition at line 314 of file ME0SegmentAlgorithm.cc.

References ME0RecHit::clone(), ME0Chamber::computeDeltaPhi(), minHitsPerSegment, funct::pow(), ME0RecHit::setPosition(), sfit_, mathSSE::sqrt(), and createJobs::tmp.

Referenced by run().

                                                                         {
   if (rechits.size() < minHitsPerSegment)
     return;
 
 #ifdef EDM_ML_DEBUG  // have lines below only compiled when in debug mode
   edm::LogVerbatim("ME0SegmentAlgorithm")
       << "[ME0SegmentAlgorithm::buildSegments] will now try to fit a ME0Segment from collection of " << rechits.size()
       << " ME0 RecHits";
   for (auto rh = rechits.begin(); rh != rechits.end(); ++rh) {
     auto me0id = (*rh)->rh->me0Id();
     auto rhLP = (*rh)->lp;
     edm::LogVerbatim("ME0SegmentAlgorithm")
         << "[RecHit :: Loc x = " << std::showpos << std::setw(9) << rhLP.x() << " Loc y = " << std::showpos
         << std::setw(9) << rhLP.y() << " Time = " << std::showpos << (*rh)->rh->tof() << " -- " << me0id.rawId()
         << " = " << me0id << " ]";
   }
 #endif
 
   MuonSegFit::MuonRecHitContainer muonRecHits;
   std::vector<const ME0RecHit*> bareRHs;
 
   // select hits from the ensemble and sort it
   for (auto rh = rechits.begin(); rh != rechits.end(); rh++) {
     bareRHs.push_back((*rh)->rh);
     // for segFit - using local point in chamber frame
     ME0RecHit* newRH = (*rh)->rh->clone();
     newRH->setPosition((*rh)->lp);
 
     MuonSegFit::MuonRecHitPtr trkRecHit(newRH);
     muonRecHits.push_back(trkRecHit);
   }
 
   // The actual fit on all hits of the vector of the selected Tracking RecHits:
   sfit_ = std::make_unique<MuonSegFit>(muonRecHits);
   bool goodfit = sfit_->fit();
   edm::LogVerbatim("ME0SegmentAlgorithm") << "[ME0SegmentAlgorithm::buildSegments] ME0Segment fit done";
 
   // quit function if fit was not OK
   if (!goodfit) {
     for (auto rh : muonRecHits)
       rh.reset();
     return;
   }
 
   // obtain all information necessary to make the segment:
   LocalPoint protoIntercept = sfit_->intercept();
   LocalVector protoDirection = sfit_->localdir();
   AlgebraicSymMatrix protoErrors = sfit_->covarianceMatrix();
   double protoChi2 = sfit_->chi2();
   // Calculate the central value and uncertainty of the segment time
   float averageTime = 0.;
   for (auto rh = rechits.begin(); rh != rechits.end(); ++rh) {
     averageTime += (*rh)->rh->tof();
   }
   if (!rechits.empty())
     averageTime = averageTime / (rechits.size());
   float timeUncrt = 0.;
   for (auto rh = rechits.begin(); rh != rechits.end(); ++rh) {
     timeUncrt += pow((*rh)->rh->tof() - averageTime, 2);
   }
 
   if (rechits.size() > 1)
     timeUncrt = timeUncrt / (rechits.size() - 1);
   timeUncrt = sqrt(timeUncrt);
 
   const float dPhi = chamber->computeDeltaPhi(protoIntercept, protoDirection);
 
   // save all information inside GEMCSCSegment
   edm::LogVerbatim("ME0SegmentAlgorithm")
       << "[ME0SegmentAlgorithm::buildSegments] will now try to make ME0Segment from collection of " << rechits.size()
       << " ME0 RecHits";
   ME0Segment tmp(bareRHs, protoIntercept, protoDirection, protoErrors, protoChi2, averageTime, timeUncrt, dPhi);
 
   edm::LogVerbatim("ME0SegmentAlgorithm") << "[ME0SegmentAlgorithm::buildSegments] ME0Segment made";
   edm::LogVerbatim("ME0SegmentAlgorithm") << "[ME0SegmentAlgorithm::buildSegments] " << tmp;
 
   for (auto rh : muonRecHits)
     rh.reset();
   me0segs.push_back(tmp);
   return;
 }

ME0SegmentAlgorithm::ProtoSegments ME0SegmentAlgorithm::chainHits	(	const ME0Chamber *	chamber,
		const HitAndPositionContainer &	rechits
	)

private

Definition at line 214 of file ME0SegmentAlgorithm.cc.

References funct::abs(), SplitLinear::begin, dTimeChainBoxMax, dataset::end, mps_fire::i, isGoodToMerge(), and DetachedQuadStep_cff::seeds.

Referenced by run().

                                                                                                           {
   ProtoSegments rechits_chains;
   ProtoSegments seeds;
   seeds.reserve(rechits.size());
   std::vector<bool> usedCluster(rechits.size(), false);
 
   // split rechits into subvectors and return vector of vectors:
   // Loop over rechits
   // Create one seed per hit
   for (unsigned int i = 0; i < rechits.size(); ++i) {
     if (std::abs(rechits[i].rh->tof()) > dTimeChainBoxMax)
       continue;
     seeds.push_back(HitAndPositionPtrContainer(1, &rechits[i]));
   }
 
   // merge chains that are too close ("touch" each other)
   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
     for (size_t MMM = NNN + 1; MMM < seeds.size(); ++MMM) {
       if (usedCluster[MMM] || usedCluster[NNN]) {
         continue;
       }
       // all is in the way we define "good";
       // try not to "cluster" the hits but to "chain" them;
       // it does the clustering but also does a better job
       // for inclined tracks (not clustering them together;
       // crossed tracks would be still clustered together)
       // 22.12.09: In fact it is not much more different
       // than the "clustering", we just introduce another
       // variable in the game - Z. And it makes sense
       // to re-introduce Y (or actually wire group mumber)
       // in a similar way as for the strip number - see
       // the code below.
       bool goodToMerge = isGoodToMerge(chamber, seeds[NNN], seeds[MMM]);
       if (goodToMerge) {
         // merge chains!
         // merge by adding seed NNN to seed MMM and erasing seed NNN
 
         // add seed NNN to MMM (lower to larger number)
         seeds[MMM].insert(seeds[MMM].end(), seeds[NNN].begin(), seeds[NNN].end());
 
         // mark seed NNN as used
         usedCluster[NNN] = true;
         // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to
         // next seed (NNN+1)
         break;
       }
     }
   }
 
   // hand over the final seeds to the output
   // would be more elegant if we could do the above step with
   // erasing the merged ones, rather than the
 
   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
     if (usedCluster[NNN])
       continue;  //skip seeds that have been marked as used up in merging
     rechits_chains.push_back(seeds[NNN]);
   }
 
   //***************************************************************
 
   return rechits_chains;
 }

ME0SegmentAlgorithm::ProtoSegments ME0SegmentAlgorithm::clusterHits ( const HitAndPositionContainer & rechits )

private

Utility functions.

Definition at line 102 of file ME0SegmentAlgorithm.cc.

References SplitLinear::begin, dXclusBoxMax, dYclusBoxMax, dataset::end, mps_fire::i, LogDebug, running_max, DetachedQuadStep_cff::seeds, and findQualityFiles::size.

Referenced by run().

                                                                                                         {
   ProtoSegments rechits_clusters;  // this is a collection of groups of rechits
 
   float dXclus_box = 0.0;
   float dYclus_box = 0.0;
 
   ProtoSegments seeds;
   seeds.reserve(rechits.size());
 
   std::vector<float> running_meanX;
   running_meanX.reserve(rechits.size());
   std::vector<float> running_meanY;
   running_meanY.reserve(rechits.size());
 
   std::vector<float> seed_minX;
   seed_minX.reserve(rechits.size());
   std::vector<float> seed_maxX;
   seed_maxX.reserve(rechits.size());
   std::vector<float> seed_minY;
   seed_minY.reserve(rechits.size());
   std::vector<float> seed_maxY;
   seed_maxY.reserve(rechits.size());
 
   // split rechits into subvectors and return vector of vectors:
   // Loop over rechits
   // Create one seed per hit
   for (unsigned int i = 0; i < rechits.size(); ++i) {
     seeds.push_back(HitAndPositionPtrContainer(1, &rechits[i]));
 
     // First added hit in seed defines the mean to which the next hit is compared
     // for this seed.
 
     running_meanX.push_back(rechits[i].lp.x());
     running_meanY.push_back(rechits[i].lp.y());
 
     // set min/max X and Y for box containing the hits in the precluster:
     seed_minX.push_back(rechits[i].lp.x());
     seed_maxX.push_back(rechits[i].lp.x());
     seed_minY.push_back(rechits[i].lp.y());
     seed_maxY.push_back(rechits[i].lp.y());
   }
 
   // merge clusters that are too close
   // measure distance between final "running mean"
   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
     for (size_t MMM = NNN + 1; MMM < seeds.size(); ++MMM) {
       if (running_meanX[MMM] == running_max || running_meanX[NNN] == running_max) {
         LogDebug("ME0SegmentAlgorithm") << "[ME0SegmentAlgorithm::clusterHits]: ALARM! Skipping used seeds, this "
                                            "should not happen - inform developers!";
         continue;  //skip seeds that have been used
       }
 
       // calculate cut criteria for simple running mean distance cut:
       //dXclus = fabs(running_meanX[NNN] - running_meanX[MMM]);
       //dYclus = fabs(running_meanY[NNN] - running_meanY[MMM]);
       // calculate minmal distance between precluster boxes containing the hits:
       if (running_meanX[NNN] > running_meanX[MMM])
         dXclus_box = seed_minX[NNN] - seed_maxX[MMM];
       else
         dXclus_box = seed_minX[MMM] - seed_maxX[NNN];
       if (running_meanY[NNN] > running_meanY[MMM])
         dYclus_box = seed_minY[NNN] - seed_maxY[MMM];
       else
         dYclus_box = seed_minY[MMM] - seed_maxY[NNN];
 
       if (dXclus_box < dXclusBoxMax && dYclus_box < dYclusBoxMax) {
         // merge clusters!
         // merge by adding seed NNN to seed MMM and erasing seed NNN
 
         // calculate running mean for the merged seed:
         if (seeds[NNN].size() + seeds[MMM].size() != 0) {
           running_meanX[MMM] = (running_meanX[NNN] * seeds[NNN].size() + running_meanX[MMM] * seeds[MMM].size()) /
                                (seeds[NNN].size() + seeds[MMM].size());
           running_meanY[MMM] = (running_meanY[NNN] * seeds[NNN].size() + running_meanY[MMM] * seeds[MMM].size()) /
                                (seeds[NNN].size() + seeds[MMM].size());
         }
 
         // update min/max X and Y for box containing the hits in the merged cluster:
         if (seed_minX[NNN] < seed_minX[MMM])
           seed_minX[MMM] = seed_minX[NNN];
         if (seed_maxX[NNN] > seed_maxX[MMM])
           seed_maxX[MMM] = seed_maxX[NNN];
         if (seed_minY[NNN] < seed_minY[MMM])
           seed_minY[MMM] = seed_minY[NNN];
         if (seed_maxY[NNN] > seed_maxY[MMM])
           seed_maxY[MMM] = seed_maxY[NNN];
 
         // add seed NNN to MMM (lower to larger number)
         seeds[MMM].insert(seeds[MMM].end(), seeds[NNN].begin(), seeds[NNN].end());
 
         // mark seed NNN as used (at the moment just set running mean to 999999.)
         running_meanX[NNN] = running_max;
         running_meanY[NNN] = running_max;
         // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to
         // next seed (NNN+1)
         break;
       }
     }
   }
 
   // hand over the final seeds to the output
   // would be more elegant if we could do the above step with
   // erasing the merged ones, rather than the
   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
     if (running_meanX[NNN] == running_max)
       continue;  //skip seeds that have been marked as used up in merging
     rechits_clusters.push_back(seeds[NNN]);
   }
 
   return rechits_clusters;
 }

bool ME0SegmentAlgorithm::isGoodToMerge	(	const ME0Chamber *	chamber,
		const HitAndPositionPtrContainer &	newChain,
		const HitAndPositionPtrContainer &	oldChain
	)

private

Definition at line 279 of file ME0SegmentAlgorithm.cc.

References funct::abs(), reco::deltaPhi(), dEtaChainBoxMax, dPhiChainBoxMax, dTimeChainBoxMax, PV3DBase< T, PVType, FrameType >::eta(), ME0Chamber::id(), phase1PixelTopology::layer, ME0DetId::nlayers(), and PV3DBase< T, PVType, FrameType >::phi().

Referenced by chainHits().

                                                                                     {
   for (size_t iRH_new = 0; iRH_new < newChain.size(); ++iRH_new) {
     GlobalPoint pos_new = newChain[iRH_new]->gp;
 
     for (size_t iRH_old = 0; iRH_old < oldChain.size(); ++iRH_old) {
       GlobalPoint pos_old = oldChain[iRH_old]->gp;
       // to be chained, two hits need to be in neighbouring layers...
       // or better allow few missing layers (upto 3 to avoid inefficiencies);
       // however we'll not make an angle correction because it
       // worsen the situation in some of the "regular" cases
       // (not making the correction means that the conditions for
       // forming a cluster are different if we have missing layers -
       // this could affect events at the boundaries )
 
       // to be chained, two hits need also to be "close" in phi and eta
       if (std::abs(reco::deltaPhi(float(pos_new.phi()), float(pos_old.phi()))) >= dPhiChainBoxMax)
         continue;
       if (std::abs(pos_new.eta() - pos_old.eta()) >= dEtaChainBoxMax)
         continue;
       // and the difference in layer index should be < (nlayers-1)
       if (std::abs(int(newChain[iRH_new]->layer) - int(oldChain[iRH_old]->layer)) >= (chamber->id().nlayers() - 1))
         continue;
       // and they should have a time difference compatible with the hypothesis
       // that the rechits originate from the same particle, but were detected in different layers
       if (std::abs(newChain[iRH_new]->rh->tof() - oldChain[iRH_old]->rh->tof()) >= dTimeChainBoxMax)
         continue;
 
       return true;
     }
   }
   return false;
 }

std::vector< ME0Segment > ME0SegmentAlgorithm::run	(	const ME0Chamber *	chamber,
		const HitAndPositionContainer &	rechits
	)

overridevirtual

Build segments for all desired groups of hits

Implements ME0SegmentAlgorithmBase.

Definition at line 51 of file ME0SegmentAlgorithm.cc.

References buildSegments(), chainHits(), clusterHits(), ME0Chamber::id(), preClustering, and preClustering_useChaining.

                                                                                                                 {
 #ifdef EDM_ML_DEBUG  // have lines below only compiled when in debug mode
   ME0DetId chId(chamber->id());
   edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegmentAlgorithm::run] build segments in chamber " << chId
                                    << " which contains " << rechits.size() << " rechits";
   for (auto rh = rechits.begin(); rh != rechits.end(); ++rh) {
     auto me0id = rh->rh->me0Id();
     auto rhLP = rh->lp;
     edm::LogVerbatim("ME0SegmentAlgorithm")
         << "[RecHit :: Loc x = " << std::showpos << std::setw(9) << rhLP.x() << " Loc y = " << std::showpos
         << std::setw(9) << rhLP.y() << " Time = " << std::showpos << rh->rh->tof() << " -- " << me0id.rawId() << " = "
         << me0id << " ]";
   }
 #endif
 
   // pre-cluster rechits and loop over all sub clusters separately
   std::vector<ME0Segment> segments_temp;
   std::vector<ME0Segment> segments;
   ProtoSegments rechits_clusters;  // this is a collection of groups of rechits
 
   if (preClustering) {
     // run a pre-clusterer on the given rechits to split obviously separated segment seeds:
     if (preClustering_useChaining) {
       // it uses X,Y,Z information; there are no configurable parameters used;
       // the X, Y, Z "cuts" are just (much) wider than reasonable high pt segments
       rechits_clusters = this->chainHits(chamber, rechits);
     } else {
       // it uses X,Y information + configurable parameters
       rechits_clusters = this->clusterHits(rechits);
     }
     // loop over the found clusters:
     for (auto sub_rechits = rechits_clusters.begin(); sub_rechits != rechits_clusters.end(); ++sub_rechits) {
       // clear the buffer for the subset of segments:
       segments_temp.clear();
       // build the subset of segments:
       this->buildSegments(chamber, (*sub_rechits), segments_temp);
       // add the found subset of segments to the collection of all segments in this chamber:
       segments.insert(segments.end(), segments_temp.begin(), segments_temp.end());
     }
     return segments;
   } else {
     HitAndPositionPtrContainer ptrRH;
     ptrRH.reserve(rechits.size());
     for (const auto& rh : rechits)
       ptrRH.push_back(&rh);
     this->buildSegments(chamber, ptrRH, segments);
     return segments;
   }
 }