GEMSegment/plugins/ME0SegmentAlgorithm.cc

0001 /**
0002  * \file ME0SegmentAlgorithm.cc
0003  *  based on CSCSegAlgoST.cc
0004  *
0005  *  \authors: Marcello Maggi, Jason Lee
0006  */
0007
0008 #include "ME0SegmentAlgorithm.h"
0009 #include "MuonSegFit.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0014 #include "DataFormats/Math/interface/deltaPhi.h"
0015 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0016
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <iostream>
0020 #include <string>
0021
0022 /* Constructor
0023  *
0024  */
0025 ME0SegmentAlgorithm::ME0SegmentAlgorithm(const edm::ParameterSet& ps)
0026     : ME0SegmentAlgorithmBase(ps), myName("ME0SegmentAlgorithm") {
0027   debug = ps.getUntrackedParameter<bool>("ME0Debug");
0028   minHitsPerSegment = ps.getParameter<unsigned int>("minHitsPerSegment");
0029   preClustering = ps.getParameter<bool>("preClustering");
0030   dXclusBoxMax = ps.getParameter<double>("dXclusBoxMax");
0031   dYclusBoxMax = ps.getParameter<double>("dYclusBoxMax");
0032   preClustering_useChaining = ps.getParameter<bool>("preClusteringUseChaining");
0033   dPhiChainBoxMax = ps.getParameter<double>("dPhiChainBoxMax");
0034   dEtaChainBoxMax = ps.getParameter<double>("dEtaChainBoxMax");
0035   dTimeChainBoxMax = ps.getParameter<double>("dTimeChainBoxMax");
0036   maxRecHitsInCluster = ps.getParameter<int>("maxRecHitsInCluster");
0037
0038   edm::LogVerbatim("ME0SegmentAlgorithm")
0039       << "[ME0SegmentAlgorithm::ctor] Parameters to build segments :: "
0040       << "preClustering = " << preClustering << " preClustering_useChaining = " << preClustering_useChaining
0041       << " dPhiChainBoxMax = " << dPhiChainBoxMax << " dEtaChainBoxMax = " << dEtaChainBoxMax
0042       << " dTimeChainBoxMax = " << dTimeChainBoxMax << " minHitsPerSegment = " << minHitsPerSegment
0043       << " maxRecHitsInCluster = " << maxRecHitsInCluster;
0044 }
0045
0046 /* Destructor
0047  *
0048  */
0049 ME0SegmentAlgorithm::~ME0SegmentAlgorithm() {}
0050
0051 std::vector<ME0Segment> ME0SegmentAlgorithm::run(const ME0Chamber* chamber, const HitAndPositionContainer& rechits) {
0052 #ifdef EDM_ML_DEBUG  // have lines below only compiled when in debug mode
0053   ME0DetId chId(chamber->id());
0054   edm::LogVerbatim("ME0SegAlgoMM") << "[ME0SegmentAlgorithm::run] build segments in chamber " << chId
0055                                    << " which contains " << rechits.size() << " rechits";
0056   for (auto rh = rechits.begin(); rh != rechits.end(); ++rh) {
0057     auto me0id = rh->rh->me0Id();
0058     auto rhLP = rh->lp;
0059     edm::LogVerbatim("ME0SegmentAlgorithm")
0060         << "[RecHit :: Loc x = " << std::showpos << std::setw(9) << rhLP.x() << " Loc y = " << std::showpos
0061         << std::setw(9) << rhLP.y() << " Time = " << std::showpos << rh->rh->tof() << " -- " << me0id.rawId() << " = "
0062         << me0id << " ]";
0063   }
0064 #endif
0065
0066   // pre-cluster rechits and loop over all sub clusters separately
0067   std::vector<ME0Segment> segments_temp;
0068   std::vector<ME0Segment> segments;
0069   ProtoSegments rechits_clusters;  // this is a collection of groups of rechits
0070
0071   if (preClustering) {
0072     // run a pre-clusterer on the given rechits to split obviously separated segment seeds:
0073     if (preClustering_useChaining) {
0074       // it uses X,Y,Z information; there are no configurable parameters used;
0075       // the X, Y, Z "cuts" are just (much) wider than reasonable high pt segments
0076       rechits_clusters = this->chainHits(chamber, rechits);
0077     } else {
0078       // it uses X,Y information + configurable parameters
0079       rechits_clusters = this->clusterHits(rechits);
0080     }
0081     // loop over the found clusters:
0082     for (auto sub_rechits = rechits_clusters.begin(); sub_rechits != rechits_clusters.end(); ++sub_rechits) {
0083       // clear the buffer for the subset of segments:
0084       segments_temp.clear();
0085       // build the subset of segments:
0086       this->buildSegments(chamber, (*sub_rechits), segments_temp);
0087       // add the found subset of segments to the collection of all segments in this chamber:
0088       segments.insert(segments.end(), segments_temp.begin(), segments_temp.end());
0089     }
0090     return segments;
0091   } else {
0092     HitAndPositionPtrContainer ptrRH;
0093     ptrRH.reserve(rechits.size());
0094     for (const auto& rh : rechits)
0095       ptrRH.push_back(&rh);
0096     this->buildSegments(chamber, ptrRH, segments);
0097     return segments;
0098   }
0099 }
0100
0101 // ********************************************************************;
0102 ME0SegmentAlgorithm::ProtoSegments ME0SegmentAlgorithm::clusterHits(const HitAndPositionContainer& rechits) {
0103   ProtoSegments rechits_clusters;  // this is a collection of groups of rechits
0104
0105   float dXclus_box = 0.0;
0106   float dYclus_box = 0.0;
0107
0108   ProtoSegments seeds;
0109   seeds.reserve(rechits.size());
0110
0111   std::vector<float> running_meanX;
0112   running_meanX.reserve(rechits.size());
0113   std::vector<float> running_meanY;
0114   running_meanY.reserve(rechits.size());
0115
0116   std::vector<float> seed_minX;
0117   seed_minX.reserve(rechits.size());
0118   std::vector<float> seed_maxX;
0119   seed_maxX.reserve(rechits.size());
0120   std::vector<float> seed_minY;
0121   seed_minY.reserve(rechits.size());
0122   std::vector<float> seed_maxY;
0123   seed_maxY.reserve(rechits.size());
0124
0125   // split rechits into subvectors and return vector of vectors:
0126   // Loop over rechits
0127   // Create one seed per hit
0128   for (unsigned int i = 0; i < rechits.size(); ++i) {
0129     seeds.push_back(HitAndPositionPtrContainer(1, &rechits[i]));
0130
0131     // First added hit in seed defines the mean to which the next hit is compared
0132     // for this seed.
0133
0134     running_meanX.push_back(rechits[i].lp.x());
0135     running_meanY.push_back(rechits[i].lp.y());
0136
0137     // set min/max X and Y for box containing the hits in the precluster:
0138     seed_minX.push_back(rechits[i].lp.x());
0139     seed_maxX.push_back(rechits[i].lp.x());
0140     seed_minY.push_back(rechits[i].lp.y());
0141     seed_maxY.push_back(rechits[i].lp.y());
0142   }
0143
0144   // merge clusters that are too close
0145   // measure distance between final "running mean"
0146   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0147     for (size_t MMM = NNN + 1; MMM < seeds.size(); ++MMM) {
0148       if (running_meanX[MMM] == running_max || running_meanX[NNN] == running_max) {
0149         LogDebug("ME0SegmentAlgorithm") << "[ME0SegmentAlgorithm::clusterHits]: ALARM! Skipping used seeds, this "
0150                                            "should not happen - inform developers!";
0151         continue;  //skip seeds that have been used
0152       }
0153
0154       // calculate cut criteria for simple running mean distance cut:
0155       //dXclus = fabs(running_meanX[NNN] - running_meanX[MMM]);
0156       //dYclus = fabs(running_meanY[NNN] - running_meanY[MMM]);
0157       // calculate minmal distance between precluster boxes containing the hits:
0158       if (running_meanX[NNN] > running_meanX[MMM])
0159         dXclus_box = seed_minX[NNN] - seed_maxX[MMM];
0160       else
0161         dXclus_box = seed_minX[MMM] - seed_maxX[NNN];
0162       if (running_meanY[NNN] > running_meanY[MMM])
0163         dYclus_box = seed_minY[NNN] - seed_maxY[MMM];
0164       else
0165         dYclus_box = seed_minY[MMM] - seed_maxY[NNN];
0166
0167       if (dXclus_box < dXclusBoxMax && dYclus_box < dYclusBoxMax) {
0168         // merge clusters!
0169         // merge by adding seed NNN to seed MMM and erasing seed NNN
0170
0171         // calculate running mean for the merged seed:
0172         if (seeds[NNN].size() + seeds[MMM].size() != 0) {
0173           running_meanX[MMM] = (running_meanX[NNN] * seeds[NNN].size() + running_meanX[MMM] * seeds[MMM].size()) /
0174                                (seeds[NNN].size() + seeds[MMM].size());
0175           running_meanY[MMM] = (running_meanY[NNN] * seeds[NNN].size() + running_meanY[MMM] * seeds[MMM].size()) /
0176                                (seeds[NNN].size() + seeds[MMM].size());
0177         }
0178
0179         // update min/max X and Y for box containing the hits in the merged cluster:
0180         if (seed_minX[NNN] < seed_minX[MMM])
0181           seed_minX[MMM] = seed_minX[NNN];
0182         if (seed_maxX[NNN] > seed_maxX[MMM])
0183           seed_maxX[MMM] = seed_maxX[NNN];
0184         if (seed_minY[NNN] < seed_minY[MMM])
0185           seed_minY[MMM] = seed_minY[NNN];
0186         if (seed_maxY[NNN] > seed_maxY[MMM])
0187           seed_maxY[MMM] = seed_maxY[NNN];
0188
0189         // add seed NNN to MMM (lower to larger number)
0190         seeds[MMM].insert(seeds[MMM].end(), seeds[NNN].begin(), seeds[NNN].end());
0191
0192         // mark seed NNN as used (at the moment just set running mean to 999999.)
0193         running_meanX[NNN] = running_max;
0194         running_meanY[NNN] = running_max;
0195         // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to
0196         // next seed (NNN+1)
0197         break;
0198       }
0199     }
0200   }
0201
0202   // hand over the final seeds to the output
0203   // would be more elegant if we could do the above step with
0204   // erasing the merged ones, rather than the
0205   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0206     if (running_meanX[NNN] == running_max)
0207       continue;  //skip seeds that have been marked as used up in merging
0208     rechits_clusters.push_back(seeds[NNN]);
0209   }
0210
0211   return rechits_clusters;
0212 }
0213
0214 ME0SegmentAlgorithm::ProtoSegments ME0SegmentAlgorithm::chainHits(const ME0Chamber* chamber,
0215                                                                   const HitAndPositionContainer& rechits) {
0216   ProtoSegments rechits_chains;
0217   ProtoSegments seeds;
0218   seeds.reserve(rechits.size());
0219   std::vector<bool> usedCluster(rechits.size(), false);
0220
0221   // split rechits into subvectors and return vector of vectors:
0222   // Loop over rechits
0223   // Create one seed per hit
0224   for (unsigned int i = 0; i < rechits.size(); ++i) {
0225     if (std::abs(rechits[i].rh->tof()) > dTimeChainBoxMax)
0226       continue;
0227     seeds.push_back(HitAndPositionPtrContainer(1, &rechits[i]));
0228   }
0229
0230   // merge chains that are too close ("touch" each other)
0231   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0232     for (size_t MMM = NNN + 1; MMM < seeds.size(); ++MMM) {
0233       if (usedCluster[MMM] || usedCluster[NNN]) {
0234         continue;
0235       }
0236       // all is in the way we define "good";
0237       // try not to "cluster" the hits but to "chain" them;
0238       // it does the clustering but also does a better job
0239       // for inclined tracks (not clustering them together;
0240       // crossed tracks would be still clustered together)
0241       // 22.12.09: In fact it is not much more different
0242       // than the "clustering", we just introduce another
0243       // variable in the game - Z. And it makes sense
0244       // to re-introduce Y (or actually wire group mumber)
0245       // in a similar way as for the strip number - see
0246       // the code below.
0247       bool goodToMerge = isGoodToMerge(chamber, seeds[NNN], seeds[MMM]);
0248       if (goodToMerge) {
0249         // merge chains!
0250         // merge by adding seed NNN to seed MMM and erasing seed NNN
0251
0252         // add seed NNN to MMM (lower to larger number)
0253         seeds[MMM].insert(seeds[MMM].end(), seeds[NNN].begin(), seeds[NNN].end());
0254
0255         // mark seed NNN as used
0256         usedCluster[NNN] = true;
0257         // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to
0258         // next seed (NNN+1)
0259         break;
0260       }
0261     }
0262   }
0263
0264   // hand over the final seeds to the output
0265   // would be more elegant if we could do the above step with
0266   // erasing the merged ones, rather than the
0267
0268   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0269     if (usedCluster[NNN])
0270       continue;  //skip seeds that have been marked as used up in merging
0271     rechits_chains.push_back(seeds[NNN]);
0272   }
0273
0274   //***************************************************************
0275
0276   return rechits_chains;
0277 }
0278
0279 bool ME0SegmentAlgorithm::isGoodToMerge(const ME0Chamber* chamber,
0280                                         const HitAndPositionPtrContainer& newChain,
0281                                         const HitAndPositionPtrContainer& oldChain) {
0282   for (size_t iRH_new = 0; iRH_new < newChain.size(); ++iRH_new) {
0283     GlobalPoint pos_new = newChain[iRH_new]->gp;
0284
0285     for (size_t iRH_old = 0; iRH_old < oldChain.size(); ++iRH_old) {
0286       GlobalPoint pos_old = oldChain[iRH_old]->gp;
0287       // to be chained, two hits need to be in neighbouring layers...
0288       // or better allow few missing layers (upto 3 to avoid inefficiencies);
0289       // however we'll not make an angle correction because it
0290       // worsen the situation in some of the "regular" cases
0291       // (not making the correction means that the conditions for
0292       // forming a cluster are different if we have missing layers -
0293       // this could affect events at the boundaries )
0294
0295       // to be chained, two hits need also to be "close" in phi and eta
0296       if (std::abs(reco::deltaPhi(float(pos_new.phi()), float(pos_old.phi()))) >= dPhiChainBoxMax)
0297         continue;
0298       if (std::abs(pos_new.eta() - pos_old.eta()) >= dEtaChainBoxMax)
0299         continue;
0300       // and the difference in layer index should be < (nlayers-1)
0301       if (std::abs(int(newChain[iRH_new]->layer) - int(oldChain[iRH_old]->layer)) >= (chamber->id().nlayers() - 1))
0302         continue;
0303       // and they should have a time difference compatible with the hypothesis
0304       // that the rechits originate from the same particle, but were detected in different layers
0305       if (std::abs(newChain[iRH_new]->rh->tof() - oldChain[iRH_old]->rh->tof()) >= dTimeChainBoxMax)
0306         continue;
0307
0308       return true;
0309     }
0310   }
0311   return false;
0312 }
0313
0314 void ME0SegmentAlgorithm::buildSegments(const ME0Chamber* chamber,
0315                                         const HitAndPositionPtrContainer& rechits,
0316                                         std::vector<ME0Segment>& me0segs) {
0317   if (rechits.size() < minHitsPerSegment)
0318     return;
0319
0320 #ifdef EDM_ML_DEBUG  // have lines below only compiled when in debug mode
0321   edm::LogVerbatim("ME0SegmentAlgorithm")
0322       << "[ME0SegmentAlgorithm::buildSegments] will now try to fit a ME0Segment from collection of " << rechits.size()
0323       << " ME0 RecHits";
0324   for (auto rh = rechits.begin(); rh != rechits.end(); ++rh) {
0325     auto me0id = (*rh)->rh->me0Id();
0326     auto rhLP = (*rh)->lp;
0327     edm::LogVerbatim("ME0SegmentAlgorithm")
0328         << "[RecHit :: Loc x = " << std::showpos << std::setw(9) << rhLP.x() << " Loc y = " << std::showpos
0329         << std::setw(9) << rhLP.y() << " Time = " << std::showpos << (*rh)->rh->tof() << " -- " << me0id.rawId()
0330         << " = " << me0id << " ]";
0331   }
0332 #endif
0333
0334   MuonSegFit::MuonRecHitContainer muonRecHits;
0335   std::vector<const ME0RecHit*> bareRHs;
0336
0337   // select hits from the ensemble and sort it
0338   for (auto rh = rechits.begin(); rh != rechits.end(); rh++) {
0339     bareRHs.push_back((*rh)->rh);
0340     // for segFit - using local point in chamber frame
0341     ME0RecHit* newRH = (*rh)->rh->clone();
0342     newRH->setPosition((*rh)->lp);
0343
0344     MuonSegFit::MuonRecHitPtr trkRecHit(newRH);
0345     muonRecHits.push_back(trkRecHit);
0346   }
0347
0348   // The actual fit on all hits of the vector of the selected Tracking RecHits:
0349   sfit_ = std::make_unique<MuonSegFit>(muonRecHits);
0350   bool goodfit = sfit_->fit();
0351   edm::LogVerbatim("ME0SegmentAlgorithm") << "[ME0SegmentAlgorithm::buildSegments] ME0Segment fit done";
0352
0353   // quit function if fit was not OK
0354   if (!goodfit) {
0355     for (auto rh : muonRecHits)
0356       rh.reset();
0357     return;
0358   }
0359
0360   // obtain all information necessary to make the segment:
0361   LocalPoint protoIntercept = sfit_->intercept();
0362   LocalVector protoDirection = sfit_->localdir();
0363   AlgebraicSymMatrix protoErrors = sfit_->covarianceMatrix();
0364   double protoChi2 = sfit_->chi2();
0365   // Calculate the central value and uncertainty of the segment time
0366   float averageTime = 0.;
0367   for (auto rh = rechits.begin(); rh != rechits.end(); ++rh) {
0368     averageTime += (*rh)->rh->tof();
0369   }
0370   if (!rechits.empty())
0371     averageTime = averageTime / (rechits.size());
0372   float timeUncrt = 0.;
0373   for (auto rh = rechits.begin(); rh != rechits.end(); ++rh) {
0374     timeUncrt += pow((*rh)->rh->tof() - averageTime, 2);
0375   }
0376
0377   if (rechits.size() > 1)
0378     timeUncrt = timeUncrt / (rechits.size() - 1);
0379   timeUncrt = sqrt(timeUncrt);
0380
0381   const float dPhi = chamber->computeDeltaPhi(protoIntercept, protoDirection);
0382
0383   // save all information inside GEMCSCSegment
0384   edm::LogVerbatim("ME0SegmentAlgorithm")
0385       << "[ME0SegmentAlgorithm::buildSegments] will now try to make ME0Segment from collection of " << rechits.size()
0386       << " ME0 RecHits";
0387   ME0Segment tmp(bareRHs, protoIntercept, protoDirection, protoErrors, protoChi2, averageTime, timeUncrt, dPhi);
0388
0389   edm::LogVerbatim("ME0SegmentAlgorithm") << "[ME0SegmentAlgorithm::buildSegments] ME0Segment made";
0390   edm::LogVerbatim("ME0SegmentAlgorithm") << "[ME0SegmentAlgorithm::buildSegments] " << tmp;
0391
0392   for (auto rh : muonRecHits)
0393     rh.reset();
0394   me0segs.push_back(tmp);
0395   return;
0396 }